# Predictive Modeling of a Two-Stage Gearbox towards Fault Detection

Authors: Edward J. Diehl; J. Tang
Journal: Shock and Vibration
 (2016)

---

## Abstract

This paper presents a systematic approach to the modeling and analysis of a benchmark two-stage gearbox test bed to characterize gear fault signatures when processed with harmonic wavelet transform (HWT) analysis. The eventual goal of condition monitoring is to be able to interpret vibration signals from nonstationary machinery in order to identify the type and severity of gear damage. To advance towards this goal, a lumped-parameter model that can be analyzed efficiently is developed which characterizes the gearbox vibratory response at the system level. The model parameters are identified through correlated numerical and experimental investigations. The model fidelity is validated first by spectrum analysis, using constant speed experimental data, and secondly by HWT analysis, using nonstationary experimental data. Model prediction and experimental data are compared for healthy gear operation and a seeded fault gear with a missing tooth. The comparison confirms that both the frequency content and the predicted, relative response magnitudes match with physical measurements. The research demonstrates that the modeling method in combination with the HWT data analysis has the potential for facilitating successful fault detection and diagnosis for gearbox systems.


---

## Body

## 1. Introduction

Gears, commonly used to transmit torque, are critical components of many energy and power, aircraft, automobile, and marine/ship systems. The topic of gear-induced vibration has been widely studied and remains an important, fundamental topic in engineering. In such systems, tooth-meshing induces vibrations, and primary and sideband tones arise which generally span a very wide frequency range. While excessive vibrations cause reliability and durability issues, vibration signals in a gearbox are frequently used as information carriers in gear system condition monitoring and fault detection/diagnosis [1, 2].The variety of monitoring and diagnosis approaches can be loosely categorized into data-based ones, model-based ones, or a hybrid of these two [3]. In data-based approaches, signal processing efforts are employed to pursue faults by analyzing measured operational data to identify distinguishable features, indicating faults based on patterns learned from empirical experience. Fault detection based on signal processing alone can be characterized as a “black box” approach; as the specific physics of the system are not explicitly considered, only matching of processed data with known fault patterns is involved. One major issue of this type of approaches is the difficulty in extrapolating from existing databases. Data-based condition monitoring approaches of gearboxes, specifically, process vibration signals which contain a complex mixture of harmonic and subharmonic components of mechanically induced motion from gear meshing, bearing rolling, and driver and load motion and random noise. Much of the important information needed for condition monitoring is hidden by the undesired components within the signal, but decomposing a measured signal by filtering often leads to removal of desired content [4]. Other drawbacks of data-based condition monitoring approaches include their lack of physical interpretation and their specificity to a particular piece of machinery, requiring machinery experience for each kind of fault in order to distinguish the specific fault.Model-based approaches to condition monitoring, often referred to as “white-box” approaches, can in theory provide baseline data by subjecting the model to the same operating condition of the actual system, and fault detection is facilitated by comparing the model prediction with the measurement data acquired by the monitoring system. Typically, dynamic gearbox models (DGMs) are developed with lumped-parameters at the system level and can yield model prediction results within reasonable computational time [4–7]. DGMs are meant to mimic the underlying physics within a gear system and simulate the vibration response of the physical unit, but a well-known weakness of most modeling efforts is the difficulty in achieving precise, or even acceptable, correlation with measured data, since DGMs only simulate the deterministic dynamic behavior. Nonplanetary multistage models have been investigated [8–13]. Apparently, existing DGM literature shows a relatively limited amount of experimental correlation for model validation. Many DGM studies rely upon the experimental data obtained by Munro [14], Kubo et al. [15], and Blankenship and Kahraman [16]. Some other studies include their own test beds for model validation [17–23]. All these investigations use single stage gearboxes operating at multiple constant speeds to demonstrate the “jump” phenomenon at natural frequencies when accelerating or decelerating to each constant speed. Experimental test data from multistage gear test beds is rare in the literature search, as are multistage DGMs, and no validated models are found.Condition monitoring of rotating machinery has historically been performed on machines operating at constant speeds by processing vibration data with the Fast Fourier Transform (FFT) into the frequency domain. Nevertheless, FFT and other frequency domain techniques are not well suited for analyzing nonstationary signals such as machinery with changing speeds [24]. This complication can be overcome by remaining in the time domain or by applying time-frequency techniques [1, 3, 25]. Time-frequency domain analysis reveals the waveform energy distribution across both time and frequency and can be used to identify transient events as well as fault patterns. Examples of time-frequency distribution techniques are the short-time Fourier transform (STFT), Wigner-Ville distribution (WVD), and wavelet transforms (WT). The STFT segments a signal into discrete time windows to find the frequency content in that window, but window size limits the resolution, making STFT only effective for analysis on slowly changing nonstationary signals [1]. While WVD does not segment time as STFT, it is susceptible to interference between “cross terms” that can cause misleading results [24]. WT is conceptually similar to FFT except that localized functions, called “mother wavelets,” are used instead of sinusoidal functions, allowing the time of an event to be retained. WT is used in numerous signal processing applications including nonstationary vibration signals [26–28]. The harmonic wavelet transform (HWT) has shown promise in identifying faults in nonstationary gear vibration signals [29].In this research, we use a benchmark two-stage gearbox test bed with variable speed controller to gather nonstationary experimental data and compare to DGM results. The DGM proposed includes translational degrees-of-freedom (DOFs) at locations and in directions corresponding to the accelerometers placed on the bearings of the test bed for improved experimental correlation. DGM parametric investigations and test bed experimental investigations are performed in parallel in order to refine and validate the model. Mimicking the shaft speed-time profile of monitored equipment with a DGM will allow a wider range of collected vibration data for building fault pattern libraries. To facilitate effective data processing especially under nonstationary conditions, the harmonic wavelet transform is employed to analyze the features of the vibration responses. It is identified that the combination of the proposed DGM and the HWT analysis yields effective correlation of vibration features with and without faults, showing promise for feature extraction and fault identification in nonstationary equipment.


## 2. Test Bed Structure and Experimental Setup

The test bed used in this research, shown in Figure1, includes a 3-HP motor driving a constant torque magnetic brake load through a two-stage gearbox [30]. The test bed is equipped with a variable speed controller and outfitted with multiple accelerometers and a digital tachometer.Figure 1
Gearbox test bed.Accelerometers are positioned on the intermediate shaft bearing hub in vertical, horizontal, and axial orientations for data collection and on the motor for data collection triggering. The accelerations and tachometer data are converted to digital signals at 20 kHz and recorded through a dSPACE data acquisition system with an ADC board. The speed of the input shaft is computer-controlled through a USB connection to the variable speed controller and digital tachometer. The representative speed-time profile input for this investigation is shown in Figure2.Figure 2
Speed-time profile input.The gearbox transmission, as shown in Figure3, is made up of three 32-tooth pinions on the input shaft, one of which drives an 80-tooth gear on the intermediate shaft connected to a 48-tooth pinion driving a 64-tooth gear on the output shaft. The 80-tooth gear can be manually shifted along the intermediate shaft to align with any of the three 32-tooth pinions. This arrangement allows one healthy and two faulty pinions to be investigated without complete removal of the shaft with each trial.Figure 3
Gearbox test bed arrangement.


## 3. Dynamic Gearbox Model (DGM) Formulation

### 3.1. Model Structure

The number of degrees-of-freedom (DOFs) representing a two-stage gear set in a DGM can range from a minimum of three (in which the angular displacement of each gear is modeled) to very large finite element models that attempt to incorporate all the local details. In general, a model must be sufficiently complex to capture the desired physical features while being small-scale enough for practical computational solving speeds. Previously Diehl et al. [31] investigated 3-, 6- and 12-DOF DGMs. In this research a 26-DOF model is formulated, which includes translational DOF representations of bearings at each end of each shaft instead of a single bearing per shaft and translational DOF at each gear center in both the vertical and horizontal directions. These additional DOFs allow shaft bending stiffness to be represented. The DOFs at the bearings can more closely represent the measured output of horizontally and vertically oriented accelerometers placed on the bearing housing of the gearbox test bed. Another benefit of the increased DOFs is the ability to represent shaft displacement at the gears for comparison to proximity probes installed in the gearbox.The 26-DOF DGM schematically shown in Figure4 is mathematically represented by (1)Iiθ¨i=α·Ti-WjRi+ctiθ˙i-θ˙i+1+ktiθi-θi+1,(2)mnx¨n=cx˙n+kxn+kbxn-xn+1-kbxn+1-xn+2+Wj·cos⁡β-ϕ,(3)mny¨n=cy˙n+kyn+kbyn-yn+1-kbyn+1-xn+2+Wj·sin⁡β-ϕ,(4)W1=km1θ1·DTE1+cm1R1θ˙1-R2θ˙2+cos⁡15.5km1x2-x5+cm1x˙2-x˙5-sin⁡15.5km1y2-y5+cm1y˙2-y˙5,(5)W2=km2θ3·DTE2+cm2R3θ˙3-R4θ˙5+cos⁡15.5km2x6-x9+cm2x˙6-x˙9+sin⁡15.5km2y6-y9+cm2y˙6-y˙9,where (1) characterizes the 6 rotational DOFs and (2) and (3) characterize the 20 translational DOFs. The complete expressions are provided in the Appendix. This system of equations is coupled by their in-common displacements, x
n and y
n, as well as the dynamic gear mesh force of each gear pair, W
j, represented by (4) and (5), that are functions of dynamic transmission error (DTE) to be discussed later. Other variable definitions are provided in Nomenclature.Figure 4
26-DOF DGM diagram.Combining (1) through (5), we obtain(6)Mz¨+Cz˙+Kz=F,where(7)M=I100⋯00I20⋯000I3⋯0⋮⋮⋮⋱00000m10,(8)K=kt1-kt10⋯0-kt1kt1+km1R12-km1R1R2⋯00-km1R1R2kt2+km1R22⋯0⋮⋮⋮⋱00000k+kb3,(9)C=ct1-ct10⋯0-ct1ct1+cm1R12-cm1R1R2⋯00-cm1R1R2ct2+cm1R22⋯0⋮⋮⋮⋱00000c+cb3,(10)F=TD0⋯TL⋯0T.Controlling the speed of the DGM is achieved by varying the amount of drive torque,T
D, in (1) using a multiplier referred to as torque factor, α. The load, T
L, of the physical system is a magnetic brake providing a constant torque load, so the drive torque of the system operating at steady speeds is related to the load torque by the overall gear ratio as shown below:(11)TD=αθDN1N3N2N4TL,where N
1
,
N
2
,
N
3, and N
4 are the number of gear teeth 32, 80, 48, and 64, respectively, for the specific gearbox test bed studied. Values of the torque factor greater than unity accelerate the system while values below decelerate and are controlled using a call-out function based on drive position, θ
D. The call-out function changes the torque factor at drive positions, determined by the angular kinematics, so a known speed-time profile can be replicated.The mesh stiffness of each gear pair varies with angular position since the number of teeth in contact of the gear pair is not constant. Additionally, within each gear pair, the point of contact moves along the involute profile, from tip to root on one tooth and root to tip on the other, the individual tooth stiffness varying with position. The contact ratio, the average quantity of teeth engaged during a mesh cycle, usually varies between 1.6 and 2.5 and is useful in describing the variation in mesh stiffness along a gear pair’s line-of-action (LOA) when treated like a linear spring on the LOA [16]. Figure 5 shows the mesh stiffness, k
m
1 and k
m
2, of each gear pair used in this model, implemented as separate functions of each pinion position, θ
1 and θ
3, respectively. These mesh stiffness profiles are developed from finite-element-analysis (FEA) models of each entire gear with unit loads placed at discrete positions along a single tooth face. To create gear pair mesh stiffness, the individual teeth stiffness degrees are combined based on contact position and assembled as springs in series for teeth pairs or parallel for quantity of teeth pairs in contact. The gear contact ratio is confirmed through this approach, revealing a peak when three teeth are in contact and an arc of lower values when only two teeth are in contact. A noteworthy assumption made in this calculation is the omission of the contribution of contact stiffness, which is dependent on load. The magnitudes of the two teeth pairs and three teeth pairs are compared with those found using equations developed by Kiekbusch et al. [32] and Kuang and Yang [33]. These comparisons show adequate agreement but most notably that the load used in this paper resulted in a relatively insignificant contribution by the contact stiffness. Most gear faults can be modeled by altering the mesh stiffness profile as will be demonstrated.Figure 5
Mesh stiffness along the line-of-action (LOA) versus pinion angle for first- and second-stage gear sets.The dynamic transmission error (DTE) is defined as, for each respective gear pair,(12)DTE1=R1θ1-R2θ2-bR1θ1-R2θ2>b,0b≤R1θ1-R2θ2≤b,R1θ1-R2θ2+bR1θ1-R2θ2<b,DTE2=R3θ3-R4θ4-bR3θ3-R4θ4>b,0b≤R3θ3-R4θ4≤b,R3θ3-R4θ4+bR3θ3-R4θ4<b.They become a loss-of-contact nonlinear function with the inclusion of backlash, 2
b, the clearance between gear pair teeth along the line-of-action [14]. This model feature has been addressed in the past and shown to be the source of nonlinear behavior in gearboxes often related to whether the units are accelerating up to or decelerating down to the constant speeds [34].Using the stiffness matrix (see (8)) and the input torque and load torque (see (10)) as a static problem, the initial displacements are found to avoid exaggerated large initial accelerations. The mass matrix (see (7)) and stiffness matrix (see (8)) are also used to determine the natural frequencies of the system by finding the eigenvalues, which are used to explore the effects of parameter variation on the system. With the initial displacements from the static case and the initial velocities set at zero, the model is solved using Runge Kutta time marching with time steps matching the experimental data sampling frequency.


### 3.2. Parametric Identification

The load torque for the DGM matches the specific tests, 10.3 Nm for the constant speed and 6.6 Nm for varying speed, which are determined by the magnetic brake manufacturer’s torque versus input current curves and are verified using a simple scale and lever arm. All four gears have a 1.59 mm module (16 teeth per inch diametral pitch) and a 14.5° pressure angle resulting. The pitch radii are 25.4 mm, 63.5 mm, 38.1 mm, and 50.8 mm for variablesR
1
,
R
2
,
R
3, and R
4, respectively. The resulting contact ratios are 2.12 and 2.16 for gear pairs one and two, respectively. The backlash used is 0.5 mm. The DGM requires stiffness, damping, and inertia parameters to be gathered, calculated, measured, or estimated for each component associated with a DOF. Table 1 summarizes the parameter values identified for the specific gearbox test bed. The masses of the shafts and gears were measured directly, while the moments of inertia were found from Solidworks models. The bending and torsional stiffness of the shafts were extracted from FEA models. The mesh stiffness, as previously discussed, was calculated based on FEA stiffness results. Other parameters, such as bearing stiffness and the damping parameters, were estimated using a combination of techniques, described below.Table 1
Parameters used in dynamic gear model.
Degree-of-freedom
Mass
Stiffness
Damping
Variable
Value
Units
Variable
Value
Units
Variable
Value
Units
1
I
D
7.51 × 10−3
kg m2
k
t
1
247
N m/rad
c
t
1
0.0926
N s/rad
2
I
1
3.86 × 10−4
kg m2
3
I
2
2.69 × 10−3
kg m2
k
t
2
39,900
N m/rad
c
t
2
1.45
N s/rad
4
I
3
4.18 × 10−4
kg m2
5
I
4
1.21 × 10−3
kg m2
k
t
3
12,600
N m/rad
c
t
3
0.507
N s/rad
6
I
L
6.64 × 10−5
kg m2
7 and 8
m
1
0.743
kg
k
3.50 × 107
N/m
c
831
N s/m
9 and 10
m
2
0.748
kg
k
b
1
1.83 × 108
N/m
c
b
1
0
N s/m
11 and 12
m
3
0.643
kg
k
3.50 × 107
N/m
c
831
N s/m
13 and 14
m
4
0.443
kg
k
3.50 × 107
N/m
c
831
N s/m
15 and 16
m
5
1.466
kg
k
b
2
a
2.43 × 108
N/m
c
b
2
a
0
N s/m
17 and 18
m
6
0.550
kg
k
b
2
1.32 × 108
N/m
c
b
2
0
N s/m
k
b
2
b
2.43 × 108
N/m
c
b
2
b
0
N s/m
19 and 20
m
7
0.443
kg
k
3.50 × 107
N/m
c
831
N s/m
21 and 22
m
8
0.651
kg
k
3.50 × 107
N/m
c
831
N s/m
23 and 24
m
9
0.942
kg
k
b
3
3.49 × 108
N/m
c
b
3
0
N s/m
25 and 26
m
10
0.648
kg
k
3.50 × 107
N/m
c
831
N s/mAccurately representing the dynamic behavior of bearings in the DGM is particularly challenging due to their stiffness nonlinearity with respect to preload, radial load, and rotational speed. Many papers have discussed this challenge [35–38]. While modeling of bearing behavior continues to advance, it is most advantageous for solution speed in a DGM to assume linear behavior by using a single value for stiffness. A series of investigations were performed to identify the most accurate single value for bearing stiffness. These investigations included comparing the load and displacement data supplied by the manufacturer (Rexnord MB) against published estimation techniques, parametric sensitivity of the DGM to bearing stiffness and comparison to constant speed experimental data, and hammer test frequency response analysis comparison to FEA modal results where bearing stiffness was included between shaft and casing. Damping was estimated from the log decrement of hammer tests on the in situ intermediate shaft, and the resulting damping ratio from this test was applied to each DOF, treating damping as being proportional to the associated inertia and stiffness values. Other assumptions, besides proportional damping and bearing stiffness linearity, include shaft mass lumped at the bearings and gears and ignoring shaft gyroscopic resonance.Features included in other DGMs but omitted in this model include manufacturing indexing errors [39], gear profile and runout errors [40], gear teeth friction [41], and the influence of the gearbox housing [42]. These omitted features will be included in future research.


### 3.3. Fault Modeling

A variety of gear faults can be simulated in a DGM by altering the time-varying mesh stiffness to reflect the reduction caused by the defect. Many previous studies have contributed to the development of fault models such as tooth spalls and cracks [43], gear wear [44], and tooth chips and breakage [45]. A missing tooth fault, represented in Figure 6(a), is analyzed in this paper as an illustrative example, and other localized faults will be studied in future work. The mesh stiffness pattern is reduced as shown in Figure 6(b). For this high contact ratio gear pair, a missing tooth results in two short durations when only one pair of teeth is in contact interrupted by an even shorter period when two teeth pairs are in contact. The faulty mesh stiffness profile is developed similarly to the healthy meshing pattern, except that the stiffness contribution of one tooth is omitted for the duration of its engagement. This distinctive pattern will influence the acceleration at the model bearing as well as the accelerometer signal in the physical unit the model represents.Figure 6
Missing tooth fault (a) model representation and (b) stiffness profile of first-stage gear.
(a)
(b)


## 3.1. Model Structure

The number of degrees-of-freedom (DOFs) representing a two-stage gear set in a DGM can range from a minimum of three (in which the angular displacement of each gear is modeled) to very large finite element models that attempt to incorporate all the local details. In general, a model must be sufficiently complex to capture the desired physical features while being small-scale enough for practical computational solving speeds. Previously Diehl et al. [31] investigated 3-, 6- and 12-DOF DGMs. In this research a 26-DOF model is formulated, which includes translational DOF representations of bearings at each end of each shaft instead of a single bearing per shaft and translational DOF at each gear center in both the vertical and horizontal directions. These additional DOFs allow shaft bending stiffness to be represented. The DOFs at the bearings can more closely represent the measured output of horizontally and vertically oriented accelerometers placed on the bearing housing of the gearbox test bed. Another benefit of the increased DOFs is the ability to represent shaft displacement at the gears for comparison to proximity probes installed in the gearbox.The 26-DOF DGM schematically shown in Figure4 is mathematically represented by (1)Iiθ¨i=α·Ti-WjRi+ctiθ˙i-θ˙i+1+ktiθi-θi+1,(2)mnx¨n=cx˙n+kxn+kbxn-xn+1-kbxn+1-xn+2+Wj·cos⁡β-ϕ,(3)mny¨n=cy˙n+kyn+kbyn-yn+1-kbyn+1-xn+2+Wj·sin⁡β-ϕ,(4)W1=km1θ1·DTE1+cm1R1θ˙1-R2θ˙2+cos⁡15.5km1x2-x5+cm1x˙2-x˙5-sin⁡15.5km1y2-y5+cm1y˙2-y˙5,(5)W2=km2θ3·DTE2+cm2R3θ˙3-R4θ˙5+cos⁡15.5km2x6-x9+cm2x˙6-x˙9+sin⁡15.5km2y6-y9+cm2y˙6-y˙9,where (1) characterizes the 6 rotational DOFs and (2) and (3) characterize the 20 translational DOFs. The complete expressions are provided in the Appendix. This system of equations is coupled by their in-common displacements, x
n and y
n, as well as the dynamic gear mesh force of each gear pair, W
j, represented by (4) and (5), that are functions of dynamic transmission error (DTE) to be discussed later. Other variable definitions are provided in Nomenclature.Figure 4
26-DOF DGM diagram.Combining (1) through (5), we obtain(6)Mz¨+Cz˙+Kz=F,where(7)M=I100⋯00I20⋯000I3⋯0⋮⋮⋮⋱00000m10,(8)K=kt1-kt10⋯0-kt1kt1+km1R12-km1R1R2⋯00-km1R1R2kt2+km1R22⋯0⋮⋮⋮⋱00000k+kb3,(9)C=ct1-ct10⋯0-ct1ct1+cm1R12-cm1R1R2⋯00-cm1R1R2ct2+cm1R22⋯0⋮⋮⋮⋱00000c+cb3,(10)F=TD0⋯TL⋯0T.Controlling the speed of the DGM is achieved by varying the amount of drive torque,T
D, in (1) using a multiplier referred to as torque factor, α. The load, T
L, of the physical system is a magnetic brake providing a constant torque load, so the drive torque of the system operating at steady speeds is related to the load torque by the overall gear ratio as shown below:(11)TD=αθDN1N3N2N4TL,where N
1
,
N
2
,
N
3, and N
4 are the number of gear teeth 32, 80, 48, and 64, respectively, for the specific gearbox test bed studied. Values of the torque factor greater than unity accelerate the system while values below decelerate and are controlled using a call-out function based on drive position, θ
D. The call-out function changes the torque factor at drive positions, determined by the angular kinematics, so a known speed-time profile can be replicated.The mesh stiffness of each gear pair varies with angular position since the number of teeth in contact of the gear pair is not constant. Additionally, within each gear pair, the point of contact moves along the involute profile, from tip to root on one tooth and root to tip on the other, the individual tooth stiffness varying with position. The contact ratio, the average quantity of teeth engaged during a mesh cycle, usually varies between 1.6 and 2.5 and is useful in describing the variation in mesh stiffness along a gear pair’s line-of-action (LOA) when treated like a linear spring on the LOA [16]. Figure 5 shows the mesh stiffness, k
m
1 and k
m
2, of each gear pair used in this model, implemented as separate functions of each pinion position, θ
1 and θ
3, respectively. These mesh stiffness profiles are developed from finite-element-analysis (FEA) models of each entire gear with unit loads placed at discrete positions along a single tooth face. To create gear pair mesh stiffness, the individual teeth stiffness degrees are combined based on contact position and assembled as springs in series for teeth pairs or parallel for quantity of teeth pairs in contact. The gear contact ratio is confirmed through this approach, revealing a peak when three teeth are in contact and an arc of lower values when only two teeth are in contact. A noteworthy assumption made in this calculation is the omission of the contribution of contact stiffness, which is dependent on load. The magnitudes of the two teeth pairs and three teeth pairs are compared with those found using equations developed by Kiekbusch et al. [32] and Kuang and Yang [33]. These comparisons show adequate agreement but most notably that the load used in this paper resulted in a relatively insignificant contribution by the contact stiffness. Most gear faults can be modeled by altering the mesh stiffness profile as will be demonstrated.Figure 5
Mesh stiffness along the line-of-action (LOA) versus pinion angle for first- and second-stage gear sets.The dynamic transmission error (DTE) is defined as, for each respective gear pair,(12)DTE1=R1θ1-R2θ2-bR1θ1-R2θ2>b,0b≤R1θ1-R2θ2≤b,R1θ1-R2θ2+bR1θ1-R2θ2<b,DTE2=R3θ3-R4θ4-bR3θ3-R4θ4>b,0b≤R3θ3-R4θ4≤b,R3θ3-R4θ4+bR3θ3-R4θ4<b.They become a loss-of-contact nonlinear function with the inclusion of backlash, 2
b, the clearance between gear pair teeth along the line-of-action [14]. This model feature has been addressed in the past and shown to be the source of nonlinear behavior in gearboxes often related to whether the units are accelerating up to or decelerating down to the constant speeds [34].Using the stiffness matrix (see (8)) and the input torque and load torque (see (10)) as a static problem, the initial displacements are found to avoid exaggerated large initial accelerations. The mass matrix (see (7)) and stiffness matrix (see (8)) are also used to determine the natural frequencies of the system by finding the eigenvalues, which are used to explore the effects of parameter variation on the system. With the initial displacements from the static case and the initial velocities set at zero, the model is solved using Runge Kutta time marching with time steps matching the experimental data sampling frequency.


## 3.2. Parametric Identification

The load torque for the DGM matches the specific tests, 10.3 Nm for the constant speed and 6.6 Nm for varying speed, which are determined by the magnetic brake manufacturer’s torque versus input current curves and are verified using a simple scale and lever arm. All four gears have a 1.59 mm module (16 teeth per inch diametral pitch) and a 14.5° pressure angle resulting. The pitch radii are 25.4 mm, 63.5 mm, 38.1 mm, and 50.8 mm for variablesR
1
,
R
2
,
R
3, and R
4, respectively. The resulting contact ratios are 2.12 and 2.16 for gear pairs one and two, respectively. The backlash used is 0.5 mm. The DGM requires stiffness, damping, and inertia parameters to be gathered, calculated, measured, or estimated for each component associated with a DOF. Table 1 summarizes the parameter values identified for the specific gearbox test bed. The masses of the shafts and gears were measured directly, while the moments of inertia were found from Solidworks models. The bending and torsional stiffness of the shafts were extracted from FEA models. The mesh stiffness, as previously discussed, was calculated based on FEA stiffness results. Other parameters, such as bearing stiffness and the damping parameters, were estimated using a combination of techniques, described below.Table 1
Parameters used in dynamic gear model.
Degree-of-freedom
Mass
Stiffness
Damping
Variable
Value
Units
Variable
Value
Units
Variable
Value
Units
1
I
D
7.51 × 10−3
kg m2
k
t
1
247
N m/rad
c
t
1
0.0926
N s/rad
2
I
1
3.86 × 10−4
kg m2
3
I
2
2.69 × 10−3
kg m2
k
t
2
39,900
N m/rad
c
t
2
1.45
N s/rad
4
I
3
4.18 × 10−4
kg m2
5
I
4
1.21 × 10−3
kg m2
k
t
3
12,600
N m/rad
c
t
3
0.507
N s/rad
6
I
L
6.64 × 10−5
kg m2
7 and 8
m
1
0.743
kg
k
3.50 × 107
N/m
c
831
N s/m
9 and 10
m
2
0.748
kg
k
b
1
1.83 × 108
N/m
c
b
1
0
N s/m
11 and 12
m
3
0.643
kg
k
3.50 × 107
N/m
c
831
N s/m
13 and 14
m
4
0.443
kg
k
3.50 × 107
N/m
c
831
N s/m
15 and 16
m
5
1.466
kg
k
b
2
a
2.43 × 108
N/m
c
b
2
a
0
N s/m
17 and 18
m
6
0.550
kg
k
b
2
1.32 × 108
N/m
c
b
2
0
N s/m
k
b
2
b
2.43 × 108
N/m
c
b
2
b
0
N s/m
19 and 20
m
7
0.443
kg
k
3.50 × 107
N/m
c
831
N s/m
21 and 22
m
8
0.651
kg
k
3.50 × 107
N/m
c
831
N s/m
23 and 24
m
9
0.942
kg
k
b
3
3.49 × 108
N/m
c
b
3
0
N s/m
25 and 26
m
10
0.648
kg
k
3.50 × 107
N/m
c
831
N s/mAccurately representing the dynamic behavior of bearings in the DGM is particularly challenging due to their stiffness nonlinearity with respect to preload, radial load, and rotational speed. Many papers have discussed this challenge [35–38]. While modeling of bearing behavior continues to advance, it is most advantageous for solution speed in a DGM to assume linear behavior by using a single value for stiffness. A series of investigations were performed to identify the most accurate single value for bearing stiffness. These investigations included comparing the load and displacement data supplied by the manufacturer (Rexnord MB) against published estimation techniques, parametric sensitivity of the DGM to bearing stiffness and comparison to constant speed experimental data, and hammer test frequency response analysis comparison to FEA modal results where bearing stiffness was included between shaft and casing. Damping was estimated from the log decrement of hammer tests on the in situ intermediate shaft, and the resulting damping ratio from this test was applied to each DOF, treating damping as being proportional to the associated inertia and stiffness values. Other assumptions, besides proportional damping and bearing stiffness linearity, include shaft mass lumped at the bearings and gears and ignoring shaft gyroscopic resonance.Features included in other DGMs but omitted in this model include manufacturing indexing errors [39], gear profile and runout errors [40], gear teeth friction [41], and the influence of the gearbox housing [42]. These omitted features will be included in future research.


## 3.3. Fault Modeling

A variety of gear faults can be simulated in a DGM by altering the time-varying mesh stiffness to reflect the reduction caused by the defect. Many previous studies have contributed to the development of fault models such as tooth spalls and cracks [43], gear wear [44], and tooth chips and breakage [45]. A missing tooth fault, represented in Figure 6(a), is analyzed in this paper as an illustrative example, and other localized faults will be studied in future work. The mesh stiffness pattern is reduced as shown in Figure 6(b). For this high contact ratio gear pair, a missing tooth results in two short durations when only one pair of teeth is in contact interrupted by an even shorter period when two teeth pairs are in contact. The faulty mesh stiffness profile is developed similarly to the healthy meshing pattern, except that the stiffness contribution of one tooth is omitted for the duration of its engagement. This distinctive pattern will influence the acceleration at the model bearing as well as the accelerometer signal in the physical unit the model represents.Figure 6
Missing tooth fault (a) model representation and (b) stiffness profile of first-stage gear.
(a)
(b)


## 4. Correlation of Predictive Modeling and Experimental Investigation

### 4.1. Model Validation with Constant Speed Experimental Comparisons

First, preliminary demonstration of the principles of damage detection and of the model capabilities to replicate many of the frequencies is made using constant speed operation of the gear test bed with a healthy pinion and a known faulty pinion. A model running with conditions matching the operating features and gear damage was run to simulate gearbox test data collected with a constant 820 r/min. Plots of accelerometer data and model acceleration versus time are shown in Figure7 for healthy and damaged gear conditions and arranged for comparison. While the acceleration magnitude range is not the same, the relative increases in magnitude due to the fault and the fault duration are quite similar. When the experimental data and model results are processed with a basic spectrum analysis (i.e., FFT), the frequencies expected to be pronounced are the shaft speed times number of teeth, the gear mesh frequency [46]. As such, the first-stage gear mesh frequency (GMF1) is expected to be 437 Hz and the second-stage gear mesh frequency (GMF2) 262 Hz. These expected frequencies and their integer multiples are listed in Table 2.Table 2
Expected frequencies (Hz) from gearbox operating at 820 r/min.
ΔGMFa
GMF2
GMF1
2 × GMF2
GMF1 + ΔGMF
3 × GMF2
2 × GMF1
4 × GMF2
3 × GMF1 and 5 × GMF2
4 × GMF1
175
262
437
525
612
787
875
1050
1312
1750
GMF: gear mesh frequency.
aThe difference between gear mesh frequencies of the first and second gear pairs.Figure 7
Experimental data and model results for 820 rpm constant speed test: (a) healthy data, (b) healthy model, (c) missing tooth data, and (d) missing tooth model.
(a)
(b)
(c)
(d)Sets of test data for this operating condition were recorded with a healthy gear set and a first-stage pinion missing a tooth. The healthy and faulty 26-DOF models were solved with torque and speed conditions matching the test. The spectrum results of unfiltered experimental data and model results are presented in Figure8. Peak frequencies are labeled and show many of the same frequency changes and features indicating gear damage.Figure 8
Spectrum analysis of experimental data and model results for 820 rpm constant speed test: (a) healthy data, (b) healthy model, (c) damaged (missing tooth) data, and (d) damaged model. ①: 1st-stage gear mesh frequency; ②: 2nd-stage gear mesh frequency.
(a)
(b)
(c)
(d)The two primary fault indicators in a gear system spectrum analysis are the increased magnitude at multiples of the gear mesh frequency and the sideband prevalence near the primary gear mesh frequencies and integer multiples [47]. The expected magnitude increases at the gear mesh frequencies and multiples, as well as the presence of sideband frequency content, are present in both test and model results. These frequency plots, while validating that the model captures many of the effects of gear damage, also reaffirm the difficulty achieving perfect agreement between the data and model. Many of the peaks are disproportionate, with the low frequency content being particularly underrepresented in the model, requiring a zoomed window to show this detail. In spite of these shortcomings, this data and model comparison demonstrates the ability of the DGM to capture the behavior of the physical gearbox and represent a damaged gear.


### 4.2. Data Processing Using Harmonic Wavelet Transform

Wavelets are useful in condition monitoring for their ability to characterize frequency and time content within a signal to find faults in rotating equipment as documented in [2, 48, 49]. In particular, the harmonic wavelet transform (HWT) [50] has several advantages over other wavelets including well-defined bandwidths and a structure similar to the natural occurring structure of vibration waves [51]. The HWT mother wavelet appears as a step function in the frequency domain, with a well-defined bandwidth, so when transformed from the frequency domain to the time domain, the resulting harmonic wavelet has little to no bandwidth leakage. This approach to mother wavelet creation is logical. Another advantage of the HWT over other wavelet transforms is the relative ease in its implementation, since it can be constructed with an algorithm which uses the Fast Fourier Transform (FFT) and Inverse Fast Fourier (IFF) Transform. The result of HWT is a matrix of coefficients, and for vibration data these coefficients represent the distribution of energy over time and frequency which are represented by contour plots.Because HWT are formed in the frequency domain, there is minimal leakage between frequency bands. An important decision needed to implement the HWT effectively is the choice of bandwidths (levels) that will efficiently balance resolution in time and frequency to represent the desired features. The divisions within a partition do not necessarily need to form a specific pattern [26–28], but a consistent partition between analyses is needed to make detailed comparisons. Every partition is a compromise between time resolution and frequency resolution, so a variety of patterns were explored for this present research. A partition that favors frequency to time by 4 : 1 was chosen for this present effort because while the experimental data sampling rate, 20 kHz, which also matches the DGM time steps, has a folding frequency of 10 kHz and therefore a range of 0–10 kHz, the range of interest was found to be typically 0–2 kHz. Processing a 4-second test with this HWT partition creates a 400 row, 100 column matrix of coefficients as seen in Figure 9.Figure 9
HWT of experimental data and model results for 820 rpm constant speed test: (a) healthy data, (b) healthy model, (c) damaged (missing tooth) data, and (d) damaged model.
(a)
(b)
(c)
(d)Indeed, Figure9 shows color contour plots, orwavelet maps, of the HWT coefficients for the constant speed tests previously shown in Figures 7 and 8. Wavelet maps indicate the energy distribution with frequency on the vertical axis and time on the horizontal axis. A contour color scale was chosen so white-to-blueish colors represent regions of low energy and reddish-to-black regions represent higher energy. The magnitude of the coefficients was adjusted to use the same contour scales for both processed data and model results. The fault’s influence on the vibration in both time and frequency is easily recognizable, especially in the DGM simulation. The disruption due to the missing tooth appears as a pulse in the HWT of the damaged gear at approximately the rate of input shaft rotation, 820 rev/min or 13.7 Hz, for a period of 0.073 seconds. This figure demonstrates the capability of HWT and DGM to characterize a fault in the signal, but it also highlights a major benefit of using a DGM to build fault libraries, as the fault signature is more easily distinguishable in the signal than it is in the damaged experimental data.


### 4.3. Comparison of Experimental Data and Model Results under Nonstationary Conditions

The time-frequency analysis capability is especially important for condition monitoring when machinery speed is not constant. One important example is the operation of a wind turbine whose rotational speed is, in part, dependent on wind speed. Data-driven condition monitoring approaches generally require periods of constant speed to build libraries and assess machinery, but DGM can mimic the nonstationary speed-time profile of the equipment and can be solved repeatedly with healthy gear and numerous simulated faults.To demonstrate the nonstationary capability, the 4-second motor speed-time profile used by [29] and shown in Figure 2 is used. The speed increases from stop to 10 Hz after one second and then increases up to 25 Hz at two seconds and maintains this speed to three seconds before decelerating to a stop. It is important to note that the speed-time profile described here represents the control message to the motor and that actual motor speed varies with each run. This variation can be mimicked by the DGM, but data-driven approaches are generally unable to account for this variation. Figure 10 presents a comparison of HWT processed test data with healthy gear and missing tooth and HWT processed model results with the respective measured speed-time profiles. Note that the speed-profile of the healthy test was used for the healthy DGM, and likewise the damaged gear test speed-profile was used for the damaged DGM. The measured motor speed multiplied by the number of teeth in the first-stage pinion is the anticipated gear mesh frequency (GMF1), and the second-stage pinion times the intermediate shaft speed is GMF2. These GMF frequency-time profiles and their integer multiples are the kinematic behavior of the system and are superimposed on the contour plots of Figure 10 for both gear stages. The energy distribution in the HWT plots matches many of these GMFs and multiples.Figure 10
HWT of experimental data and model results for varying speed test: (a) healthy data, (b) healthy model, (c) damaged (missing tooth) data, and (d) damaged model, with gear mesh frequencies (GMF) and integer harmonics superimposed.
(a)
(b)
(c)
(d)Visual comparisons of the HWT contour plots in Figure10 reveal obvious differences between the healthy and damaged cases in both experimental data and model results. While the healthy-damaged differences are not entirely consistent between data and model, there are fault feature trends in common, such as the increase in magnitude in GMF2 and widening of the regions near and around GMF1. While these contour plots aid in visualizing the wave map, the subtle numerical differences within the coefficient matrix contain feature patterns to be used to build libraries within pattern recognition algorithms.To demonstrate the effectiveness of this approach, Figure11 presents a contour plot of the subtraction of a damaged gear DGM HWT coefficient matrix from a healthy one when both use the same speed-profile. The ability to control the speed of the DGM precisely enables isolation of the effects of the damaged gear on the energy distribution in frequency and time. A noteworthy aspect is the similarity in the relatively constant speed region to the pulses found in the constant speed region of Figure 9. These pulses represent a potential feature that a pattern recognition algorithm can use when developing fault libraries.Figure 11
Contour plot of difference between HWT of healthy and damaged model results, with gear mesh frequencies (GMF) and integer harmonics superimposed.


## 4.1. Model Validation with Constant Speed Experimental Comparisons

First, preliminary demonstration of the principles of damage detection and of the model capabilities to replicate many of the frequencies is made using constant speed operation of the gear test bed with a healthy pinion and a known faulty pinion. A model running with conditions matching the operating features and gear damage was run to simulate gearbox test data collected with a constant 820 r/min. Plots of accelerometer data and model acceleration versus time are shown in Figure7 for healthy and damaged gear conditions and arranged for comparison. While the acceleration magnitude range is not the same, the relative increases in magnitude due to the fault and the fault duration are quite similar. When the experimental data and model results are processed with a basic spectrum analysis (i.e., FFT), the frequencies expected to be pronounced are the shaft speed times number of teeth, the gear mesh frequency [46]. As such, the first-stage gear mesh frequency (GMF1) is expected to be 437 Hz and the second-stage gear mesh frequency (GMF2) 262 Hz. These expected frequencies and their integer multiples are listed in Table 2.Table 2
Expected frequencies (Hz) from gearbox operating at 820 r/min.
ΔGMFa
GMF2
GMF1
2 × GMF2
GMF1 + ΔGMF
3 × GMF2
2 × GMF1
4 × GMF2
3 × GMF1 and 5 × GMF2
4 × GMF1
175
262
437
525
612
787
875
1050
1312
1750
GMF: gear mesh frequency.
aThe difference between gear mesh frequencies of the first and second gear pairs.Figure 7
Experimental data and model results for 820 rpm constant speed test: (a) healthy data, (b) healthy model, (c) missing tooth data, and (d) missing tooth model.
(a)
(b)
(c)
(d)Sets of test data for this operating condition were recorded with a healthy gear set and a first-stage pinion missing a tooth. The healthy and faulty 26-DOF models were solved with torque and speed conditions matching the test. The spectrum results of unfiltered experimental data and model results are presented in Figure8. Peak frequencies are labeled and show many of the same frequency changes and features indicating gear damage.Figure 8
Spectrum analysis of experimental data and model results for 820 rpm constant speed test: (a) healthy data, (b) healthy model, (c) damaged (missing tooth) data, and (d) damaged model. ①: 1st-stage gear mesh frequency; ②: 2nd-stage gear mesh frequency.
(a)
(b)
(c)
(d)The two primary fault indicators in a gear system spectrum analysis are the increased magnitude at multiples of the gear mesh frequency and the sideband prevalence near the primary gear mesh frequencies and integer multiples [47]. The expected magnitude increases at the gear mesh frequencies and multiples, as well as the presence of sideband frequency content, are present in both test and model results. These frequency plots, while validating that the model captures many of the effects of gear damage, also reaffirm the difficulty achieving perfect agreement between the data and model. Many of the peaks are disproportionate, with the low frequency content being particularly underrepresented in the model, requiring a zoomed window to show this detail. In spite of these shortcomings, this data and model comparison demonstrates the ability of the DGM to capture the behavior of the physical gearbox and represent a damaged gear.


## 4.2. Data Processing Using Harmonic Wavelet Transform

Wavelets are useful in condition monitoring for their ability to characterize frequency and time content within a signal to find faults in rotating equipment as documented in [2, 48, 49]. In particular, the harmonic wavelet transform (HWT) [50] has several advantages over other wavelets including well-defined bandwidths and a structure similar to the natural occurring structure of vibration waves [51]. The HWT mother wavelet appears as a step function in the frequency domain, with a well-defined bandwidth, so when transformed from the frequency domain to the time domain, the resulting harmonic wavelet has little to no bandwidth leakage. This approach to mother wavelet creation is logical. Another advantage of the HWT over other wavelet transforms is the relative ease in its implementation, since it can be constructed with an algorithm which uses the Fast Fourier Transform (FFT) and Inverse Fast Fourier (IFF) Transform. The result of HWT is a matrix of coefficients, and for vibration data these coefficients represent the distribution of energy over time and frequency which are represented by contour plots.Because HWT are formed in the frequency domain, there is minimal leakage between frequency bands. An important decision needed to implement the HWT effectively is the choice of bandwidths (levels) that will efficiently balance resolution in time and frequency to represent the desired features. The divisions within a partition do not necessarily need to form a specific pattern [26–28], but a consistent partition between analyses is needed to make detailed comparisons. Every partition is a compromise between time resolution and frequency resolution, so a variety of patterns were explored for this present research. A partition that favors frequency to time by 4 : 1 was chosen for this present effort because while the experimental data sampling rate, 20 kHz, which also matches the DGM time steps, has a folding frequency of 10 kHz and therefore a range of 0–10 kHz, the range of interest was found to be typically 0–2 kHz. Processing a 4-second test with this HWT partition creates a 400 row, 100 column matrix of coefficients as seen in Figure 9.Figure 9
HWT of experimental data and model results for 820 rpm constant speed test: (a) healthy data, (b) healthy model, (c) damaged (missing tooth) data, and (d) damaged model.
(a)
(b)
(c)
(d)Indeed, Figure9 shows color contour plots, orwavelet maps, of the HWT coefficients for the constant speed tests previously shown in Figures 7 and 8. Wavelet maps indicate the energy distribution with frequency on the vertical axis and time on the horizontal axis. A contour color scale was chosen so white-to-blueish colors represent regions of low energy and reddish-to-black regions represent higher energy. The magnitude of the coefficients was adjusted to use the same contour scales for both processed data and model results. The fault’s influence on the vibration in both time and frequency is easily recognizable, especially in the DGM simulation. The disruption due to the missing tooth appears as a pulse in the HWT of the damaged gear at approximately the rate of input shaft rotation, 820 rev/min or 13.7 Hz, for a period of 0.073 seconds. This figure demonstrates the capability of HWT and DGM to characterize a fault in the signal, but it also highlights a major benefit of using a DGM to build fault libraries, as the fault signature is more easily distinguishable in the signal than it is in the damaged experimental data.


## 4.3. Comparison of Experimental Data and Model Results under Nonstationary Conditions

The time-frequency analysis capability is especially important for condition monitoring when machinery speed is not constant. One important example is the operation of a wind turbine whose rotational speed is, in part, dependent on wind speed. Data-driven condition monitoring approaches generally require periods of constant speed to build libraries and assess machinery, but DGM can mimic the nonstationary speed-time profile of the equipment and can be solved repeatedly with healthy gear and numerous simulated faults.To demonstrate the nonstationary capability, the 4-second motor speed-time profile used by [29] and shown in Figure 2 is used. The speed increases from stop to 10 Hz after one second and then increases up to 25 Hz at two seconds and maintains this speed to three seconds before decelerating to a stop. It is important to note that the speed-time profile described here represents the control message to the motor and that actual motor speed varies with each run. This variation can be mimicked by the DGM, but data-driven approaches are generally unable to account for this variation. Figure 10 presents a comparison of HWT processed test data with healthy gear and missing tooth and HWT processed model results with the respective measured speed-time profiles. Note that the speed-profile of the healthy test was used for the healthy DGM, and likewise the damaged gear test speed-profile was used for the damaged DGM. The measured motor speed multiplied by the number of teeth in the first-stage pinion is the anticipated gear mesh frequency (GMF1), and the second-stage pinion times the intermediate shaft speed is GMF2. These GMF frequency-time profiles and their integer multiples are the kinematic behavior of the system and are superimposed on the contour plots of Figure 10 for both gear stages. The energy distribution in the HWT plots matches many of these GMFs and multiples.Figure 10
HWT of experimental data and model results for varying speed test: (a) healthy data, (b) healthy model, (c) damaged (missing tooth) data, and (d) damaged model, with gear mesh frequencies (GMF) and integer harmonics superimposed.
(a)
(b)
(c)
(d)Visual comparisons of the HWT contour plots in Figure10 reveal obvious differences between the healthy and damaged cases in both experimental data and model results. While the healthy-damaged differences are not entirely consistent between data and model, there are fault feature trends in common, such as the increase in magnitude in GMF2 and widening of the regions near and around GMF1. While these contour plots aid in visualizing the wave map, the subtle numerical differences within the coefficient matrix contain feature patterns to be used to build libraries within pattern recognition algorithms.To demonstrate the effectiveness of this approach, Figure11 presents a contour plot of the subtraction of a damaged gear DGM HWT coefficient matrix from a healthy one when both use the same speed-profile. The ability to control the speed of the DGM precisely enables isolation of the effects of the damaged gear on the energy distribution in frequency and time. A noteworthy aspect is the similarity in the relatively constant speed region to the pulses found in the constant speed region of Figure 9. These pulses represent a potential feature that a pattern recognition algorithm can use when developing fault libraries.Figure 11
Contour plot of difference between HWT of healthy and damaged model results, with gear mesh frequencies (GMF) and integer harmonics superimposed.


## 5. Concluding Remarks

In this research, a comprehensive dynamic gearbox model (DGM) is developed for a benchmark two-stage gearbox test bed. The central premise of the proposed approach in this research is that gear faults can be included in DGMs and then processed with HWT to develop and improve condition monitoring pattern recognition libraries for fault classification and severity identification, ultimately leading to diagnosis and even prognosis. The primary advantage of using HWT and DGM to simulate gear faults is the ability to mimic the nonstationary operation of the physical gear unit, specifically torque load and changing speeds, in order to characterize the influence of the fault on the vibration signal. The obstacle this overcomes is seen even when attempting to keep test conditions identical; there is variation in the speed of the gear test unit, making direct comparisons between healthy and faulty gears difficult. Using the DGM as a go-between and the fault detection sensitivity of HWT, a technique is developed which promises to eliminate a portion of the uncertainty inherent in comparing nonstationary signals.Follow-on research will include investigations into other fault types and fault severity to develop the libraries of pattern recognition algorithms, thereby further merging the data-driven and model-driven approaches.



---

*Source: 9638325-2016-06-23.xml*