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.
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 Section 1, 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
Figure 1 shows the structure of a two-stage planetary gear set. All the gears are standard involute spur tooth. The first stage () is a simple planetary gear set with planets (). The second stage () is a meshed planet planetary gear set with planet-planet mesh pairs (). The planets of each stage distribute equidistantly and have the same parameters. The ring gear and the carrier are shared by both stages. In this structure, the sun gear is fixed. The sun gear and carrier are identified as input member and output member. and are the input motor torque and output load torque.
The dynamic model of this two-stage planetary gear set is shown in Figure 2. 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 and viscous damping coefficient 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 and direction, denoted by and , and rotation around the axis ( represents the center components , , , and and the planets ).
(a)
(b)
2.1. Equations of Motion
Define as the rotational coordinates, where ( represents the center components , , and and the planets ) are the base radii of the gears. According to the meshing relation, the relative gear mesh displacements along the line of action of stage are expressed as follows:where , for internal mesh pair and external mesh pair, and for internal mesh pair and external and mesh pairs. is the spacing angle of the planets. is the pressure angle.
The dynamic gear mesh force is defined as follows:the bearing forces are defined as follows:where is the mesh damping coefficient, is the periodically time-varying mesh stiffness, is the damping coefficient in radial direction, is the radial stiffness of the bearing, is the damping coefficient in torsional direction, and is the torsional stiffness of the bearing.
The equations of motion for sun gear can be derived as follows:The equations of motion for ring gear are written as follows:The equations of motion for planets are given as follows:where and represent the mass and mass moment of inertia of the gears.
For the carrier-planet subsystem, is defined as the coordinate in place of , where is the radius from the carrier center to the planet center. The relative displacements between carrier and planet in and directions are given as follows:The differential equations of motion for carrier and planets are expressed as follows:where and 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:where In these equations, is the inertia matrix, is the displacement matrix, is the mesh damping matrix, is the bearing damping matrix, is the periodically time-varying mesh stiffness matrix, the submatrices of are given in Appendix, is the bearing stiffness matrix, and 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 Table 1.
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.
Thus, for the single-tooth pair meshing duration, the total mesh stiffness can be calculated as follows:where subscripts 1 and 2 represent the driving gear and driven gear. , , , and 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:where represents the first pair of meshing teeth and 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 and are set to be 0. and denote the theoretical meshing line and practical action line. is the pitch point of each mesh pair. As shown in Figure 1, the sun gear is fixed and is defined as the input member. Assuming the sun gear rotates in the counter-clockwise direction, the planets , , and are in clockwise, counter-clockwise, and counter-clockwise rotation directions, respectively.
(a)
(b)
Table 2 shows the relative phases of this gear set. The superscripts represent the first stage and the second stage of the gear set, respectively. and denote the relative phase between the th sun-planet mesh pair and the first sun-planet mesh pair. and are the relative phase between the th ring-planet mesh pair and the first ring-planet mesh pair. is the relative phase between the th planet-planet mesh pair and the first planet-planet mesh pair. and are the relative phase between the th ring-planet mesh pair and sun-planet mesh pair. is the relative phase between the th planet-planet mesh pair and sun-planet mesh pair. “+” dictate the phase lag, and “−” dictate the phase lead.
We use MATLAB programs to obtain the numerical values of total mesh stiffness as a function of time , 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. indicates the th 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 th planet-planet mesh pair and sun-planet mesh pair is in time, where is the mesh period of the second stage. Similarly, the phase difference between the th ring-planet mesh pair and sun-planet mesh pair can be expressed as in time.
(a)
(b)
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 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 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 .
The total mesh stiffness curves of the sun-planet pair in second stage with different crack lengths are shown in Figure 8. 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.
The mesh damping coefficient is set to be proportional to the total mesh stiffness; that is,where is the scale constant measured in seconds. For this gear set, s, s, s, and 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 Table 1, 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 gears and 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 is equal to 0.00606 s, while the impulse duration of 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 , the fundamental meshing frequency of second stage , and their harmonics , in which the frequency component 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:in this equation, is the number of sun gear teeth and is the rotation speed of sun gear relative to the carrier.
(a)
(b)
Supposing a 70% level (, ) crack appears on the sun gear of the first stage, the dynamic response of sun gear 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 , at the frequencies of , and is equal to 16.8 Hz. So the side frequency components indicate that the gear set have fault, and the fault exists on the sun gear .
(a)
(b)
(c)
Similarly, when the 70% level (, ) 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 s, as shown in Figure 12(a). The Fourier spectrum of the gear set with cracked sun gear , as shown in Figure 12(b), has almost the same frequency components with Figure 11(b): the fundamental meshing frequency of the first stage , the fundamental meshing frequency of the second stage , their harmonics , 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 , and is 18.2 Hz. As expected, the side frequency components contain the fault information, and the fault exists on the sun gear .
(a)
(b)
(c)
As described in Section 2, 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.
(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: According to (14), the equation can also be written as follows: where 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 gear . The third one is the signal collected from the same gearbox with naturally broken sun gear .
4.1. Experimental Specification
The planetary gearbox test rig is shown in Figure 14 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.
(a)
(b)
(a)
(b)
(c)
For the damaged sun gear experiment, the perfect sun gear of the first stage () 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 are in simulation, these useful fault features are expected to be detected. For the damaged sun gear experiment, the perfect sun gear of the second stage () is replaced by a broken one, whereas all the other parts are normal. The fault characteristic frequencies of , 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 Figure 16. The measured spectrum contains the fundamental meshing frequency of the first stage (), about 167.3 Hz, the fundamental meshing frequency of the second stage (), about 180.6 Hz, and their harmonics (). 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.
(a)
(b)
4.3. Damaged Sun Gear Signal Analysis
Figure 17 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 , 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 (), the fundamental meshing frequency of the second stage (), their harmonics (), 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 (). 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 (). So we can conclude that the gear has fault. The conclusion matches the results obtained from the simulation.
(a)
(b)
(c)
(a)
(b)
Figure 20 shows the waveform and Fourier spectrum of the signal sample (the sample length is 1 s) from the broken sun gear 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 (). 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 (). 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.
(a)
(b)
(c)
(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.
Appendix
Consider the following:
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
This research is supported by the National High Technology Research and Development Program of China (no. 2013AA040204) and the National Key Basic Research Development Plan of China (no. 2011CB013403). The authors thank Dr. X H. Liang for the help in validating the analytical mesh stiffness results. Comments and suggestions from anonymous referees and editors are quit valuable in the improvement of the quality of this paper.