# Fault Diagnosis for a Multistage Planetary Gear Set Using Model-Based Simulation and Experimental Investigation

Authors: Guoyan Li; Fangyi Li; Yifan Wang; Dehao Dong
Journal: Shock and Vibration
 (2016)

---

## Abstract

The gear damage will induce modulation effects in vibration signals. A thorough analysis of modulation sidebands spectral structure is necessary for fault diagnosis of planetary gear set. However, the spectral characteristics are complicated in practice, especially for a multistage planetary gear set which contains close frequency components. In this study, a coupled lateral and torsional dynamic model is established to predict the modulation sidebands of a two-stage compound planetary gear set. An improved potential energy method is used to calculate the time-varying mesh stiffness of each gear pair, and the influence of crack propagation on the mesh stiffness is analyzed. The simulated signals of the gear set are obtained by using Runge-Kutta numerical analysis method. Meanwhile, the sidebands characteristics are summarized to exhibit the modulation effects caused by sun gear damage. At the end, the experimental signals collected from an industrial SD16 planetary gearbox are analyzed to verify the theoretical derivations. The results of experiment agree well with the simulated analysis.


---

## Body

## 1. Introduction

Remanufacturing is one of the best solutions to avoid resource shortage in the construction machinery industry. However, the quality of remanufactured products has long been an issue of concern. Planetary gearbox is widely used in large-scale and complex mechanical equipment such as wind turbines, helicopters, and construction machinery. The remanufacturing of planetary gearbox will create huge value for the economy development. The fault diagnosis of recycled gearboxes plays a vital role in guaranteeing the quality of remanufactured products. However, the recycled objects have experienced a service cycle, so the fault forms are complicated and various due to different service conditions and geometric structures. Therefore, the fault detection and damaged status assessment of recycled gearboxes become important research topics.The vibration information can effectively reflect the running status of gearbox. Therefore, vibration-based fault diagnosis is the most effective technique in practice. The majority of research is focused on vibration signals analysis in recent years. The advanced signal processing methods [1, 2] have been used to extract useful fault features [3] from vibration signals and assess the damage status. So far, the analysis of modulation sidebands has been an important and successful application in fault feature extraction. Feng et al. [4, 5] proposed the explicit equations to describe the signal models considering both the amplitude and frequency modulation due to gear damage and time-varying transfer paths and calculated the characteristic frequency of planetary gear set with local damage and distributed damage. Inalpolat et al. [6, 7] studied the amplitude modulation and frequency modulation effect due to manufacturing errors of planetary gear set through a mathematical model and dynamic model, and the amplitude modulation effect by carrier rotation was also taken account of.Based on these researches, the vibration signals of planetary gear set have more complex spectral characteristics than fixed shaft gear trains. The damage on gear tooth surface such as wear, pitting, chipping, and cracking will induce the multiplicative amplitude modulation and frequency modulation to the meshing frequency and its harmonics. The time-varying transfer path will also result in a multiplicative amplitude modulation effect on the vibration signals. Furthermore, the measured spectra will contain additional sidebands activity due to certain manufacturing errors, bearings damage, shafts deformations, and strong background noise. The varying operational conditions often lead to nonstationary and nonlinear signals, which make it more difficult to diagnose faults via vibration signals analysis.Model-based dynamic analysis and simulation have been an effective method to analyze the vibration signals. More accurate dynamic models of planetary gear set have been developed in recent years considering more degrees of freedom, different nonlinear factors, damage, and so on. Lin and Parker [8] developed an analytical model of planetary gear set to investigate the natural frequency and vibration modes, which considers two translations and one rotation freedom of each compound. In another paper, Lin and Parker [9] studied the parametric instability caused by the mesh stiffness, which was modeled as rectangular waveform. The gear mesh stiffness is one of the main nonlinear parameters of the dynamic model. The sudden change of gear mesh stiffness due to gear damage will cause periodic impulses in time domain and induce the modulation phenomenon in frequency domain, which is important for fault diagnosis. Yang and Lin [10] proposed the potential energy method to calculate the time-varying mesh stiffness of a one-stage fixed shaft spur gear set, which was assumed to include three components: Hertzian contact energy, bending energy, and axial compressive energy. Later, Tian [11] redefined the potential energy method by taking the shear energy into consideration and investigated the influence of different types of gear damage on the mesh stiffness of a fixed shaft gear pair. Zhou et al. [12] and Chen and Shao [13] made the model more realistic by taking into account the deformation of gear body. Zhou et al. [12], Wu et al. [14, 15], and Tian et al. [16] investigated the effect of crack levels on the mesh stiffness based on the works of Tian [11] and used several statistical indicators to evaluate the changes induced by crack propagation in the simulated vibration signals. Liang et al. [17] extended the model to evaluate the mesh stiffness of an external-internal gear pair of the planetary gear set. The mesh stiffness model mentioned above was simplified as a cantilever beam on the base circle. Actually, the gear tooth starts from the root circle, and there is no uniform equation to describe the tooth curve between the basic circle and the root circle. Later, Liang et al. [18, 19] used straight lines to simplify this part and derived the new equations to calculate the mesh stiffness components. Wan et al. [20] also considered the potential energy stored in this part and analyzed the dynamic response of a single-stage fixed shaft spur gear system considering the influence of gear crack. More works about dynamic analysis of gear set have been reported considering more factors [20, 21], such as different gear damage, gear transmission error, bearing stiffness, and friction.These results can provide useful information for dynamic analysis of damaged gear set. However, most of the works were limited to the simple and one-stage gear set. The gearbox in practice always has compound and multistage gear sets. In the vibration signals collected from a multistage gearbox, there are many close frequency components including different meshing frequencies and their harmonics. So the modulation sidebands will be more complicated and the fault features will be difficult to identify.This paper is organized as follows: In Section1, a review for model-based vibration signal analysis is introduced. In Section 2, a nonlinear dynamic model of a two-stage compound planetary gear set is established. Then, the mesh stiffness of each mesh pair is analytically evaluated using improved potential energy method. Further, the influence of crack propagation on the mesh stiffness is investigated when a crack appears in the sun gear of the second stage. In Section 3, the Runge-Kutta numerical analysis method is used to simulate the vibration response of the gear set in time domain. And also, the corresponding spectral analysis is applied on the simulated signals to obtain the sidebands caused by sun gear damage. The side frequencies are extracted as the useful fault features for fault diagnosis. Further, the influence of crack propagation on the vibration signal is analyzed. In Section 4, the experimental signals of an industrial SD16 planetary gearbox are analyzed to validate the dynamic model and the simulation results. Finally, conclusion and discussion are provided in Section 5. Further improvement will be provided in the future research.


## 2. Nonlinear Dynamic Model of a Two-Stage Planetary Gear Set

Figure1 shows the structure of a two-stage planetary gear set. All the gears are standard involute spur tooth. The first stage (n=1) is a simple planetary gear set with N1 planets p˙i  (N1=3). The second stage (n=2) is a meshed planet planetary gear set with N2 planet-planet mesh pairs p-i-p^i  (N2=3). The planets of each stage distribute equidistantly and have the same parameters. The ring gear r and the carrier c are shared by both stages. In this structure, the sun gear s1 is fixed. The sun gear s2 and carrier c are identified as input member and output member. Tin and Tout are the input motor torque and output load torque.Figure 1
Transmission diagram of the two-stage planetary gear set.The dynamic model of this two-stage planetary gear set is shown in Figure2. A lumped-parameter model is considered in this study. The gears and carrier are treated as perfectly mounted rigid bodies with ideal geometries. The intertooth friction is ignored for simplicity. The gear mesh interactions are modeled by spring-damping structures acting along the line of action, and the periodically time-varying mesh stiffness kjpi and viscous damping coefficient cjpi are considered. The bearings are represented by spring-damping structures, and the bearing stiffness and damping coefficient are set to be constant. In this model, each component has three degrees of freedom: translations in x and y direction, denoted by xj and yj, and rotation around the axis θj (j represents the center components c, r, s1, and s2 and the planets pi(p˙i,p-i,p^i)).Figure 2
Dynamic models for (a) the first stage and (b) the second stage.
(a)
(b)


### 2.1. Equations of Motion

Defineuj=rjθj as the rotational coordinates, where rj (j represents the center components r, s1, and s2 and the planets pi(p˙i,p-i,p^i)) are the base radii of the gears. According to the meshing relation, the relative gear mesh displacements δjpi along the line of action of stage n are expressed as follows:(1)δs1p˙i=-xs1sin⁡φs1p˙i+xp˙isin⁡φs1p˙i+ys1cos⁡φs1p˙i-yp˙icos⁡φs1p˙i+us1+up˙i,δrp˙i=-xrsin⁡φrp˙i+xp˙isin⁡φrp˙i+yrcos⁡φrp˙i-yp˙icos⁡φrp˙i+ur-up˙i,δs2p-i=-xs2sin⁡φs2p-i+xp-isin⁡φs2p-i+ys2cos⁡φs2p-i-yp-icos⁡φs2p-i+us2+up-i,δp-ip^i=xp-isin⁡φp-ip^i-xp^isin⁡φp-ip^i-yp-icos⁡φp-ip^i+yp^icos⁡φp-ip^i-up-i-up^i,δrp^i=xrsin⁡φrp^i-xp^isin⁡φrp^i-yrcos⁡φrp^i+yp^icos⁡φrp^i-ur+up^i,where φjpi=φpi-δjα,δj=-1 for internal r-p˙i mesh pair and external p-i-p^i mesh pair, and δj=1 for internal r-p^i mesh pair and external s1-p˙i and s2-p-i mesh pairs. φpi is the spacing angle of the planets. α is the pressure angle.The dynamic gear mesh force is defined as follows:(2)Fjpi=cjpiδ˙jpi+kjpiδjpi;the bearing forces are defined as follows:(3)Fbjx=cbjx˙j+kbjxj,Fbjy=cbjy˙j+kbjyj,Fuj=cuju˙j+kujuj,where cjpi is the mesh damping coefficient, kjpi is the periodically time-varying mesh stiffness, cbj is the damping coefficient in radial direction, kbj is the radial stiffness of the bearing, cuj is the damping coefficient in torsional direction, and kuj is the torsional stiffness of the bearing.The equations of motion for sun gear can be derived as follows:(4)ms1x¨s1-∑i=1Nsin⁡φs1p˙iFs1p˙i+Fbs1x=0,ms1y¨s1+∑i=1Ncos⁡φs1p˙iFs1p˙i+Fbs1y=0,Is1rs12u¨s1+∑i=1NFs1p˙i+Fus1=0,ms2x¨s2-∑i=1Nsin⁡φs2p-iFs2p-i+Fbs2x=0,ms2y¨s2+∑i=1Ncos⁡φs2p-iFs2p-i+Fbs2y=0,Is2rs22u¨s2+∑i=1NFs2p-i+Fus2=Tinrs2.The equations of motion for ring gear are written as follows:(5)mrx¨r+∑i=1Nsin⁡φrp^iFrp^i-∑i=1Nsin⁡φrp˙iFrp˙i+Fbrx=0,mry¨r-∑i=1Ncos⁡φrp^iFrp^i+∑i=1Ncos⁡φrp˙iFrp˙i+Fbry=0,Irrr2u¨r-∑i=1NFrp^i+∑i=1NFrp˙i+Fur=0.The equations of motion for planets are given as follows:(6)mp˙ix¨p˙i+sin⁡φs1p˙iFs1p˙i+sin⁡φrp˙iFrp˙i=0,mp˙iy¨p˙i-cos⁡φs1p˙iFs1p˙i-cos⁡φrp˙iFrp˙i=0,Ip˙irp˙i2u¨p˙i+Fs1p˙i-Frp˙i=0,mp-ix¨p-i+sin⁡φp-ip^iFp-ip^i+sin⁡φs2p-iFs2p-i=0,mp-iy¨p-i-cos⁡φp-ip^iFp-ip^i-cos⁡φs2p-iFs2p-i=0,Ip-irp-i2u¨p-i+Fs2p-i-Fp-ip^i=0,mp^ix¨p^i-sin⁡φp-ip^iFp-ip^i-sin⁡φrp^iFrp^i=0,mp^iy¨p^i+cos⁡φp-ip^iFp-ip^i+cos⁡φrp^iFrp^i=0,Ip^irp^i2u¨p^i-Fp-ip^i+Frp^i=0,where mj and Ij represent the mass and mass moment of inertia of the gears.For the carrier-planet subsystem,uc=rcθc is defined as the coordinate in place of θc, where rc is the radius from the carrier center to the planet center. The relative displacements between carrier and planet in xc and yc directions are given as follows:(7)δcpix=xc-xpi-sin⁡φpiuc,δcpiy=yc-ypi+cos⁡φpiuc.The differential equations of motion for carrier and planets are expressed as follows:(8)mcx¨c+∑i=1Ncbpiδ˙cpix+kbpiδcpix+cbcx˙c+kbcxc=0,mcy¨c+∑i=1Ncbpiδ˙cpiy+kbpiδcpiy+cbcx˙c+kbcxc=0,Icrc2-∑i=1Ncbpiδ˙cpix+kbpiδcpixsin⁡φpi+∑i=1Ncbpiδ˙cpiy+kbpiδcpiycos⁡φpi+cucu˙c+kucuc=-Toutrc,mpix¨pi-cbpiδ˙cpix+kbpiδcpix=0,mpiy¨pi-cbpiδ˙cpiy+kbpiδbpiy=0,where mc and Ic represent the mass and mass moment of inertia of the carrier.Thus, the overall equations of motion for the two-stage planetary gear set can be constructed systematically in matrix form as follows:(9)MQ¨t+Cm+CbQ˙t+Km+KbQt=T,where (10)M=diagmjmjIjrj2,Q=xjyjujT,Km=Kcp000Kcp˙1Kcp˙2Kcp˙3Kcp-1Kcp-2Kcp-3Kcp^1Kcp^2Kcp^3Krp00Krp˙1Krp˙2Krp˙3000Krp^1Krp^2Krp^3Ks1p0Ks1p˙1Ks1p˙2Ks1p˙3000000Ks2p000Ks2p-1Ks2p-2Ks2p-3000Kp˙100000000Kp˙20000000Kp˙3000000Kp-100Kp-1p^100Kp-200Kp-2p^20Kp-300Kp-3p^3Kp^100Kp^20symmetricKp^3,Kb=diag⁡KbcKbrKbs1Kbs20⋯0,Cb=diag⁡CbcCbrCbs1Cbs20⋯0.In these equations, M is the inertia matrix, Q is the displacement matrix, Cm is the mesh damping matrix, Cb is the bearing damping matrix, Km is the periodically time-varying mesh stiffness matrix, the submatrices of Km are given in Appendix, Kb is the bearing stiffness matrix, and T is the applied external torque matrix.All the parameters of this gear set are provided by the cooperative enterprise according to an industrial SD16 planetary gearbox and are listed in Table1.Table 1
Main parameters of the gear set.
Parameters
Sun gear
Planet
Ring gear
Carrier
s
1
s
2
p
i
p
˙
i
,
p
-
i
,
p
^
i
r
c
Number of teeth
30
30
21
72
78
—
Module (mm)
3.5
3.5
3.5
3.5
—
Pressure angle (deg)
20
20
20
20
—
Pitch diameter (mm)
105
105
73.5
252
273
—
Base diameter (mm)
98.668
98.668
69.067
236.803
256.536
—
Addendum circle diameter (mm)
112
112
80.5
245.734
266.677
—
Root diameter (mm)
96.25
96.25
64.75
260.75
281.75
Width of teeth (mm)
40.5
40.5
38.5
92
—
Material
20 NiCrMoH
20 NiCrMoH
20 NiCrMoH
42 CrMo
ZG40Mn2
Young’s modulus (Pa)
2.06 × 1011
2.06 × 1011
2.06 × 1011
2.12 × 1011
2.02 × 1011
Poisson’s ratio
0.30
0.30
0.30
0.28
0.30
Mass (kg)
2.11
2.11
0.92
11.89
15.54
Mass moment of inertia (kg⋅mm2)
0.360 × 104
0.360 × 104
0.081 × 104
24.095 × 104
12.729 × 104
Input rotation speed (rpm)
700
Input motor torque (N⋅m)
96
Output load torque (N⋅m)
200
Transmission ratio
2.08
Radial stiffness of the bearing (N/m)
108
Damping coefficient of the bearing (N⋅s/m)
105


### 2.2. Mesh Stiffness Evaluation

In this section, the mesh stiffness of each mesh pair for this two-stage gear set will be calculated, respectively, based on the improved model proposed by Liang et al. [18]. The mesh stiffness model is simplified as a cantilever beam on the root circle. The potential energy method is used to analyze the mesh stiffness, which is assumed to include four components: Hertzian contact energy, bending energy, shear energy, and axial compressive energy. According to the design parameters of this gear set, the root circles of gears are smaller than the base circles for the external-external mesh pairs. The tooth curve between the base circle and the root circle is not an involute curve and there is no uniform equation to describe this part. So this part is simplified as straight lines [18], as shown in Figure 3.Figure 3
The mesh stiffness beam model for healthy gear.Thus, for the single-tooth pair meshing duration, the total mesh stiffness can be calculated as follows:(11)kjpi=11/kh+1/kb1+1/ks1+1/ka1+1/kb2+1/ks2+1/ka2,where subscripts 1 and 2 represent the driving gear and driven gear. kh, kb, ks, and ka represent the Hertzian contact stiffness, bending stiffness, shear stiffness, and axial compressive stiffness, respectively.For the double-tooth pair mesh duration, the total mesh stiffness can be expressed as follows:(12)kjpi=∑i=1211/kh,i+1/kb1,i+1/ks1,i+1/ka1,i+1/kb2,i+1/ks2,i+1/ka2,i,where i=1 represents the first pair of meshing teeth and i=2 represents the second. The detail of derivations of each stiffness component was given in Liang et al. [18].It is essential to correctly define the phasing relationships between each mesh pair. References [17–19, 22] provided an analytical calculation of the relative phase for a one-stage and simple planetary gear set. Extending the methods proposed in [17–19, 22] and combining the parameters listed in Table 1, the relative phases of this two-stage planetary gear set can be calculated. Figure 4 shows a meshing sketch for this gear set and only one of the planet branches is shown for simplicity purpose. The spacing angles of p˙1 and p-1 are set to be 0. NiMi¯ and BiEi¯(i=1,2,3,…) denote the theoretical meshing line and practical action line. Pi(i=1,2,3,…) is the pitch point of each mesh pair. As shown in Figure 1, the sun gear s1 is fixed and s2 is defined as the input member. Assuming the sun gear s2 rotates in the counter-clockwise direction, the planets p-i, p^i, and p˙i are in clockwise, counter-clockwise, and counter-clockwise rotation directions, respectively.Figure 4
A meshing sketch for the two-stage planetary gear set: (a) the first stage and (b) the second stage.
(a)
(b)Table2 shows the relative phases of this gear set. The superscripts represent the first stage and the second stage of the gear set, respectively. γsn′ and γsn′′ denote the relative phase between the nth sun-planet mesh pair and the first sun-planet mesh pair. γrn′ and γrn′′ are the relative phase between the nth ring-planet mesh pair and the first ring-planet mesh pair. γpn′′ is the relative phase between the nth planet-planet mesh pair and the first planet-planet mesh pair. γrs′ and γrs′′ are the relative phase between the nth ring-planet mesh pair and sun-planet mesh pair. γps′′ is the relative phase between the nth planet-planet mesh pair and sun-planet mesh pair. “+” dictate the phase lag, and “−” dictate the phase lead.Table 2
Relative phases of the gear set.
Stage 1
γ
s
1
′
γ
s
2
′
γ
s
3
′
γ
r
1
′
γ
r
2
′
γ
r
3
′
γ
r
s
′
0
0
0
0
0
0
0
Stage 2
γ
s
1
′
′
γ
s
2
′
′
γ
s
3
′
′
γ
p
1
′
′
γ
p
2
′
′
γ
p
3
′
′
γ
r
1
′
′
γ
r
2
′
′
γ
r
3
′
′
γ
p
s
′
′
γ
r
s
′
′
0
0
0
0
0
0
0
0
0
-
0.13
+0.38We use MATLAB programs to obtain the numerical values of total mesh stiffness as a function of timet, which are plotted in Figures 5 and 6. It is observed that the mesh stiffness varies periodically. There exists a discontinuous point in a mesh period, which is the transition point from double-tooth pair to single-tooth pair. For the first stage, the mean stiffness values are 1.676 × 109 N/m for sun-planet pair and 1.677 × 109 N/m for ring-planet pair. γrs′=0 indicates the nth sun-planet mesh pair and the ring-planet mesh pair mesh at the pitch point simultaneously. For the second stage, the mean stiffness values are 1.676 × 109 N/m for sun-planet pair, 1.567 × 109 N/m for planet-planet pair, and 1.573 × 109 N/m for ring-planet pair. The phase difference between the nth planet-planet mesh pair and sun-planet mesh pair is γps′′Tm in time, where Tm is the mesh period of the second stage. Similarly, the phase difference between the nth ring-planet mesh pair and sun-planet mesh pair can be expressed as γrs′′Tm in time.Figure 5
Total mesh stiffness of the perfect first stage versus timet within 0.1 s: (a) sun-planet pair and (b) ring-planet pair.
(a)
(b)Figure 6
Total mesh stiffness of the perfect second stage versus timet within 0.1 s: (a) sun-planet pair, (b) planet-planet pair, and (c) ring-planet pair.Liang et al. [18] assumed the crack propagation path to be a straight line starting from the critical area of the tooth. As shown in Figure 7, the crack propagates along the straight line until it reaches the central line. The intersection angle v between the crack path and the central line is set to 45°. Then, the crack changes its direction until a sudden breakage appears. The Hertzian and axial compressive stiffness remain the same; however, the bending and shear stiffness will change due to the crack. In this section, the influence of crack on the mesh stiffness of a sun/planet pair in the second stage is investigated when the crack appears in the sun gear. Based on the crack model provided above, five crack levels are evaluated: 10%, 30%, 50%, 70%, and 80% cracks. When the crack line q1 reaches the central line, the crack level is set to 50%, and when a sudden break occurs, the corresponding crack level is 80%. Table 3 shows the five crack levels and corresponding crack length in the sun gear s2.Table 3
Crack levels and corresponding crack length in the sun gear of second stage.
Crack levels
Crack length (mm)
1
10%
q
1
=
0.942
2
30%
q
1
=
2.826
3
50%
q
1
=
4.71
4
70%
q
1
=
4.71, q2=1.884
5
80%
Tooth missingFigure 7
The cantilever beam model for local damaged gear.The total mesh stiffness curves of the sun-planet pair in second stage with different crack lengths are shown in Figure8. It can be observed that the mesh stiffness greatly decreases when the cracked tooth is in meshing. The influences spread over the mating duration of the cracked gear, which corresponds to 0.0089 s. As the size of crack grows, the mesh stiffness reduces correspondingly. When the sudden break appears, the mean mesh stiffness reduces to 42% of the original. This is important for dynamic analysis of the gear set.Figure 8
The mesh stiffness of sun-planet pair with different crack levels.The mesh damping coefficient is set to be proportional to the total mesh stiffness; that is,(13)cjpi=μjpikjpi,where μjpi is the scale constant measured in seconds. For this gear set, μs1p˙i=μs2p-i=2.70×10-6 s,  μp-ip^i=2.40×10-6 s,  μrp˙i=3.16×10-6 s, and  μrp^i=3.26×10-6 s.


## 2.1. Equations of Motion

Defineuj=rjθj as the rotational coordinates, where rj (j represents the center components r, s1, and s2 and the planets pi(p˙i,p-i,p^i)) are the base radii of the gears. According to the meshing relation, the relative gear mesh displacements δjpi along the line of action of stage n are expressed as follows:(1)δs1p˙i=-xs1sin⁡φs1p˙i+xp˙isin⁡φs1p˙i+ys1cos⁡φs1p˙i-yp˙icos⁡φs1p˙i+us1+up˙i,δrp˙i=-xrsin⁡φrp˙i+xp˙isin⁡φrp˙i+yrcos⁡φrp˙i-yp˙icos⁡φrp˙i+ur-up˙i,δs2p-i=-xs2sin⁡φs2p-i+xp-isin⁡φs2p-i+ys2cos⁡φs2p-i-yp-icos⁡φs2p-i+us2+up-i,δp-ip^i=xp-isin⁡φp-ip^i-xp^isin⁡φp-ip^i-yp-icos⁡φp-ip^i+yp^icos⁡φp-ip^i-up-i-up^i,δrp^i=xrsin⁡φrp^i-xp^isin⁡φrp^i-yrcos⁡φrp^i+yp^icos⁡φrp^i-ur+up^i,where φjpi=φpi-δjα,δj=-1 for internal r-p˙i mesh pair and external p-i-p^i mesh pair, and δj=1 for internal r-p^i mesh pair and external s1-p˙i and s2-p-i mesh pairs. φpi is the spacing angle of the planets. α is the pressure angle.The dynamic gear mesh force is defined as follows:(2)Fjpi=cjpiδ˙jpi+kjpiδjpi;the bearing forces are defined as follows:(3)Fbjx=cbjx˙j+kbjxj,Fbjy=cbjy˙j+kbjyj,Fuj=cuju˙j+kujuj,where cjpi is the mesh damping coefficient, kjpi is the periodically time-varying mesh stiffness, cbj is the damping coefficient in radial direction, kbj is the radial stiffness of the bearing, cuj is the damping coefficient in torsional direction, and kuj is the torsional stiffness of the bearing.The equations of motion for sun gear can be derived as follows:(4)ms1x¨s1-∑i=1Nsin⁡φs1p˙iFs1p˙i+Fbs1x=0,ms1y¨s1+∑i=1Ncos⁡φs1p˙iFs1p˙i+Fbs1y=0,Is1rs12u¨s1+∑i=1NFs1p˙i+Fus1=0,ms2x¨s2-∑i=1Nsin⁡φs2p-iFs2p-i+Fbs2x=0,ms2y¨s2+∑i=1Ncos⁡φs2p-iFs2p-i+Fbs2y=0,Is2rs22u¨s2+∑i=1NFs2p-i+Fus2=Tinrs2.The equations of motion for ring gear are written as follows:(5)mrx¨r+∑i=1Nsin⁡φrp^iFrp^i-∑i=1Nsin⁡φrp˙iFrp˙i+Fbrx=0,mry¨r-∑i=1Ncos⁡φrp^iFrp^i+∑i=1Ncos⁡φrp˙iFrp˙i+Fbry=0,Irrr2u¨r-∑i=1NFrp^i+∑i=1NFrp˙i+Fur=0.The equations of motion for planets are given as follows:(6)mp˙ix¨p˙i+sin⁡φs1p˙iFs1p˙i+sin⁡φrp˙iFrp˙i=0,mp˙iy¨p˙i-cos⁡φs1p˙iFs1p˙i-cos⁡φrp˙iFrp˙i=0,Ip˙irp˙i2u¨p˙i+Fs1p˙i-Frp˙i=0,mp-ix¨p-i+sin⁡φp-ip^iFp-ip^i+sin⁡φs2p-iFs2p-i=0,mp-iy¨p-i-cos⁡φp-ip^iFp-ip^i-cos⁡φs2p-iFs2p-i=0,Ip-irp-i2u¨p-i+Fs2p-i-Fp-ip^i=0,mp^ix¨p^i-sin⁡φp-ip^iFp-ip^i-sin⁡φrp^iFrp^i=0,mp^iy¨p^i+cos⁡φp-ip^iFp-ip^i+cos⁡φrp^iFrp^i=0,Ip^irp^i2u¨p^i-Fp-ip^i+Frp^i=0,where mj and Ij represent the mass and mass moment of inertia of the gears.For the carrier-planet subsystem,uc=rcθc is defined as the coordinate in place of θc, where rc is the radius from the carrier center to the planet center. The relative displacements between carrier and planet in xc and yc directions are given as follows:(7)δcpix=xc-xpi-sin⁡φpiuc,δcpiy=yc-ypi+cos⁡φpiuc.The differential equations of motion for carrier and planets are expressed as follows:(8)mcx¨c+∑i=1Ncbpiδ˙cpix+kbpiδcpix+cbcx˙c+kbcxc=0,mcy¨c+∑i=1Ncbpiδ˙cpiy+kbpiδcpiy+cbcx˙c+kbcxc=0,Icrc2-∑i=1Ncbpiδ˙cpix+kbpiδcpixsin⁡φpi+∑i=1Ncbpiδ˙cpiy+kbpiδcpiycos⁡φpi+cucu˙c+kucuc=-Toutrc,mpix¨pi-cbpiδ˙cpix+kbpiδcpix=0,mpiy¨pi-cbpiδ˙cpiy+kbpiδbpiy=0,where mc and Ic represent the mass and mass moment of inertia of the carrier.Thus, the overall equations of motion for the two-stage planetary gear set can be constructed systematically in matrix form as follows:(9)MQ¨t+Cm+CbQ˙t+Km+KbQt=T,where (10)M=diagmjmjIjrj2,Q=xjyjujT,Km=Kcp000Kcp˙1Kcp˙2Kcp˙3Kcp-1Kcp-2Kcp-3Kcp^1Kcp^2Kcp^3Krp00Krp˙1Krp˙2Krp˙3000Krp^1Krp^2Krp^3Ks1p0Ks1p˙1Ks1p˙2Ks1p˙3000000Ks2p000Ks2p-1Ks2p-2Ks2p-3000Kp˙100000000Kp˙20000000Kp˙3000000Kp-100Kp-1p^100Kp-200Kp-2p^20Kp-300Kp-3p^3Kp^100Kp^20symmetricKp^3,Kb=diag⁡KbcKbrKbs1Kbs20⋯0,Cb=diag⁡CbcCbrCbs1Cbs20⋯0.In these equations, M is the inertia matrix, Q is the displacement matrix, Cm is the mesh damping matrix, Cb is the bearing damping matrix, Km is the periodically time-varying mesh stiffness matrix, the submatrices of Km are given in Appendix, Kb is the bearing stiffness matrix, and T is the applied external torque matrix.All the parameters of this gear set are provided by the cooperative enterprise according to an industrial SD16 planetary gearbox and are listed in Table1.Table 1
Main parameters of the gear set.
Parameters
Sun gear
Planet
Ring gear
Carrier
s
1
s
2
p
i
p
˙
i
,
p
-
i
,
p
^
i
r
c
Number of teeth
30
30
21
72
78
—
Module (mm)
3.5
3.5
3.5
3.5
—
Pressure angle (deg)
20
20
20
20
—
Pitch diameter (mm)
105
105
73.5
252
273
—
Base diameter (mm)
98.668
98.668
69.067
236.803
256.536
—
Addendum circle diameter (mm)
112
112
80.5
245.734
266.677
—
Root diameter (mm)
96.25
96.25
64.75
260.75
281.75
Width of teeth (mm)
40.5
40.5
38.5
92
—
Material
20 NiCrMoH
20 NiCrMoH
20 NiCrMoH
42 CrMo
ZG40Mn2
Young’s modulus (Pa)
2.06 × 1011
2.06 × 1011
2.06 × 1011
2.12 × 1011
2.02 × 1011
Poisson’s ratio
0.30
0.30
0.30
0.28
0.30
Mass (kg)
2.11
2.11
0.92
11.89
15.54
Mass moment of inertia (kg⋅mm2)
0.360 × 104
0.360 × 104
0.081 × 104
24.095 × 104
12.729 × 104
Input rotation speed (rpm)
700
Input motor torque (N⋅m)
96
Output load torque (N⋅m)
200
Transmission ratio
2.08
Radial stiffness of the bearing (N/m)
108
Damping coefficient of the bearing (N⋅s/m)
105


## 2.2. Mesh Stiffness Evaluation

In this section, the mesh stiffness of each mesh pair for this two-stage gear set will be calculated, respectively, based on the improved model proposed by Liang et al. [18]. The mesh stiffness model is simplified as a cantilever beam on the root circle. The potential energy method is used to analyze the mesh stiffness, which is assumed to include four components: Hertzian contact energy, bending energy, shear energy, and axial compressive energy. According to the design parameters of this gear set, the root circles of gears are smaller than the base circles for the external-external mesh pairs. The tooth curve between the base circle and the root circle is not an involute curve and there is no uniform equation to describe this part. So this part is simplified as straight lines [18], as shown in Figure 3.Figure 3
The mesh stiffness beam model for healthy gear.Thus, for the single-tooth pair meshing duration, the total mesh stiffness can be calculated as follows:(11)kjpi=11/kh+1/kb1+1/ks1+1/ka1+1/kb2+1/ks2+1/ka2,where subscripts 1 and 2 represent the driving gear and driven gear. kh, kb, ks, and ka represent the Hertzian contact stiffness, bending stiffness, shear stiffness, and axial compressive stiffness, respectively.For the double-tooth pair mesh duration, the total mesh stiffness can be expressed as follows:(12)kjpi=∑i=1211/kh,i+1/kb1,i+1/ks1,i+1/ka1,i+1/kb2,i+1/ks2,i+1/ka2,i,where i=1 represents the first pair of meshing teeth and i=2 represents the second. The detail of derivations of each stiffness component was given in Liang et al. [18].It is essential to correctly define the phasing relationships between each mesh pair. References [17–19, 22] provided an analytical calculation of the relative phase for a one-stage and simple planetary gear set. Extending the methods proposed in [17–19, 22] and combining the parameters listed in Table 1, the relative phases of this two-stage planetary gear set can be calculated. Figure 4 shows a meshing sketch for this gear set and only one of the planet branches is shown for simplicity purpose. The spacing angles of p˙1 and p-1 are set to be 0. NiMi¯ and BiEi¯(i=1,2,3,…) denote the theoretical meshing line and practical action line. Pi(i=1,2,3,…) is the pitch point of each mesh pair. As shown in Figure 1, the sun gear s1 is fixed and s2 is defined as the input member. Assuming the sun gear s2 rotates in the counter-clockwise direction, the planets p-i, p^i, and p˙i are in clockwise, counter-clockwise, and counter-clockwise rotation directions, respectively.Figure 4
A meshing sketch for the two-stage planetary gear set: (a) the first stage and (b) the second stage.
(a)
(b)Table2 shows the relative phases of this gear set. The superscripts represent the first stage and the second stage of the gear set, respectively. γsn′ and γsn′′ denote the relative phase between the nth sun-planet mesh pair and the first sun-planet mesh pair. γrn′ and γrn′′ are the relative phase between the nth ring-planet mesh pair and the first ring-planet mesh pair. γpn′′ is the relative phase between the nth planet-planet mesh pair and the first planet-planet mesh pair. γrs′ and γrs′′ are the relative phase between the nth ring-planet mesh pair and sun-planet mesh pair. γps′′ is the relative phase between the nth planet-planet mesh pair and sun-planet mesh pair. “+” dictate the phase lag, and “−” dictate the phase lead.Table 2
Relative phases of the gear set.
Stage 1
γ
s
1
′
γ
s
2
′
γ
s
3
′
γ
r
1
′
γ
r
2
′
γ
r
3
′
γ
r
s
′
0
0
0
0
0
0
0
Stage 2
γ
s
1
′
′
γ
s
2
′
′
γ
s
3
′
′
γ
p
1
′
′
γ
p
2
′
′
γ
p
3
′
′
γ
r
1
′
′
γ
r
2
′
′
γ
r
3
′
′
γ
p
s
′
′
γ
r
s
′
′
0
0
0
0
0
0
0
0
0
-
0.13
+0.38We use MATLAB programs to obtain the numerical values of total mesh stiffness as a function of timet, which are plotted in Figures 5 and 6. It is observed that the mesh stiffness varies periodically. There exists a discontinuous point in a mesh period, which is the transition point from double-tooth pair to single-tooth pair. For the first stage, the mean stiffness values are 1.676 × 109 N/m for sun-planet pair and 1.677 × 109 N/m for ring-planet pair. γrs′=0 indicates the nth sun-planet mesh pair and the ring-planet mesh pair mesh at the pitch point simultaneously. For the second stage, the mean stiffness values are 1.676 × 109 N/m for sun-planet pair, 1.567 × 109 N/m for planet-planet pair, and 1.573 × 109 N/m for ring-planet pair. The phase difference between the nth planet-planet mesh pair and sun-planet mesh pair is γps′′Tm in time, where Tm is the mesh period of the second stage. Similarly, the phase difference between the nth ring-planet mesh pair and sun-planet mesh pair can be expressed as γrs′′Tm in time.Figure 5
Total mesh stiffness of the perfect first stage versus timet within 0.1 s: (a) sun-planet pair and (b) ring-planet pair.
(a)
(b)Figure 6
Total mesh stiffness of the perfect second stage versus timet within 0.1 s: (a) sun-planet pair, (b) planet-planet pair, and (c) ring-planet pair.Liang et al. [18] assumed the crack propagation path to be a straight line starting from the critical area of the tooth. As shown in Figure 7, the crack propagates along the straight line until it reaches the central line. The intersection angle v between the crack path and the central line is set to 45°. Then, the crack changes its direction until a sudden breakage appears. The Hertzian and axial compressive stiffness remain the same; however, the bending and shear stiffness will change due to the crack. In this section, the influence of crack on the mesh stiffness of a sun/planet pair in the second stage is investigated when the crack appears in the sun gear. Based on the crack model provided above, five crack levels are evaluated: 10%, 30%, 50%, 70%, and 80% cracks. When the crack line q1 reaches the central line, the crack level is set to 50%, and when a sudden break occurs, the corresponding crack level is 80%. Table 3 shows the five crack levels and corresponding crack length in the sun gear s2.Table 3
Crack levels and corresponding crack length in the sun gear of second stage.
Crack levels
Crack length (mm)
1
10%
q
1
=
0.942
2
30%
q
1
=
2.826
3
50%
q
1
=
4.71
4
70%
q
1
=
4.71, q2=1.884
5
80%
Tooth missingFigure 7
The cantilever beam model for local damaged gear.The total mesh stiffness curves of the sun-planet pair in second stage with different crack lengths are shown in Figure8. It can be observed that the mesh stiffness greatly decreases when the cracked tooth is in meshing. The influences spread over the mating duration of the cracked gear, which corresponds to 0.0089 s. As the size of crack grows, the mesh stiffness reduces correspondingly. When the sudden break appears, the mean mesh stiffness reduces to 42% of the original. This is important for dynamic analysis of the gear set.Figure 8
The mesh stiffness of sun-planet pair with different crack levels.The mesh damping coefficient is set to be proportional to the total mesh stiffness; that is,(13)cjpi=μjpikjpi,where μjpi is the scale constant measured in seconds. For this gear set, μs1p˙i=μs2p-i=2.70×10-6 s,  μp-ip^i=2.40×10-6 s,  μrp˙i=3.16×10-6 s, and  μrp^i=3.26×10-6 s.


## 3. Analysis of the Simulated Signals

To illustrate the dynamic response of the system, the MATLAB’s ode15s function is used to solve the equations of motion derived above. Then, the vibration responses of each component of this gear set as a function of time are numerically simulated. Applying the fft function based on the Fast Fourier Transform algorithm provided by MATLAB toolbox, the corresponding Fourier spectrum can be obtained. According to the parameters provided in Table1, the dominant frequencies of this gear set are estimated within 2000 Hz. Considering the sample frequency options set in the experimental equipment, the sample frequency used in this simulation is selected as 5120 Hz.The vibration velocity responses of healthy sun gearss1 and s2 are shown, respectively, in Figure 9. The periodical impulses caused by time-varying mesh stiffness are obvious in the vibration signals. And the duration between every two impulses of s1 is equal to 0.00606 s, while the impulse duration of s2 is 0.00547 s. Figure 10 shows the Fourier spectrum of healthy gear set. As expected, the spectrum comprises four frequency components: the fundamental meshing frequency of first stage  f1=168Hz, the fundamental meshing frequency of second stage f2=182Hz, and their harmonics  kf1,kf2k=1,2,…, in which the frequency component 8f2 is dominant. So, for the two-stage planetary gear set, the meshing frequencies of each stage are very close to each other, which can be expressed as follows:(14)fmesh=zsunnsun-carrier60;in this equation, zsun is the number of sun gear teeth and nsun-carrier is the rotation speed of sun gear relative to the carrier.Figure 9
Vibration velocity responses of healthy sun gears: (a) sun gears1 and (b) sun gear s2.
(a)
(b)Figure 10
Fourier spectrum of healthy gear set.Supposing a 70% level (q1=4.71, q2=1.884)  crack appears on the sun gear of the first stage, the dynamic response of sun gear s1 and the corresponding Fourier spectrum of the gear set are shown in Figure 11, respectively. Compared with Figure 9(a), when the cracked sun gear contacts with the mating perfect planet, a sudden change will occur in time domain signal. As the meshing of damaged sun gear, the impulse appears at a time period of 0.060 s. Compared with Figure 10, many spectral peaks appear as sidebands around the meshing frequency components. And the amplitudes of these sidebands are small compared with the meshing frequencies. In order to analyze the details, the spectrum is enlarged in Figure 11(c). It is observed that the sidebands center around the frequency components kf1k=1,2,…, at the frequencies of kf1±mfs1k=1,2,…,m=1,2,…, and fs1 is equal to 16.8 Hz. So the side frequency components mfs1m=1,2,… indicate that the gear set have fault, and the fault exists on the sun gear s1.Figure 11
Vibration response of the gear set with 70% level crack on sun gears1: (a) waveform of cracked sun gear s1, (b) Fourier spectrum, and (c) enlarged Fourier spectrum.
(a)
(b)
(c)Similarly, when the 70% level (q1=4.71, q2=1.884) crack exists on the surface of sun gear in the second stage, the contact of cracked sun gear tooth with mating planets leads to obvious shocks at a repeating period, which is ΔT=0.055 s, as shown in Figure 12(a). The Fourier spectrum of the gear set with cracked sun gear s2, as shown in Figure 12(b), has almost the same frequency components with Figure 11(b): the fundamental meshing frequency of the first stage f1=168Hz, the fundamental meshing frequency of the second stage f2=182Hz, their harmonics kf1,kf2k=1,2,…, and the side frequencies. To illustrate the details, the spectrum is enlarged in Figure 12(c). It is observed that almost all the side peaks appear at the frequencies of kf2±mfs2k=1,2,…,m=1,2,…, and fs2 is 18.2 Hz. As expected, the side frequency components mfs2(m=1,2,…) contain the fault information, and the fault exists on the sun gear s2.Figure 12
Vibration response of the gear set with 70% level crack on sun gears2: (a) waveform of cracked sun gear s2, (b) Fourier spectrum, and (c) enlarged Fourier spectrum.
(a)
(b)
(c)As described in Section2, the crack grows from 10% (in the early stage) to 70% (close to a sudden breakage) with an increment percentage of 20% (1.884 mm). The enlarged Fourier spectrum of the gear set with different crack levels are shown in Figure 13. Clearly, the influence of the damaged sun gear is not obvious, and the amplitudes of sidebands are too small to be identified easily when the crack is in the early stage (the crack level is lower than 10%). As the crack level increases, the amplitude of sidebands becomes stronger, but the sideband structures remain the same.Figure 13
Enlarged Fourier spectrum of the gear set when the crack levels are 10%, 30%, 50%, and 70%.
(a)
(b)
(c)
(d)According to the analysis provided above, the contact of damaged gear tooth with mating gears will lead to periodical impulses in time domain, and the impulse period of each gear component is different and related with the system parameters. Further, the shocks will induce the multiplicative modulation to the meshing frequencies in frequency domain and lead to the sidebands center around the corresponding meshing frequency components. These side frequency components can provide useful information for accurately locating and reliably assessing the damage status of the system. Given the parameters and transmission structure of the planetary gear set, the fault characteristic frequency of sun gear can be derived as follows:(15)fsun=Nplanetnsun-carrier60. According to (14), the equation can also be written as follows:(16)fsun=Nplanetfmeshzsun, where Nplanet is the number of planets mating with corresponding sun gear.


## 4. Experimental Results

In this section, we use three experimental signals to validate the dynamic model provided above and fault characteristics of the multistage planetary gear set obtained from the simulation. One is the vibration signal collected from a new industrial SD16 planetary gearbox, which is considered as the baseline signal. The second one is the signal collected from the same gearbox with naturally broken sun gears1. The third one is the signal collected from the same gearbox with naturally broken sun gear s2.


### 4.1. Experimental Specification

The planetary gearbox test rig is shown in Figure14 and the SD16 planetary gearbox used in the experiments is shown in Figure 15(a). The fault experiments are done on the second row of this gearbox, as shown in Figure 15(b), which has the same structure and parameters as the two-stage model provided in Section 2. A three-axis accelerometer of model KD1002S is mounted on the gearbox casing. The data and vibration signals are measured by CRAS vibration signal analysis system. A 12.8 s span data is recorded with a sampling frequency of 5120 Hz under each experimental condition. During the experiments, the input speed is set to 700 rpm, and a load of 200 N·m is applied to the output shaft.Figure 14
(a) Planetary gearbox test rig and (b) structure diagram of the test rig.
(a)
(b)Figure 15
(a) SD16 planetary gearbox, (b) the second row of this gearbox, and (c) sun gear with broken tooth.
(a)
(b)
(c)For the damaged sun gears1 experiment, the perfect sun gear of the first stage (s1) is replaced by the broken one to mimic the common fault status of the gearbox, as shown in Figure 15(c). Because the fault characteristic frequencies of s1 are mfs1m=1,2,3,…  in simulation, these useful fault features are expected to be detected. For the damaged sun gear s2 experiment, the perfect sun gear of the second stage (s2) is replaced by a broken one, whereas all the other parts are normal. The fault characteristic frequencies of s2, mfs2m=1,2,3,… in simulation, are expected to be extracted.


### 4.2. Normal Vibration Signal Analysis

For a long time span, the input speed and load will fluctuate around the set value during the experiment. This will make the measured Fourier spectral lines deviate the theoretical values dramatically. However, from a short time view of point, the signal can be assumed to be relative stationary [4]. For this reason, we extract 1 s long signal samples from the 12.8 s span data to analyze the Fourier spectral structure, so as to obtain more accurate spectral lines.The waveform and corresponding Fourier spectrum of the baseline signal sample collected from new SD16 planetary gearbox are shown in Figure16. The measured spectrum contains the fundamental meshing frequency of the first stage (f1), about 167.3 Hz, the fundamental meshing frequency of the second stage (f2), about 180.6 Hz, and their harmonics (kf1,kf2(k=1,2,3,…)). Some sidebands appear around the main frequency components, but most of them are irregular and the amplitudes are small. So this does not indicate that the gearbox is damaged.Figure 16
Waveform and corresponding Fourier spectrum of the baseline signal collected from the new SD16 planetary gearbox.
(a)
(b)


### 4.3. Damaged Sun Gear Signal Analysis

Figure17 shows the waveform and Fourier spectrum of the signal sample (the sample length is 1 s) from the same gearbox with naturally broken sun gear s1, respectively. Compared with Figure 16(a), the kurtosis indicator of the signal sample increases from 3.01 to 4.04. Several impulses are visible in the time domain, but the accurate impulse period is difficult to evaluate due to strong back noise. The dominant spectral peaks appear at the frequency components: the fundamental meshing frequency of the first stage (f1), the fundamental meshing frequency of the second stage (f2), their harmonics (kf1,kf2(k=1,2,3,…)), and the sidebands, as the baseline spectrum, but the magnitude of the spectral peaks increases significantly. All the signs indicate the anomalies of the gearbox. As described in the simulation, the sidebands are associated with the transient impulses and contain useful information of damage location. So the zoomed-in Fourier spectrum is shown in Figure 17(c). It can be observed that the sidebands appear around the meshing frequencies of the first stage, at the frequencies of kf1±mfs1 (k=1,2,…,m=1,2,…). To further validate the fault features, the spectral kurtosis signal processing method [23] is employed in this section to extract the transient impulses. Figure 18 presents the spectral kurtosis for the signal sample (the sample length is 5 s). The maximum spectral kurtosis appears at level 5.5 for the central frequency of 800 Hz. The filtered fault signal from the optimal band-pass filter around 800 Hz and the corresponding squared envelope spectrum of the filtered signal are shown in Figure 19. The periodic impulses are obvious in time domain, and a pattern of 0.062 s is confirmed. Frequency of 16.23 Hz and its harmonics exist in the frequency domain, which correspond with the calculation of (16) for the fault characteristic frequency of the sun gear in the first stage mfs1  (m=1,2,3,…). So we can conclude that the gear s1 has fault. The conclusion matches the results obtained from the simulation.Figure 17
Vibration signal from the gearbox with naturally broken sun gears1: (a) waveform, (b) Fourier spectrum, and (c) zoomed-in Fourier spectrum.
(a)
(b)
(c)Figure 18
Spectral kurtosis for the signal sample in Figure17.Figure 19
(a) Filtered fault signal from the optimal band-pass filter in the case of the damaged sun gears1 and (b) squared envelope spectrum of the filtered signal.
(a)
(b)Figure20 shows the waveform and Fourier spectrum of the signal sample (the sample length is 1 s) from the broken sun gear s2 case. The kurtosis value of the signal sample increases from 3.01 to 5.24. More sidebands appear in the frequency band from 500 to 1000 Hz, and most of them have bigger amplitudes than the baseline signal. As shown in Figure 20(c), the sidebands center around the meshing frequencies of the second stage, at the frequencies of kf2±mfs2 (k=1,2,…,m=1,2,…). The corresponding spectral kurtosis is presented in Figure 21. The maximum spectral kurtosis appears at level 6.5 for the central frequency of 80 Hz. The filtered fault signal from the optimal band-pass filter around 80 Hz and the corresponding squared envelope spectrum of the filtered signal are shown in Figure 22. The periodic impulses are evident, and a pattern of 0.056 s is confirmed. Frequency of 17.67 Hz and its harmonics exist in the frequency domain, which correspond with the calculation of (16) for the fault characteristic frequency of the sun gear in the second stage mfs2  (m=1,2,3,…). All the signs indicate that the gearbox is fault, and the damage occurs on the sun gear of the second stage. The results match well with the simulation.Figure 20
Vibration signal from the gearbox with naturally broken sun gears2: (a) waveform, (b) Fourier spectrum, and (c) zoomed-in Fourier spectrum.
(a)
(b)
(c)Figure 21
Spectral kurtosis for the signal sample in Figure20.Figure 22
(a) Filtered fault signal from the optimal band-pass filter in the case of the damaged sun gears2 and (b) squared envelope spectrum of the filtered signal.
(a)
(b)


## 4.1. Experimental Specification

The planetary gearbox test rig is shown in Figure14 and the SD16 planetary gearbox used in the experiments is shown in Figure 15(a). The fault experiments are done on the second row of this gearbox, as shown in Figure 15(b), which has the same structure and parameters as the two-stage model provided in Section 2. A three-axis accelerometer of model KD1002S is mounted on the gearbox casing. The data and vibration signals are measured by CRAS vibration signal analysis system. A 12.8 s span data is recorded with a sampling frequency of 5120 Hz under each experimental condition. During the experiments, the input speed is set to 700 rpm, and a load of 200 N·m is applied to the output shaft.Figure 14
(a) Planetary gearbox test rig and (b) structure diagram of the test rig.
(a)
(b)Figure 15
(a) SD16 planetary gearbox, (b) the second row of this gearbox, and (c) sun gear with broken tooth.
(a)
(b)
(c)For the damaged sun gears1 experiment, the perfect sun gear of the first stage (s1) is replaced by the broken one to mimic the common fault status of the gearbox, as shown in Figure 15(c). Because the fault characteristic frequencies of s1 are mfs1m=1,2,3,…  in simulation, these useful fault features are expected to be detected. For the damaged sun gear s2 experiment, the perfect sun gear of the second stage (s2) is replaced by a broken one, whereas all the other parts are normal. The fault characteristic frequencies of s2, mfs2m=1,2,3,… in simulation, are expected to be extracted.


## 4.2. Normal Vibration Signal Analysis

For a long time span, the input speed and load will fluctuate around the set value during the experiment. This will make the measured Fourier spectral lines deviate the theoretical values dramatically. However, from a short time view of point, the signal can be assumed to be relative stationary [4]. For this reason, we extract 1 s long signal samples from the 12.8 s span data to analyze the Fourier spectral structure, so as to obtain more accurate spectral lines.The waveform and corresponding Fourier spectrum of the baseline signal sample collected from new SD16 planetary gearbox are shown in Figure16. The measured spectrum contains the fundamental meshing frequency of the first stage (f1), about 167.3 Hz, the fundamental meshing frequency of the second stage (f2), about 180.6 Hz, and their harmonics (kf1,kf2(k=1,2,3,…)). Some sidebands appear around the main frequency components, but most of them are irregular and the amplitudes are small. So this does not indicate that the gearbox is damaged.Figure 16
Waveform and corresponding Fourier spectrum of the baseline signal collected from the new SD16 planetary gearbox.
(a)
(b)


## 4.3. Damaged Sun Gear Signal Analysis

Figure17 shows the waveform and Fourier spectrum of the signal sample (the sample length is 1 s) from the same gearbox with naturally broken sun gear s1, respectively. Compared with Figure 16(a), the kurtosis indicator of the signal sample increases from 3.01 to 4.04. Several impulses are visible in the time domain, but the accurate impulse period is difficult to evaluate due to strong back noise. The dominant spectral peaks appear at the frequency components: the fundamental meshing frequency of the first stage (f1), the fundamental meshing frequency of the second stage (f2), their harmonics (kf1,kf2(k=1,2,3,…)), and the sidebands, as the baseline spectrum, but the magnitude of the spectral peaks increases significantly. All the signs indicate the anomalies of the gearbox. As described in the simulation, the sidebands are associated with the transient impulses and contain useful information of damage location. So the zoomed-in Fourier spectrum is shown in Figure 17(c). It can be observed that the sidebands appear around the meshing frequencies of the first stage, at the frequencies of kf1±mfs1 (k=1,2,…,m=1,2,…). To further validate the fault features, the spectral kurtosis signal processing method [23] is employed in this section to extract the transient impulses. Figure 18 presents the spectral kurtosis for the signal sample (the sample length is 5 s). The maximum spectral kurtosis appears at level 5.5 for the central frequency of 800 Hz. The filtered fault signal from the optimal band-pass filter around 800 Hz and the corresponding squared envelope spectrum of the filtered signal are shown in Figure 19. The periodic impulses are obvious in time domain, and a pattern of 0.062 s is confirmed. Frequency of 16.23 Hz and its harmonics exist in the frequency domain, which correspond with the calculation of (16) for the fault characteristic frequency of the sun gear in the first stage mfs1  (m=1,2,3,…). So we can conclude that the gear s1 has fault. The conclusion matches the results obtained from the simulation.Figure 17
Vibration signal from the gearbox with naturally broken sun gears1: (a) waveform, (b) Fourier spectrum, and (c) zoomed-in Fourier spectrum.
(a)
(b)
(c)Figure 18
Spectral kurtosis for the signal sample in Figure17.Figure 19
(a) Filtered fault signal from the optimal band-pass filter in the case of the damaged sun gears1 and (b) squared envelope spectrum of the filtered signal.
(a)
(b)Figure20 shows the waveform and Fourier spectrum of the signal sample (the sample length is 1 s) from the broken sun gear s2 case. The kurtosis value of the signal sample increases from 3.01 to 5.24. More sidebands appear in the frequency band from 500 to 1000 Hz, and most of them have bigger amplitudes than the baseline signal. As shown in Figure 20(c), the sidebands center around the meshing frequencies of the second stage, at the frequencies of kf2±mfs2 (k=1,2,…,m=1,2,…). The corresponding spectral kurtosis is presented in Figure 21. The maximum spectral kurtosis appears at level 6.5 for the central frequency of 80 Hz. The filtered fault signal from the optimal band-pass filter around 80 Hz and the corresponding squared envelope spectrum of the filtered signal are shown in Figure 22. The periodic impulses are evident, and a pattern of 0.056 s is confirmed. Frequency of 17.67 Hz and its harmonics exist in the frequency domain, which correspond with the calculation of (16) for the fault characteristic frequency of the sun gear in the second stage mfs2  (m=1,2,3,…). All the signs indicate that the gearbox is fault, and the damage occurs on the sun gear of the second stage. The results match well with the simulation.Figure 20
Vibration signal from the gearbox with naturally broken sun gears2: (a) waveform, (b) Fourier spectrum, and (c) zoomed-in Fourier spectrum.
(a)
(b)
(c)Figure 21
Spectral kurtosis for the signal sample in Figure20.Figure 22
(a) Filtered fault signal from the optimal band-pass filter in the case of the damaged sun gears2 and (b) squared envelope spectrum of the filtered signal.
(a)
(b)


## 5. Conclusions

Model-based dynamic analysis has been an important application for identifying and extracting useful fault features from vibration signals. In this study, the nonlinear dynamic model of a two-stage planetary gear set is established. Based on the improved potential energy method, the time-varying mesh stiffness of each mesh pair is calculated, and the influence of crack propagation on the mesh stiffness of sun-planet mesh pair is presented. Then, the simulated signals of the gear set are analyzed in time and frequency domain. According to the sidebands structure, the equations of fault characteristic frequencies of sun gear components are derived. These results provide a theoretical reference for diagnosing faults via modulation sidebands analysis in practice. Further, the experimental signals collected from a SD16 planetary gearbox are analyzed, and the spectral kurtosis is employed to extract the fault characteristic frequencies of sun gear from the processed signals. In future research, more suitable signal processing methods will be developed to improve the vibration test signal from the field, and more useful fault features will be extracted to accurately locate the gear damage and reliably assess the damage status.



---

*Source: 9263298-2015-12-27.xml*