Gear fault diagnosis using gear meshing stiffness identified by gearbox housing vibration signals

Abstract

Gearbox fault diagnosis based on vibration sensing has drawn much attention for a long time. For highly integrated complicated mechanical systems, the intercoupling of structure transfer paths results in a great reduction or even change of signal characteristics during the process of original vibration transmission. Therefore, using gearbox housing vibration signal to identify gear meshing excitation signal is of great significance to eliminate the influence of structure transfer paths, but accompanied by huge scientific challenges. This paper establishes an analytical mathematical description of the whole transfer process from gear meshing excitation to housing vibration. The gear meshing stiffness (GMS) identification approach is proposed by using housing vibration signals for two stages of inversion based on the mathematical description. Specifically, the linear system equations of transfer path analysis are first inverted to identify the bearing dynamic forces. Then the dynamic differential equations are inverted to identify the GMS. Numerical simulation and experimental results demonstrate the proposed method can realize gear fault diagnosis better than the original housing vibration signal and has the potential to be generalized to other speeds and loads. Some interesting properties are discovered in the identified GMS spectra, and the results also validate the rationality of using meshing stiffness to describe the actual gear meshing process. The identified GMS has a clear physical meaning and is thus very useful for fault diagnosis of the complicated equipment.

Key words： gearbox fault diagnosis; meshing stiffness; identification; transfer path; signal processing

1 Introduction

As key machine components of the transmission system, gearboxes play an essential role in the aerospace, automotive, wind turbine, and power generation industries [1]. To avoid unwanted downtime, expensive repair procedures, and even human casualties, gearbox condition monitoring and fault diagnosis based on vibration sensing have drawn much attention for a long time.

One type of mainstream methods extracts fault features [2] from the time domain or frequency domain of the vibration signal, and then determines whether the machine faulty according to feature changing and fault mechanism [3]. Furthermore, the fault features can be used to describe the machine degradation process [4] with the help of some intuitive properties (e.g., monotonicity [5] and divisibility [6]). To ensure the extracted features are more sensitive to the fault, signal decomposition methods [7-9] decompose the fault characteristic components from the original vibration signal. Signal decomposition effectively avoids the interference of signal components unrelated to the fault and thus succeeds in the field of fault diagnosis. In addition, some signal decomposition methods (e.g., nonlinear chirp mode decomposition [10], variational mode decomposition [11], and synchrosqueezing transform [12]) are able to deal with non-stationary conditions, and have been widely used in fault diagnosis of flywheel bearing [13], wind turbine gearbox [14], and aviation piston engine [15].

The other type of mainstream methods is to design different neural networks to let the computer automatically learn the implicit fault mode in the vibration signal, such as deep belief network [16], sparse filtering [17], deep transfer learning [18], and stacked denoising autoencoders [19]. This type of methods weakens fault mechanism because the artfully designed neural networks can obtain complex fault features [20] beyond human understanding by learning [21,22]. However, for the life cycle [23] of an equipment, that fault data can be acquired only at the last moment when the fault occurs, causing it very difficult to obtain fault data in the actual industry [24]. Therefore, small samples [25] and generalization [26] are the focus of attention.

The above two types of mainstream methods depend on the acquisition of vibration signals with good quality. In other words, they assume the vibration signals in the healthy state and fault state can be distinguished. In terms of mechanism, gearbox vibration originates mainly from the meshing action of the gear teeth, then is transmitted through the shafts and bearings to the housing [27]. The structure transfer paths have less influence for simple rotary machines working under high-speed and heavy-duty conditions, and vibration transducers can capture the fault characteristic. The above assumption is established. However, for some highly integrated complicated mechanical systems, e.g., armored vehicles, installing the transducers in a location like bearing support is impossible. Unlike the laboratory vibration test, the intercoupling of structure transfer paths results in a great reduction or even change of signal characteristics during the process of original vibration transmission. Therefore, using gearbox housing vibration signal to inversely identify gear meshing excitation signal for eliminating the influence of structure transfer paths is of great significance, but accompanied by huge scientific challenges.

For the inverse problem of excitation identification in mechanical systems, Moore−Penrose pseudo-inverse [28], truncated generalized singular value decomposition [29], and regularization method [30] are popularly used. Research results of bridge structure’s dynamic load identification [31], architecture’s wind load identification [32], propeller blade’s dynamic load identification [33], and helicopter rotor’s load identification [34] have been reported. However, it is almost undeveloped for gear meshing excitation identification due to the lack of related theories about the forward modeling and the reverse identification of gear meshing vibration transfer. The latest research result realizes meshing force measurement and identification of gear ring using the dynamic model and vibration signal [35,36]. The modeling ignores the influence of structure transfer paths and depends on the accurate calibration of the test process. Thus, one cannot use the method to deal with housing vibration signals in real mechanical systems. Our previous work used the bearing dynamic force as a bridge connecting gear meshing excitation and housing vibration and proved that the fault features of bearing dynamic force are more sensitive to fault than the original vibration signal [37]. Under the proposed modeling framework, the structure transfer path is modeled by frequency response functions (FRFs) in the frequency domain [38]. It can be used to identify bearing dynamic forces by housing vibration signals. Nevertheless, gear meshing is parametric excitation in the view of dynamics. Therefore, the transfer path equation cannot describe the vibration transfer process from gear meshing excitation to bearing dynamic force.

In this paper, an analytical mathematical description of the entire transfer process from gear meshing excitation to housing vibration is established by using joint modeling of dynamics and transfer path analysis. We propose a two-stage inversion approach to identify gear meshing stiffness (GMS) using housing vibration signals based on the analytical model. Specifically, the linear system equations of transfer path analysis are first inverted to identify the bearing dynamic forces. Then the dynamic differential equations are inverted to identify the GMS. Numerical simulation and experimental results show that the identified GMS has some interesting properties. For example, signal components unrelated to faults are suppressed, and the spectrum noise level is related to fault severity. The results demonstrate the proposed approach can realize gear fault diagnosis and has the potential to generalize to other speeds and loads. Since the frequency components of the identified GMS are very similar to theoretical modeling, the research results also validate the rationality of using time-varying stiffness to describe the actual gear meshing process. The identified GMS has a clear physical meaning and thus is beneficial for fault diagnosis and classification of complicated equipment.

The following content of this paper is divided into five parts. In Section 2, joint modeling of gearbox vibration based on dynamics and transfer path is proposed to describe the vibration transfer process from gear meshing excitation to housing vibration. Section 3 elaborates the proposed gear meshing stiffness identification (GMSI) approach with detailed formula derivation. In Section 4, the effectiveness of the proposed method is demonstrated through a simulation example, with classifying healthy, pitting, wear, and fracture gears. The experimental verification on gearbox fault classification is performed in Section 5. The influence of measuring points’ location and the generalization ability of the approach are also discussed. Finally, related conclusions are summarized in Section 6.

2 Joint modeling of gearbox vibration based on dynamics and transfer path

The transmission process from gear meshing to gearbox housing vibration needs to be investigated to realize GMSI based on gearbox housing vibration signal. We modeled the gearbox vibration by combining dynamics and transfer path analysis, as shown in Fig.1. Specifically, the vibration transmission process from gear meshing to bearing dynamic force is described by dynamic modeling of the gear−shaft−bearing system. Then the bearing dynamic force is seen as excitation on the housing interface, and transfer path analysis is used to describe the vibration transmission process from bearing dynamic force to gearbox housing vibration.

Fig.1 Gearbox vibration modeling by combining dynamics and transfer path analysis.

Full size|PPT slide

2.1 Dynamic modeling of gear−shaft−bearing system

The shaft in a gear system is usually treated as a beam component in the dynamic modeling. Based on the idea of element division in the finite element method (FEM). The shaft is divided into many sections, with each section simulated by a two-node Timoshenko beam. Considering lateral deformation and torsional deformation, the degree of freedom (DOF) of the beam model can be represented as below:

(1){\mathit{u}}^{\mathrm{e}}={\left[\begin{array}{cccccccccccc}{x}_A& y_A& {\mathit{\text{z}}}_A& {\theta }_{xA}& {\theta }_{yA}& {\theta }_{\mathit{\text{z}}A}& x_B& y_B& {\mathit{\text{z}}}_B& {\theta }_{xB}& {\theta }_{yB}& {\theta }_{\mathit{\text{z}}B}\end{array}\right]}^{\mathrm{T}},

where {\mathit{u}}^{\mathrm{e}} denotes the displacement vectors of nodes A and B on the Timoshenko beam, x, y, and z denote the translational DOFs, and {\theta }_x, {\theta }_y, and {\theta }_{\mathit{\text{z}}} denote the rotational DOFs.

Then the mass matrix, stiffness matrix, and gyroscopic torque matrix of each DOF can be calculated respectively based on the classical theory of beam element [39]. Note that the speed in the gyroscopic torque matrix is the speed of the corresponding shaft.

The bearing in a gear system is the supporting part of the shaft. An uncoupled elastic support model is used to simulate each bearing direction as an elastic element with equivalent stiffness, without considering the coupling effect between each DOF. And the motion differential equation of the bearing used in this paper is as below:

(2){\mathit{k}}_{\mathrm{B}}{\mathit{x}}_{\mathrm{B}}=\text{diag}\left(\begin{array}{cccccc}{k}_{xx}& k_{yy}& k_{\mathit{\text{zz}}}& k_{{\theta }_x{\theta }_x}& k_{{\theta }_y{\theta }_y}& 0\end{array}\right)\left[\begin{array}{c}{x}_{\mathrm{b}}\\ y_{\mathrm{b}}\\ {\mathit{\text{z}}}_{\mathrm{b}}\\ {\theta }_{x\mathrm{b}}\\ {\theta }_{y\mathrm{b}}\\ {\theta }_{\mathit{\text{z}}\mathrm{b}}\end{array}\right]={\mathit{F}}_{\mathrm{B}},

where {\mathit{F}}_{\mathrm{B}} denotes the bearing force, {\mathit{k}}_{\mathrm{B}} and {\mathit{x}}_{\mathrm{B}} denote the stiffness vector and the displacement vector, k_{xx}, k_{yy}, k_{\mathit{\text{zz}}}, k_{{\theta }_x{\theta }_x}, and k_{{\theta }_y{\theta }_y} denote different translational and rotational DOFs of bearing stiffness, and x_{\mathrm{b}}, y_{\mathrm{b}}, {\mathit{\text{z}}}_{\mathrm{b}}, {\theta }_{x\mathrm{b}}, {\theta }_{y\mathrm{b}}, and {\theta }_{\mathit{\text{z}}\mathrm{b}} denote the DOFs of bearing.

The gear in a gear system is modeled as a lumped mass model with different mass and moment of inertia in different directions. Six DOFs in translation and rotation are considered, same as a node. And the mass matrix {\mathit{M}}_{\mathrm{g}} is as below:

(3){\mathit{M}}_{\mathrm{g}}=\text{diag}\left(\begin{array}{cccccc}{m}_{\mathrm{g}}& m_{\mathrm{g}}& m_{\mathrm{g}}& I_x& I_y& I_z\end{array}\right),

where {\mathit{M}}_{\mathrm{g}} denotes the mass matrix of a lumped mass model, m_{\mathrm{g}} is gear mass, and I_x, I_y, and I_z denote the moment of inertia in different directions.

The process of gear meshing is modeled as a spring. And the stiffness of the spring is the GMS. Suppose the gears are spur gears with zero helix angle, the generalized coordinate of the gear pair can be defined as

(4){\mathit{X}}_{12}={\left[\begin{array}{cccccccccccc}{x}_1& y_1& {\mathit{\text{z}}}_1& {\theta }_{x1}& {\theta }_{y1}& {\theta }_{\mathit{\text{z}}1}& x_2& y_2& {\mathit{\text{z}}}_2& {\theta }_{x2}& {\theta }_{y2}& {\theta }_{\mathit{\text{z}}2}\end{array}\right]}^{\mathrm{T}},

where {\mathit{X}}_{12} denotes the generalized coordinate of the gear pair, and x_i, y_i, {\mathit{\text{z}}}_i, {\theta }_{xi}, {\theta }_{yi}, and {\theta }_{\mathit{\text{z}}i} denote the DOFs of gear i (i = 1,2).

The motion differential equation of the gear pair can be established by ignoring the influence of tooth backlash, friction, and meshing damping:

(5){\mathit{M}}_{12}{\ddot{\mathit{X}}}_{12}+{\mathit{K}}_{12}{\mathit{X}}_{12}={\mathit{F}}_{12}+{\mathit{F}}_{\mathrm{e}},

where {\mathit{M}}_{12} is the mass matrix of the gear pair, {\mathit{K}}_{12} is the meshing stiffness matrix of the gear pair, {\mathit{F}}_{12} is the external load, {\mathit{F}}_{\mathrm{e}} is the excitation caused by no-load transmission error of the gear pair, and {\ddot{\mathit{X}}}_{12} denotes the second derivative of {\mathit{X}}_{12}. The formulas of them can be written as below:

(6){\mathit{M}}_{12}=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\left(\begin{array}{cccccccccccc}{m}_1& m_1& m_1& I_1& I_1& J_1& m_2& m_2& m_2& I_2& I_2& J_2\end{array}\right),

(7){\mathit{K}}_{12}=k_{12}{\mathit{\alpha }}_{12}^{\mathrm{T}}{\mathit{\alpha }}_{12},

(8)\begin{array}{rl}{\mathit{\alpha }}_{12}=& [\begin{array}{cccccc}-\sin {\psi }_{12}& \cos {\psi }_{12}& 0& 0& 0& \mathrm{s}\mathrm{g}\mathrm{n}\left(r_{\mathrm{b}1}\right)\end{array}\\ & \begin{array}{cccccc}\sin {\psi }_{12}& -\cos {\psi }_{12}& 0& 0& 0& \mathrm{s}\mathrm{g}\mathrm{n}\left(r_{\mathrm{b}2}\right)\end{array}],\end{array}

(9){\mathit{F}}_{12}={\left[\begin{array}{cccccccccccc}0& 0& 0& 0& 0& \mathrm{s}\mathrm{g}\mathrm{n}\left(T_1\right)& 0& 0& 0& 0& 0& \mathrm{s}\mathrm{g}\mathrm{n}\left(T_2\right)\end{array}\right]}^{\mathrm{T}},

where m_1 and m_2 are the masses of gears 1 and 2, respectively, I_1 and I_2 are the transverse moments of inertia of gears 1 and 2, respectively, J_1 and J_2 is the polar moments of inertia of gears 1 and 2, respectively, k_{12} is the meshing stiffness of the gear pair, {\mathit{\alpha }}_{12} is the projection vector of the gear pair, sgn(∙) denotes the influence of gear rotation direction with the value 1 for driving gear rotates anticlockwise and −1 for driving gear rotates clockwise, r_{\mathrm{b}1} and r_{\mathrm{b}2} denote radii of base circle, T_1 and T_2 denote the torques applied to the input shaft and output shaft, respectively, and {\psi }_{12} is the angle between the positive y-axis and the meshing surface.

Assemble the mass matrix, stiffness matrix, and damping matrix of the gear‒shaft‒bearing system based on the finite element theory, as shown in Fig.2. n_i denotes the node number and {\mathit{K}}_i^{\text{s}} denotes the unit stiffness matrix of the system. The motion differential equation of the whole gear‒shaft‒bearing system can be represented as below:

Fig.2 Assembly method of the gear−shaft−bearing system matrix.

Full size|PPT slide

(10)\mathit{M}\ddot{\mathit{u}}+\mathit{D}\dot{\mathit{u}}+\mathit{K}\mathit{u}=\mathit{F},

where M, D, and K denote the system mass matrix, the system damping matrix, and the system stiffness matrix, respectively, u denotes generalized coordinate of the system, also known as nodes displacement, \dot{\mathit{u}} and \ddot{\mathit{u}} denote the first and second derivative of \mathit{u}, respectively, and F is the external excitation. The mass matrix M contains shaft mass, gear mass, and lumped mass. The damping matrix D contains viscous damping, bearing damping, and gyroscopic torque. The stiffness matrix K contains shaft stiffness, meshing stiffness, and bearing stiffness. The viscous damping in the damping matrix D is described by Rayleigh damping [40].

2.2 Transfer path modeling of gearbox vibration

Without loss of generality, the gearbox to be analyzed is seen as a linear time-invariant system. Consider ignoring the axial vibration of the rotor, and the bearing dynamic force can be applied on the housing interface as an excitation vector to achieve coupling transfer path modeling [37,41] in the frequency domain:

(11)\mathit{x}=\mathit{H}{\mathit{F}}_{\mathrm{b}},

where {\mathit{F}}_{\mathrm{b}}={\left[f_1\left(\omega \right),f_2\left(\omega \right),...,f_m\left(\omega \right)\right]}^{\mathrm{T}} represents the column vector composed of m bearing dynamic force excitation spectra, \mathit{x}={\left[x_1\left(\omega \right),x_2\left(\omega \right),...,x_n\left(\omega \right)\right]}^{\mathrm{T}} represents the column vector composed of n response channels spectra, \omega represents angular frequency, and H represents the transfer function matrix composed of n\times m FRFs.

To ensure that structure transfer paths are decoupled in different directions and eliminate the mutual coupling effect among structure transfer paths on real and imaginary parts of the bearing dynamic force excitation signal, Eq. (11) can be modified as below:

(12){\mathit{x}}_{\mathrm{d}}={\mathit{H}}_{\mathrm{d}}{\mathit{F}}_{\mathrm{d}},

where {\mathit{x}}_{\mathrm{d}}, {\mathit{H}}_{\mathrm{d}}, and {\mathit{F}}_{\mathrm{d}} denote the decoupled response channels spectra vector, the decoupled transfer function matrix, and the decoupled bearing dynamic force excitation spectra vector, respectively.

(13){\mathit{H}}_{\mathrm{d}}=\left[\begin{array}{cccc}{\mathit{H}}_{11}\left(\omega \right)& {\mathit{H}}_{12}\left(\omega \right)& \cdots & {\mathit{H}}_{1m}\left(\omega \right)\\ {\mathit{H}}_{21}\left(\omega \right)& {\mathit{H}}_{22}\left(\omega \right)& \cdots & {\mathit{H}}_{2m}\left(\omega \right)\\ ⋮& ⋮& ⋮& ⋮\\ {\mathit{H}}_{n1}\left(\omega \right)& {\mathit{H}}_{n2}\left(\omega \right)& \cdots & {\mathit{H}}_{nm}\left(\omega \right)\end{array}\right],

(14){\mathit{H}}_{ij}\left(\omega \right)=\left[\begin{array}{cccccc}{H}_{ix\mathrm{R}jx\mathrm{R}}\left(\omega \right)& H_{ix\mathrm{I}jx\mathrm{R}}\left(\omega \right)& H_{iy\mathrm{R}jx\mathrm{R}}\left(\omega \right)& H_{iy\mathrm{I}jx\mathrm{R}}\left(\omega \right)& H_{iz\mathrm{R}jx\mathrm{R}}\left(\omega \right)& H_{iz\mathrm{I}jx\mathrm{R}}\left(\omega \right)\\ H_{ix\mathrm{R}jx\mathrm{I}}\left(\omega \right)& H_{ix\mathrm{I}jx\mathrm{I}}\left(\omega \right)& H_{iy\mathrm{R}jx\mathrm{I}}\left(\omega \right)& H_{iy\mathrm{I}jx\mathrm{I}}\left(\omega \right)& H_{iz\mathrm{R}jx\mathrm{I}}\left(\omega \right)& H_{iz\mathrm{I}jx\mathrm{I}}\left(\omega \right)\\ H_{ix\mathrm{R}jy\mathrm{R}}\left(\omega \right)& H_{ix\mathrm{I}jy\mathrm{R}}\left(\omega \right)& H_{iy\mathrm{R}jy\mathrm{R}}\left(\omega \right)& H_{iy\mathrm{I}jy\mathrm{R}}\left(\omega \right)& H_{iz\mathrm{R}jy\mathrm{R}}\left(\omega \right)& H_{iz\mathrm{I}jy\mathrm{R}}\left(\omega \right)\\ H_{ix\mathrm{R}jy\mathrm{I}}\left(\omega \right)& H_{ix\mathrm{I}jx\mathrm{I}}\left(\omega \right)& H_{iy\mathrm{R}jy\mathrm{I}}\left(\omega \right)& H_{iy\mathrm{I}jx\mathrm{I}}\left(\omega \right)& H_{iz\mathrm{R}jy\mathrm{I}}\left(\omega \right)& H_{iz\mathrm{I}jx\mathrm{I}}\left(\omega \right)\\ H_{ix\mathrm{R}jz\mathrm{R}}\left(\omega \right)& H_{ix\mathrm{I}jz\mathrm{R}}\left(\omega \right)& H_{iy\mathrm{R}jz\mathrm{R}}\left(\omega \right)& H_{iy\mathrm{I}jz\mathrm{R}}\left(\omega \right)& H_{iz\mathrm{R}jz\mathrm{R}}\left(\omega \right)& H_{iz\mathrm{I}jz\mathrm{R}}\left(\omega \right)\\ H_{ix\mathrm{R}jz\mathrm{I}}\left(\omega \right)& H_{ix\mathrm{I}jz\mathrm{I}}\left(\omega \right)& H_{iy\mathrm{R}jz\mathrm{I}}\left(\omega \right)& H_{iy\mathrm{I}jz\mathrm{I}}\left(\omega \right)& H_{iz\mathrm{R}jz\mathrm{I}}\left(\omega \right)& H_{iz\mathrm{I}jz\mathrm{I}}\left(\omega \right)\end{array}\right],

(15)\begin{array}{rl}& {\mathit{F}}_{\mathrm{d}}={\left[{\mathit{F}}_1,{\mathit{F}}_2,...,{\mathit{F}}_m\right]}^{\mathrm{T}},\\ & {\mathit{F}}_i={\left[f_{ix\mathrm{R}}\left(\omega \right),f_{ix\mathrm{I}}\left(\omega \right),f_{iy\mathrm{R}}\left(\omega \right),f_{iy\mathrm{I}}\left(\omega \right),f_{iz\mathrm{R}}\left(\omega \right),f_{iz\mathrm{I}}\left(\omega \right)\right]}^{\mathrm{T}},\end{array}

(16)\begin{array}{rl}& {\mathit{x}}_{\mathrm{d}}={\left[{\mathit{x}}_1,{\mathit{x}}_2,...,{\mathit{x}}_n\right]}^{\mathrm{T}},\\ & {\mathit{x}}_j={\left[x_{jx\mathrm{R}}\left(\omega \right),x_{jx\mathrm{I}}\left(\omega \right),x_{jy\mathrm{R}}\left(\omega \right),x_{jy\mathrm{I}}\left(\omega \right),x_{jz\mathrm{R}}\left(\omega \right),x_{jz\mathrm{I}}\left(\omega \right)\right]}^{\mathrm{T}},\end{array}

where H_{ix\mathrm{R}jx\mathrm{R}}\left(\omega \right) represents the FRF from the real part of x component at the excitation channel i to the real part of x component at the response channel j, and can be calculated by

(17)H_{ix\mathrm{R}jx\mathrm{R}}\left(\omega \right)=x_{jx\mathrm{R}}\left(\omega \right)/\phantom{x_{jx\mathrm{R}}\left(\omega \right)f_{ix\mathrm{R}}\left(\omega \right)}{f}_{ix\mathrm{R}}\left(\omega \right),

where x_{jx\mathrm{R}}\left(\omega \right) represents the real part of the x component of housing vibration response at location j, f_{ix\mathrm{R}}\left(\omega \right) represents the real part of the x component of bearing dynamic force at bearing i, and the other elements are similar. The FRFs can be measured by hammer experiment or numerically simulated by FEM.

3 GMSI approach

In this section, we propose a GMSI method using gearbox housing vibration signal.

3.1 Method

Under the modeling framework described in Section 2, the bearing dynamic forces need to be identified by the housing vibration first. In most cases the vibration response channels are less than bearing dynamic force excitation channels, thus the transfer function matrix is not a square matrix, and Eq. (12) needs to be modified as

(18){\mathit{H}}_{\mathrm{d}}^{\mathrm{T}}{\mathit{x}}_{\mathrm{d}}={\mathit{H}}_{\mathrm{d}}^{\mathrm{T}}{\mathit{H}}_{\mathrm{d}}{\mathit{F}}_{\mathrm{d}}.

Using least square method to solve Eq. (18) will cause a large error due to the ill-condition, especially for underdetermined conditions. Therefore, a weighted iteration approach is proposed to control the error. Performing weighted decomposition to {\mathit{F}}_{\mathrm{d}}, namely {\mathit{F}}_{\mathrm{d}}=\mathit{W}{\mathit{f}}_{\mathrm{d}}, and substitute it into Eq. (18), then the basic idea of the weighted iteration approach is derived as follows:

(19){\mathit{H}}_{\mathrm{d}}^{\mathrm{T}}{\mathit{x}}_{\mathrm{d}}={\mathit{H}}_{\mathrm{d}}^{\mathrm{T}}{\mathit{H}}_{\mathrm{d}}\mathit{W}{\mathit{f}}_{\mathrm{d}},

where W is the weighted matrix used to regulates the ill-condition degree of {\mathit{H}}_{\mathrm{d}}^{\mathrm{T}}{\mathit{H}}_{\mathrm{d}}, the upper bound of identification error can then be derived as below, similar with Ref. [42]:

(20)\begin{array}{rl}\frac{\Vert \delta {\mathit{f}}_{\mathrm{d}}\Vert }{\Vert {\mathit{f}}_{\mathrm{d}}\Vert }⩽& k\left(\widetilde{{\mathit{H}}_{\mathrm{d}}^{\mathrm{T}}{\mathit{H}}_{\mathrm{d}}}\mathit{W}\right)\frac{\Vert \mathit{W}\Vert \Vert {\mathit{H}}_{\mathrm{d}}^{\mathrm{T}}{\mathit{H}}_{\mathrm{d}}\Vert \Vert \delta {\mathit{H}}_{\mathrm{d}}^{\mathrm{T}}{\mathit{x}}_{\mathrm{d}}\Vert }{\Vert \widetilde{{\mathit{H}}_{\mathrm{d}}^{\mathrm{T}}{\mathit{H}}_{\mathrm{d}}}\mathit{W}\Vert \Vert {\mathit{H}}_{\mathrm{d}}^{\mathrm{T}}{\mathit{x}}_{\mathrm{d}}\Vert }\\ & +k\left(\widetilde{{\mathit{H}}_{\mathrm{d}}^{\mathrm{T}}{\mathit{H}}_{\mathrm{d}}}\mathit{W}\right)\frac{\Vert \delta {\mathit{H}}_{\mathrm{d}}^{\mathrm{T}}{\mathit{H}}_{\mathrm{d}}\Vert }{\Vert \widetilde{{\mathit{H}}_{\mathrm{d}}^{\mathrm{T}}{\mathit{H}}_{\mathrm{d}}}\Vert },\end{array}

where k\left(\widetilde{{\mathit{H}}_{\mathrm{d}}^{\mathrm{T}}{\mathit{H}}_{\mathrm{d}}}\mathit{W}\right)=\Vert \widetilde{{\mathit{H}}_{\mathrm{d}}^{\mathrm{T}}{\mathit{H}}_{\mathrm{d}}}\mathit{W}\Vert \Vert {\left(\widetilde{{\mathit{H}}_{\mathrm{d}}^{\mathrm{T}}{\mathit{H}}_{\mathrm{d}}}\mathit{W}\right)}^{†}\Vert is the generalized condition number of \widetilde{{\mathit{H}}_{\mathrm{d}}^{\mathrm{T}}{\mathit{H}}_{\mathrm{d}}}\mathit{W}, the terms without superscript indicate analytical values, the superscript ~ indicates the term containing a numerical error, and the numerical error is represented by the prefix \delta. The upper bound of identification error without weighted decomposition can be derived similarly:

(21)\frac{\Vert \delta {\mathit{F}}_{\mathrm{d}}\Vert }{\Vert {\mathit{F}}_{\mathrm{d}}\Vert }⩽k\left(\widetilde{{\mathit{H}}_{\mathrm{d}}^{\mathrm{T}}{\mathit{H}}_{\mathrm{d}}}\right)\frac{\Vert {\mathit{H}}_{\mathrm{d}}^{\mathrm{T}}{\mathit{H}}_{\mathrm{d}}\Vert \Vert \delta {\mathit{H}}_{\mathrm{d}}^{\mathrm{T}}{\mathit{x}}_{\mathrm{d}}\Vert }{\Vert \widetilde{{\mathit{H}}_{\mathrm{d}}^{\mathrm{T}}{\mathit{H}}_{\mathrm{d}}}\Vert \Vert {\mathit{H}}_{\mathrm{d}}^{\mathrm{T}}{\mathit{x}}_{\mathrm{d}}\Vert }+k\left(\widetilde{{\mathit{H}}_{\mathrm{d}}^{\mathrm{T}}{\mathit{H}}_{\mathrm{d}}}\right)\frac{\Vert \delta {\mathit{H}}_{\mathrm{d}}^{\mathrm{T}}{\mathit{H}}_{\mathrm{d}}\Vert }{\Vert \widetilde{{\mathit{H}}_{\mathrm{d}}^{\mathrm{T}}{\mathit{H}}_{\mathrm{d}}}\Vert }.

Comparing Eq. (20) with Eq. (21), to ensure the identification error is reduced after performing weighted decomposition, the following constraints should be satisfied:

(22)

Using A to represent {\mathit{H}}_{\mathrm{d}}^{\mathrm{T}}{\mathit{H}}_{\mathrm{d}} and a_{ij} to represent the elements, a diagonal weighted matrix is suggested as below:

(23)\begin{array}{rl}& \mathit{W}=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\left(\begin{array}{cccc}{w}_1& w_2& \cdots & w_n\end{array}\right),w_i=\sqrt{{\left(\sum _{j=1}^n{\left|a_{ij}\right|}^2\right)}^{-1}},\\ & i=1,2,\dots ,n.\end{array}

Then the main diagonal elements of matrix {\mathit{H}}_{\mathrm{d}}^{\mathrm{T}}{\mathit{H}}_{\mathrm{d}}\mathit{W} are weighted by adding \alpha {\mathit{f}}_{\mathrm{d}} on both sides of Eq. (19), and the equation of the weighted iteration approach is derived as below:

(24)\left({\mathit{H}}_{\mathrm{d}}^{\mathrm{T}}{\mathit{H}}_{\mathrm{d}}\mathit{W}+\alpha \mathit{E}\right){\mathit{f}}_{\mathrm{d}}={\mathit{H}}_{\mathrm{d}}^{\mathrm{T}}{\mathit{x}}_{\mathrm{d}}+\alpha {\mathit{f}}_{\mathrm{d}},

where the weighted factor \alpha \left(\alpha >0\right) is similar to the regularization parameter in Tikhonov regularization [43], and can be chosen by generalized cross validation [44]. Equation (24) can be solved in an iterative form as below:

(25)\left({\mathit{H}}_{\mathrm{d}}^{\mathrm{T}}{\mathit{H}}_{\mathrm{d}}\mathit{W}+\alpha \mathit{E}\right){\mathit{f}}_{\mathrm{d}}^{\left(k+1\right)}={\mathit{H}}_{\mathrm{d}}^{\mathrm{T}}{\mathit{x}}_{\mathrm{d}}+\alpha {\mathit{f}}_{\mathrm{d}}^{\left(k\right)}.

Suppose that after K iterations {\mathit{f}}_{\mathrm{d}} converges to {\mathit{f}}_{\mathrm{d}}^{\left(K\right)}, the bearing dynamic force can be reconstructed in the frequency domain by

(26)\widetilde{{\mathit{F}}_{\mathrm{d}}}={\mathit{W}\mathit{f}}_{\mathrm{d}}^{\left(K\right)}.

The complex spectra of the bearing dynamic force \widetilde{{\mathit{F}}_{\mathrm{b}}} can be easily obtained using \widetilde{{\mathit{F}}_{\mathrm{d}}}, for example, f_i\left(\omega \right)= f_{ix\mathrm{R}}\left(\omega \right)+f_{ix\mathrm{I}}\left(\omega \right)\cdot \mathrm{i}, i denotes the imaginary unit. Then the vibration response of all bearing nodes {\tilde{\mathit{x}}}_{\mathrm{b}} can be calculated as

(27){\tilde{\mathit{x}}}_{\mathrm{b}}=⟨\widetilde{{\mathit{F}}_{\mathrm{b}}},{\mathit{k}}_{\mathrm{b}}⟩,

where \widetilde{{\mathit{F}}_{\mathrm{b}}}={\left[{\tilde{f}}_1\left(\omega \right),{\tilde{f}}_2\left(\omega \right),...,{\tilde{f}}_m\left(\omega \right)\right]}^{\mathrm{T}}, {\mathit{k}}_{\mathrm{b}}=[1/\phantom{1k_{\mathrm{b}1}}{k}_{\mathrm{b}1},1/\phantom{1k_{\mathrm{b}2}}{k}_{\mathrm{b}2},..., 1/\phantom{1k_{\mathrm{b}m}}{k}_{\mathrm{b}m}{]}^{\mathrm{T}}, and k_{\mathrm{b}i} denotes the equivalent stiffness of bearing corresponding to {\tilde{f}}_i\left(\omega \right).

Defining transition matrix {\mathit{H}}_{\mathrm{T}} in the frequency domain to convert the vibration response of bearing nodes to all nodes as below:

(28){\mathit{H}}_{\mathrm{T}}={\mathit{u}}_{\mathrm{f}}{\mathit{F}}_{\mathrm{b}}^{\mathrm{T}}{\left({\mathit{F}}_{\mathrm{b}}{\mathit{F}}_{\mathrm{b}}^{\mathrm{T}}\right)}^{-1},

where {\mathit{u}}_{\mathrm{f}}=\mathrm{F}\mathrm{T}\left[\mathit{u}\right], \mathrm{F}\mathrm{T}[\cdot ] denotes Fourier transform (FT). The vibration response of all nodes in the frequency domain can be estimated by

(29){\tilde{\mathit{u}}}_{\mathrm{f}}={\mathit{H}}_{\mathrm{T}}{\tilde{\mathit{x}}}_{\mathrm{b}}.

Then the inverse Fourier transform (IFT) is used to convert the estimated vibration response from the frequency domain to the time domain:

(30)\tilde{\mathit{u}}=\mathrm{I}\mathrm{F}\mathrm{T}\left[{\tilde{\mathit{u}}}_{\mathrm{f}}\right].

To numerical identify the system stiffness matrix, Eq. (10) is rewritten as a discrete form as below:

(31)\mathit{M}{\mathrm{\Delta }}^2\mathit{u}+\mathit{D}\mathrm{\Delta }\mathit{u}+\mathit{K}\mathit{u}=\mathit{F}.

Calculating the first- and second-order differences of the estimated vibration response \tilde{\mathit{u}}, and doing mathematical manipulation to Eq. (31), we can yield

(32)\tilde{\mathit{u}}{\tilde{\mathit{u}}}^{\mathrm{T}}{\tilde{\mathit{K}}}^{\mathrm{T}}=\tilde{\mathit{u}}{\left(\mathit{F}-\mathit{M}{\mathrm{\Delta }}^2\tilde{\mathit{u}}-\mathit{D}\mathrm{\Delta }\tilde{\mathit{u}}\right)}^{\mathrm{T}}.

Simply, let \mathit{V}=\tilde{\mathit{u}}{\tilde{\mathit{u}}}^{\mathrm{T}}, \mathit{Y}=\tilde{\mathit{u}}{\left(\mathit{F}-\mathit{M}{\mathrm{\Delta }}^2\tilde{\mathit{u}}-\mathit{D}\mathrm{\Delta }\tilde{\mathit{u}}\right)}^{\mathrm{T}}, and suppose the system stiffness matrix has p columns. We can convert Eq. (32) into positive definite linear equations as below:

(33)\mathit{V}{\tilde{\mathit{K}}}^{\mathrm{T}}\left(:,i\right)=\mathit{Y}\left(:,i\right),i=1,2,\dots ,p,

where \left(:,i\right) indicates the ith column of the matrix. The generalized minimal residual algorithm [45] is suggested to solve the linear equations to estimate the system stiffness matrix \tilde{\mathit{K}}.

Suppose the gear−shaft−bearing system has q gear pairs, the whole meshing stiffness matrix {\tilde{\mathit{K}}}_{\mathrm{m}} can be estimated by subtracting the stiffness matrix of gears, shafts, and bearings from the estimated system stiffness matrix as below:

(34){\tilde{\mathit{K}}}_{\mathrm{m}}=\tilde{\mathit{K}}-\left(\mathit{K}-\sum _{i=1}^q{\mathit{K}}_{\mathrm{m}i}\right),

where {\mathit{K}}_{\mathrm{m}i} denotes the meshing stiffness matrix of gear pair i, for single-stage gear−shaft−bearing system, {\mathit{K}}_{\mathrm{m}}={\mathit{K}}_{\mathrm{m}1}={\mathit{K}}_{12}.

The meshing stiffness matrix can be rewritten as the following form according to Eq. (7):

(35)\left(\mathit{K}-\sum _{i=1}^q{\mathit{K}}_{\mathrm{m}i}\right)=\left[\begin{array}{ccc}\begin{array}{ccc}{k}_{\mathrm{m}1}& \cdots & 0\\ ⋮& \ddots & ⋮\\ 0& \cdots & k_{\mathrm{m}1}\end{array}& \cdots & \begin{array}{ccc}{k}_{\mathrm{m}q}& \cdots & 0\\ ⋮& \ddots & ⋮\\ 0& \cdots & k_{\mathrm{m}q}\end{array}\end{array}\right]\mathit{\psi },

(36)\mathit{\psi }=\left[\begin{array}{c}{\mathbf{\Theta }}_1\\ {\mathbf{\Theta }}_2\\ ⋮\\ {\mathbf{\Theta }}_q\end{array}\right],

where k_{\mathrm{m}i} denotes the meshing stiffness of gear pair i, ψ denotes the projection matrix of all gear nodes, and {\mathbf{\Theta }}_i represents the projection matrix obtained by assembling the projection vector {\mathit{\alpha }}_i^{\mathrm{T}}{\mathit{\alpha }}_i (can be calculated based on Eq. (8) for each gear pair) to the corresponding gear nodes. Note that the dimension of {\mathbf{\Theta }}_i is the same as the system stiffness matrix K.

Then we can derive the temporary variable containing meshing stiffness of all gear pairs by the least square method:

(37){\tilde{\mathit{k}}}_{\mathrm{m}}=\left[\begin{array}{ccc}\begin{array}{ccc}{\tilde{k}}_{\mathrm{m}1}& \cdots & 0\\ ⋮& \ddots & ⋮\\ 0& \cdots & {\tilde{k}}_{\mathrm{m}1}\end{array}& \cdots & \begin{array}{ccc}{\tilde{k}}_{\mathrm{m}q}& \cdots & 0\\ ⋮& \ddots & ⋮\\ 0& \cdots & {\tilde{k}}_{\mathrm{m}q}\end{array}\end{array}\right]={\tilde{\mathit{K}}}_{\mathrm{m}}{\mathit{\psi }}^{\mathrm{T}}{\left(\mathit{\psi }{\mathit{\psi }}^{\mathrm{T}}\right)}^{-1},

where {\tilde{\mathit{k}}}_{\mathrm{m}} denotes the temporary variable containing meshing stiffness of all gear pairs, {\tilde{\mathit{K}}}_{\mathrm{m}} denotes the estimated whole meshing stiffness matrix, and {\tilde{k}}_{\mathrm{m}i} denotes the estimated meshing stiffness of gear pair i.

The meshing stiffness of each gear pair can be finally identified as below:

(38){\tilde{k}}_{\mathrm{m}i}=\frac{1}{p}\sum _{j=\left(i-1\right)p+1}^{ip}{\tilde{\mathit{k}}}_{\mathrm{m}}\left(:,j\right),i=1,2,\dots ,q.

The trend component in the {\tilde{k}}_{\mathrm{m}i} can be eliminated by averaging the results of each 360° rotation period.

Section 4 will demonstrate that the identified GMS is different in the time and frequency domains for different gear faults. Therefore, the identified GMS can be used for gear fault diagnosis and classification based on the similarity index.

3.2 Fault diagnosis process in engineering practice

The proposed GMSI with its application can be summarized in Fig.3 where θ indicates the rotation angle of the gear pair during meshing process. Clearly, the fault diagnosis process in engineering practice using the proposed method is summarized as follows:

Fig.3 Flowchart of gear meshing stiffness identification with its application.

Full size|PPT slide

(i) According to the theory described in Section 2.1, building the dynamic model of the gear−shaft−bearing system (see Eq. (10)) by the mechanical parameters. Then constructing the system mass matrix M, the system damping matrix D, the external load vector F, and calculating the transition matrix (see Eq. (28)) in the frequency domain.

(ii) Testing the vibration signal x of the gearbox and obtaining the transfer function matrix H of the gearbox housing by hammer experiment or FEM model, then establishing the transfer path equation (see Eq. (11)) according to the theory described in Section 2.2.

(iii) Solving the inverse problem of the transfer path equation to identify bearing dynamic forces and estimate displacement responses of the bearing nodes based on Eqs. (18)−(27).

(iv) Estimating all node displacement responses of the gear−shaft−bearing system by displacement responses of the bearing nodes and transition matrix. Then calculate the inverse problem of the dynamic equation to estimate the system stiffness matrix based on Eqs. (30)−(38). The GMS can be extracted from the estimated system stiffness matrix.

(v) Combining the identified GMS with the similarity index to diagnose or classify gear faults.

4 Numerical simulation and results

This section will demonstrate the effectiveness of the proposed GMSI method through a simulated single-stage gear−shaft−bearing system fault classification example.

4.1 Housing vibration construction

For a single-stage gear−shaft−bearing system composed of two spur gears, two shafts, and four bearings, suppose k_{{\theta }_x{\theta }_x} and k_{{\theta }_y{\theta }_y} in Eq. (2) are all zero, and k_{xx}=k_{yy}=1.1837\times {10}^8 N/m, k_{\mathit{\text{zz}}}=1\times {10}^8 N/m. The parameters of gears and shafts are listed in Tab.1.

Tab.1 Parameters of gears and shafts

Gear	Teeth number	Face width/mm	Inner diameter/mm	Module	Shaft length/mm
Gear at driving shaft	25	16	25	2.5	405
Gear at driven shaft	58	16	25	2.5	405

The GMS is calculated using the loaded tooth contact analysis method [46,47]. This method is competent in stiffness calculation under various fault types. We use 20 points to describe the meshing stiffness of the gear pair in a single-double tooth mesh cycle. Then the sampling frequency is 5000 Hz when the driving shaft speed is 600 r/min without loss of generality. The torque is 60 N·m, and clearly, the excitation caused by rotor eccentric is ignored. The GMS waveforms of healthy and two typical gear faults (pitting and wear) are given in Fig.4(a). The theoretical meshing stiffness spectra of healthy and the above two typical gear faults are also calculated and shown in Fig.4(b). Hereinafter the direct current component is subtracted from the GMS signal to show the spectra clearly. Compared with healthy gear, the magnitudes at meshing frequency (f_{\mathrm{m}}) and its harmonic components of wear gear reduce obviously. The magnitudes at meshing frequency and its harmonic components of pitting gear increase slightly, accompanying unobvious sideband frequencies appearing in the spectrum. In the following, we will demonstrate that some of the above properties still exist in the identified GMS.

Fig.4 Simulated gear meshing stiffness (a) waveforms and (b) spectra of healthy and two typical gear faults.

Full size|PPT slide

The eight bearing dynamic forces at four bearings (1x, 1z, 2x, 2z, 3x, 3z, 4x, 4z) are calculated based on Eq. (10) and the above-defined dynamic parameters. Then the bearing dynamic forces are used to simulate housing vibration responses by Eqs. (12)‒(17). Herein the vibration signals of two triaxial accelerometers are constructed to simulate more difficult underdetermined conditions in engineering. The spectra of the six housing vibration signals for healthy gear are given in Fig.5, which mainly contain the meshing frequency and its harmonic components. There are sideband frequencies modulated by the rotating frequency of the driving shaft at each meshing frequency. The simulation method and results are similar for fault gears. Affected by the structure transfer paths, the magnitudes of the spectra on different housing locations are very different, although the frequency components are the same.

Fig.5 Simulated housing vibration response at two measuring points.

Full size|PPT slide

4.2 Simulated GMSI

Using the proposed GMSI method to solve the inverse problem of Eq. (12), we can eliminate the influence of housing structure transfer paths and identify all the bearing dynamic forces first. Considering the gearbox structure in engineering is usually complicated, such inverse problem is underdetermined in most cases. Therefore, controlling the error is crucial for the subsequent GMSI. The identified bearing dynamic forces compared with theoretical ones for the simulation example are given in Fig.6. They are in good agreement, which shows the inversion error has been effectively controlled.

Fig.6 Identified bearing dynamic forces compared with theoretical ones.

Full size|PPT slide

Then the identified bearing dynamic forces are used to further identify the time series of GMS by solving the inverse problem of Eq. (31). The spectra of the meshing stiffness can be obtained by performing FT to the identified time series directly. Furthermore, using the periodic invariance property of meshing stiffness, the waveform in a single-double tooth mesh cycle can be calculated by averaging the segments belonging to different periods in the identified time series. And the trend component in the waveform can also be eliminated in the process.

The identified waveforms and spectra of GMS for healthy, pitting, and wear gear are given in Fig.7. Although the absolute magnitudes of the identified spectra are distorted compared with the theoretical ones, the relative magnitude and frequency components of different fault types match well with the theory. Compared with healthy gear, the magnitudes at wear gear’s meshing frequency and its harmonic components reduce, while pitting gear increases slightly. The spectrum of pitting gear also contains unobvious sideband frequencies. In addition, from the waveforms, the meshing stiffness of pitting and wear gear compared with healthy gear is similar to the simulated ones. The above properties inspired us to use similarity index in the time or frequency domain to realize fault diagnosis and classification instead of relying on the magnitude.

Fig.7 Identified gear meshing stiffness: (a) waveforms of healthy, pitting, and wear gear; and (b) spectra of healthy, pitting, and wear gear.

Full size|PPT slide

4.3 Fault classification based on similarity index

We perform the same process to identify GMS under another condition (the torque is 50 N·m) and use the GMS shown in Fig.7 (the torque is 60 N·m) as the basis of similarity to classify healthy and different fault gears. For fault diagnosis methods that rely on the magnitude, the decrease in load means the decrease in magnitude, resulting in the wrong diagnosis conclusion. However, using the proposed GMSI method with some simple similarity indexes, such as 2-norm and energy ratio, can avoid the above limitation. Considering the properties of waveform and spectrum for different gear fault types mentioned above, 2-norm (defined as S₂) is used to classify healthy and pitting gears, while energy ratio (defined as S_e) is used to classify healthy and wear gears:

(39)S_2={\Vert {\hat{\mathit{k}}}_{{\mathrm{m}}_{\mathrm{f}}}-{\hat{\mathit{k}}}_{{\mathrm{m}}_{\mathrm{b}}}\Vert }_2,

(40)S_{\mathrm{e}}=\sum _{i=1}^q\left([{\mathit{k}}_{{\mathrm{m}}_{\mathrm{f}}}{]}_i-\overline{{\mathit{k}}_{{\mathrm{m}}_{\mathrm{f}}}}\right)/\sum _{i=1}^q\left([{\mathit{k}}_{{\mathrm{m}}_{\mathrm{b}}}{]}_i-\overline{{\mathit{k}}_{{\mathrm{m}}_{\mathrm{b}}}}\right),

where {\mathit{k}}_{{\mathrm{m}}_{\mathrm{f}}} and {\mathit{k}}_{{\mathrm{m}}_{\mathrm{b}}} denote the identified meshing stiffness spectrum and basis spectrum, {\hat{\mathit{k}}}_{{\mathrm{m}}_{\mathrm{f}}} and {\hat{\mathit{k}}}_{{\mathrm{m}}_{\mathrm{b}}} represent the normalization of {\mathit{k}}_{{\mathrm{m}}_{\mathrm{f}}} and {\mathit{k}}_{{\mathrm{m}}_{\mathrm{b}}}, {\left[{\mathit{k}}_{{\mathrm{m}}_{\mathrm{f}}}\right]}_i and {\left[{\mathit{k}}_{{\mathrm{m}}_{\mathrm{b}}}\right]}_i denote the ith element of the q-dimensional vectors {\mathit{k}}_{{\mathrm{m}}_{\mathrm{f}}} and {\mathit{k}}_{{\mathrm{m}}_{\mathrm{b}}}, respectively, and \overline{{\mathit{k}}_{{\mathrm{m}}_{\mathrm{f}}}} and \overline{{\mathit{k}}_{{\mathrm{m}}_{\mathrm{b}}}} indicate their mean value or direct current component, respectively. Note that the definition of similarity index is not unique, and one can define a more powerful similarity index to obtain more robust fault classification results in practice.

The similarity index results of the simulation example are listed in Tab.2. Compared with the basis healthy GMS, S₂ of the identified healthy GMS is the smallest, while compared with the basis pitting GMS, S₂ of the identified pitting GMS is the smallest. The above phenomenon proves that the proposed GMSI method with a 2-norm similarity index can correctly classify gear pitting fault. Furthermore, S_e of the identified GMS for wear gear is about 70% of that for healthy gear, proving that the proposed GMSI method with an energy ratio similarity index can correctly classify gear wear fault.

Tab.2 Similarity index results of the simulation example under 600 r/min

Case of combination	S2 for pitting	Se for wear
Identified healthy GMS & basis healthy GMS	0.3551	0.7748
Identified healthy GMS & basis fault GMS	0.6832	1.0605
Identified fault GMS & basis healthy GMS	0.6542	0.5202
Identified fault GMS & basis fault GMS	0.3382	0.7121

5 Experimental verification

5.1 Experiment settings

This section will verify the potential of using the proposed GMSI method to diagnose real gear faults. The experiments are conducted in a single-stage gearbox test rig, as shown in Fig.8. The pitting and wear fault gears are manufactured to simulate the working state of the test rig for different gear faults, which is achieved by replacing the healthy gear pairs with the manufactured fault gear pairs. The manufactured fault gears are shown on the bottom of Fig.8. Note that the manufactured single tooth wear gear has similar fault characteristics to pitting gear, which is different from all teeth wear gear in the simulation. The gear parameters are given in Tab.3.

Tab.3 Gear parameters of the single-stage gearbox test rig

Gear	Tooth number	Meshing frequency order	Rotating frequency order
Gear on the driving shaft	21	21	1
Gear on the driven shaft	82	21	0.2561

Fig.8 Gearbox test rig and different fault gears.

Full size|PPT slide

Let the test rig operate for a period of time under four operating conditions (including two speeds, 1800 and 2400 r/min, and two loads, 3 and 9 N∙m) and three states (healthy, gear pitting, gear wear), respectively. We use nine triaxial accelerometers to test twelve groups of vibration signals with their serial numbers and locations marked (e.g., 1#) in Fig.8. The sampling frequency is set as 8000 Hz. Without loss of generality, we uniformly use signal segment of 10 s for analysis in the following.

Obviously, the vibration signals measured by accelerometers 1# to 5# contain more complicated structure transfer paths compared with that of 6# to 9#. Therefore, the FRFs from different bearing dynamic forces to locations 1 to 5 are approximately tested by hammer experiment. We hammer four bearing supports in different directions (x, y, z) and test the hammer force and vibration signals simultaneously. Then the complex FRFs of the gearbox housing are obtained by calculating the ratio of these two types of signals in the frequency domain. Taking the FRFs from bearing dynamic force 1x to measuring points 1#‒5# as an example, Fig.9 shows the real and imaginary FRFs test results.

Fig.9 Decoupling frequency response functions (FRFs) test (take excitation at 1x as an example).

Full size|PPT slide

Fig.10 shows the spectra of the measured radial vibration signals at measuring point 1# for 1800 r/min and 9 N∙m condition. Clearly, the full-band signals in the x and y directions and the local band signals marked with fault characteristic frequencies (f_m denotes the meshing frequency, f_si denotes the driving shaft rotating frequency, and f_so denotes the driven shaft rotating frequency) are given simultaneously. From the x- and y-direction vibration signals, the spectra contain many interference frequencies that do not belong to the fault characteristic frequencies. And the signal characteristics of different faults are very similar. We also extract the magnitudes at fault characteristic frequencies in Fig.10. The results show that not all the fault characteristic frequencies can indicate pitting fault, and the fault characteristics of wear fault are almost the same as the healthy state. Therefore, it is tough to diagnose and classify pitting and wear faults only based on the vibration signals.

Fig.10 Spectra of the measured vibration signals for different fault types.

Full size|PPT slide

5.2 Experimental GMSI

Using the vibration signals of measuring points 1#, 2#, and 3#, the bearing dynamic forces are identified first, as shown in Fig.11. The meshing frequency with its harmonic component are clearly observed in each bearing dynamic force spectrum, indicating the inversion error has been effectively controlled.

Fig.11 Spectra of the identified bearing dynamic forces in the experiment.

Full size|PPT slide

To ensure the identified bearing dynamic forces can be further used to identify the time series of GMS, we need to build the dynamic model of the gear−shaft−bearing system in the gearbox of the test rig. The modeling approach is illustrated by Fig.12, where the system mass matrix and damping matrix are easily calculated using the drawings and material properties. The bearing stiffness is estimated by the tested FRFs. Finally, the meshing stiffness matrix containing the GMS information can be identified. Similar to the simulation, the spectra of the meshing stiffness can be obtained by directly performing FT to the identified time series. At the same time, the waveform in a single-double tooth mesh cycle can be calculated by averaging the segments belonging to different periods.

Fig.12 Dynamic modeling of the gear−shaft−bearing system of the test rig.

Full size|PPT slide

Fig.13 shows the identified waveforms and spectra of GMS for healthy, pitting, and wear gear. Limited by non-periodic sampling and noise, the waveforms (Fig.13(a)) have been distorted, although the relative magnitude still matches with the simulation. The pitting and wear faults in the experiment are all single-tooth faults, similar to the simulated pitting gear fault. Therefore, the spectra contain sideband frequencies, as shown in Fig.13(b). Different from simulation, the experimental identification results are interfered by background noise. Thus, only the sideband frequencies with high energy are apparent in the spectra. Interestingly, the frequency components not belonging to the meshing frequency and sideband frequencies are eliminated in the GMS spectra. Clearly, the magnitudes at fault characteristic frequencies are extracted in the following, as shown in Fig.13(c). Compared with the original housing vibration signal, the magnitudes at fault characteristic frequencies in the identified GMS spectra can better indicate the fault. Actually, after transferring through complicated structure paths and being polluted by measurement noise, some weak fault features in the gear meshing vibration have already been very difficult to be distinguished in the housing vibration signal. The GMS identified by the proposed approach reflects the original vibration of the gear pair, and thus more sensitive to gear fault. Considering the properties reflected by Fig.13, the identified GMS spectra are used to achieve fault diagnosis.

Fig.13 Identified gear meshing stiffness in the experiment: (a) waveforms of healthy, pitting, wear gear, (b) spectra of healthy, pitting, wear gear, and (c) fault characteristic frequency magnitude of healthy, pitting, wear gear.

Full size|PPT slide

As a comparison, herein, we use a traditional method to diagnose gear fault. Specifically, the intrinsic mode function (IMF) corresponding to gear fault characteristic frequencies is decomposed using variational mode decomposition (VMD) [11], which is achieved by defining the central frequency of the IMF as gear meshing frequency. Fig.14(a) shows the waveform and time-frequency distribution of the decomposed IMF, taking the Y-direction signal of measuring point 1# as an example. The normalized root mean square (RMS), defined as the ratio of root mean square to the healthy state, is calculated for the decomposed IMFs of the measuring points 1#‒5# (Y direction) and compared with that of the identified GMS, as shown in Fig.14(b). The results show that only part of the measuring points can diagnose pitting fault, and almost all the measuring points are incapable of diagnosing wear fault. In contrast, the identified GMS can simultaneously diagnose pitting and wear fault, which is superior to the original vibration signals. The above results also demonstrate that the proposed GMSI approach is effective and has advantages compared to the traditional VMD fault diagnosis method.

Fig.14 Comparison study of the proposed gear meshing stiffness identification and traditional fault diagnosis approach: (a) waveform and time−frequency distribution of the intrinsic mode function and (b) normalization RMS of different measuring points and the identified gear meshing stiffness (GMS).

Full size|PPT slide

5.3 Influence of measuring location

Based on the proposed GMSI approach, the GMS spectra are identified using the vibration signals of measuring points 2#‒4# and 3#‒5#, respectively. The magnitudes at fault characteristic frequencies are also extracted and given in Fig.15. Overall, the phenomenon reflected by Fig.15 is the same as Fig.13. The influence of measuring location on the identified GMS spectra is reflected by the magnitude of noise relative to the fault characteristic frequencies. The smaller the magnitude of noise relative to the fault characteristic frequencies is, the more sideband frequencies appear and clear. From the perspective of fault diagnosis, the properties of the identified GMS spectra using housing vibration signals at different measuring points are consistent.

Fig.15 Identified gear meshing stiffness in the experiment: fault characteristic frequency magnitude of healthy, wear, and pitting gear using (a) measuring points 2#‒4# and (b) measuring points 3#‒5#.

Full size|PPT slide

To give a more in-depth explanation of the phenomenon, we simply decompose the proposed GMSI approach into a two-step inverse problem shown as:

(41)\{\begin{array}{c}{k}_{\mathrm{m}}\to f\left(\mathit{M},\mathit{D},\mathit{F},{\mathit{F}}_{\mathrm{d}}\right),\\ {\mathit{F}}_{\mathrm{d}}\to f\left({\mathit{x}}_{\mathrm{d}},{\mathit{H}}_{\mathrm{d}}\right).\end{array}

Considering the mass matrix M, the damping matrix D, and the external excitation F are independent of the measuring location, the identified GMS depends on the accuracy of identifying bearing dynamic forces using different vibration signals. For two different sets of measuring points, the vibration transfer equation can be expressed as:

(42)\{\begin{array}{c}{\mathit{x}}_{\mathrm{d}\left(i\right)}={\mathit{H}}_{\mathrm{d}\left(i\right)}{\mathit{F}}_{\mathrm{d}},\\ {\mathit{x}}_{\mathrm{d}\left(j\right)}={\mathit{H}}_{\mathrm{d}\left(j\right)}{\mathit{F}}_{\mathrm{d}},\end{array}

where {\mathit{x}}_{\mathrm{d}\left(i\right)} and {\mathit{x}}_{\mathrm{d}\left(j\right)} denote the ith and jth set of housing vibration signals, respectively, and {\mathit{H}}_{\mathrm{d}\left(i\right)} and {\mathit{H}}_{\mathrm{d}\left(j\right)} denote the corresponding transfer function matrix, respectively. We can derive the relationship between the vibration signals of different measuring points from Eq. (42) as below:

(43){\mathit{x}}_{\mathrm{d}\left(j\right)}={\mathit{H}}_{\mathrm{d}\left(j\right)}{\mathit{H}}_{\mathrm{d}\left(i\right)}^{†}{\mathit{x}}_{\mathrm{d}\left(i\right)},

where the superscript † denotes Moore−Penrose pseudo-inverse. Equation (43) demonstrates that using different vibration signals to identify meshing stiffness is intrinsically associated with the inversion error of {\mathit{H}}_{\mathrm{d}\left(i\right)}^{†}. The inversion error is further reflected as the magnitude of noise relative to the fault characteristic frequencies in the identification results. Therefore, the gear fault can be diagnosed using vibration signals measured on any locations by controlling the inversion error with the proposed method.

5.4 Generalization ability under different speed

As stated in the numerical simulation, using the identified GMS spectra with some simple similarity index can potentially achieve fault classification further. In this subsection and next subsection, we want to discuss the generalization ability of the fault classification based on similarity index under different speeds and loads, respectively. Considering the pitting and wear faults in the experiment and the pitting fault in the simulation belong to single-tooth faults with different severity, hereinafter the similarity index S₂ (defined in Eq. (39)) is used to classify healthy, pitting, and wear.

Fig.16 shows the fault characteristic frequency magnitudes extracted from the identified GMS spectra of different gear faults under another speed (2400 r/min). The reflective phenomenon is still the same as Fig.13 and Fig.15. The amplitude at each fault characteristic frequency together indicates the gear pair has a single-tooth fault, which is better than the original housing vibration signal.

Fig.16 Identified gear meshing stiffness in the experiment: fault characteristic frequency magnitude of healthy, wear, and pitting gear under another speed.

Full size|PPT slide

Then we use S₂ to classify the two single-tooth faults. As a comparison, we calculate S₂ of the original housing vibration signal under 2400 r/min first. The housing vibration signal spectrum under 1800 r/min and 9 N∙m conditions is used as the basis. Considering the meshing frequency is different at different speeds, we perform order tracking to the vibration signals. At the same time, considering the housing vibration signals contain many interference frequencies unrelated to the fault, only [18, 24] range of order band is used to analyze in the following. Same as before, the smaller the similarity index is, the closer it is to the basis. Fairly, we also perform order tracking to the identified GMS and use [18, 24] range of order band to calculate the similarity index. The GMS identified under 1800 r/min and 9 N∙m conditions is used as the basis. The results are listed in Tab.4. It can be seen from Tab.4 that under 2400 r/min, using the original housing vibration signal can only distinguish pitting fault, while using the identified GMS can correctly classify pitting and wear faults.

Tab.4 Similarity index results under different speed using the housing vibration signals or the identified GMS under 1800 r/min and 9 N∙m conditions as the basis

Conditions of vibration signals or GMS	S2
	Healthy housing vibration signal basis	Pitting housing vibration signal basis	Wear housing vibration signal basis	Healthy GMS basis	Pitting GMS basis	Wear GMS basis
Healthy (2400 r/min)	2.8552	2.9223	2.9556	6.0981	3.9826	5.9761
Pitting (2400 r/min)	2.3120	1.4709	2.5707	6.6985	3.8285	6.0991
Wear (2400 r/min)	2.3688	2.0241	2.8034	6.4100	5.9921	5.7597

5.5 Generalization ability under different load

Fig.17 shows the fault characteristic frequency magnitudes extracted from the identified GMS spectra of gear pitting fault and healthy states under another load (3 N∙m). The amplitude at each fault characteristic frequency can still indicates the gear pair has a single-tooth fault.

Fig.17 Identified gear meshing stiffness in the experiment: fault characteristic frequency magnitude of healthy and pitting gear under another load.

Full size|PPT slide

Then we use S₂ to classify healthy gear and pitting fault gear. Similarly, we calculate S₂ of the original housing vibration signal under 3 N·m as a comparison. The housing vibration signal spectrum and the identified GMS under 1800 r/min and 9 N∙m conditions are still used as the basis. We also perform order tracking to them and use [18, 24] range of order band to calculate the similarity index. The results are listed in Tab.5. It can be seen from Tab.5 that under 3 N∙m, using the original housing vibration signal cannot distinguish healthy gear and pitting gear, while using the identified GMS can correctly classify them.

Tab.5 Similarity index results under different load using the housing vibration signals or the identified GMS under 1800 r/min and 9 N∙m conditions as the basis

Conditions of vibration signals or GMS	S2
	Healthy housing vibration signal basis	Pitting housing vibration signal basis	Healthy GMS basis	Pitting GMS basis
Healthy (T = 3 N∙m)	2.5424	2.4864	6.3155	4.0944
Pitting (T = 3 N∙m)	2.3711	1.9553	6.4963	3.9728

The above analysis shows the potential of the proposed GMSI approach to generalize to other speeds and loads. This property is beneficial for fault diagnosis and classification of complicated equipment due to the identified GMS has clear physical meaning. The results also demonstrate that the identified GMS is superior to the original housing vibration signal in fault diagnosis and classification.

5.6 Relationship between spectrum noise and fault severity

In this subsection, we want to discuss another interesting phenomenon found in the identified GMS spectrum: the spectrum noise level corresponds to the fault severity.

The spectrum noise level (defined as the root mean square value of the magnitude of the displayed frequency band) of all the identified GMS spectra used in the above experimental analysis are calculated for healthy and different faults states. Note that the pitting and wear gears used in the experiment have single-tooth faults, and the fault severity of pitting gear is more serious. Fig.18 demonstrates that the spectrum noise level matches the fault severity exactly for every condition. This property makes the identified GMS spectrum noise level have the potential to be further used to describe the machine degradation process.

Fig.18 Spectrum noise level of the identified gear meshing stiffness indicates gear fault severity. Cond1: measuring points 1#‒3#, 1800 r/min, and 9 N∙m; Cond2: measuring points 2#‒4#, 1800 r/min, and 9 N∙m; Cond3: measuring points 3#‒5#, 1800 r/min, and 9 N∙m; Cond4: measuring points 1#‒3#, 2400 r/min, and 9 N∙m; Cond5: measuring points 1#‒3#, 1800 r/min, and 3 N∙m.

Full size|PPT slide

We try to explain this phenomenon from the perspective of fault impact. Considering introducing a noise disturbed fault pulse \delta \left(t\right) with amplitude A in the GMS time series k\left(t\right), then A = 0 for healthy gear pair. Then GMS spectra of healthy gear and different fault gears can be calculated as below:

(44)\{\begin{array}{c}{k}_{\mathrm{h}}\left(\omega \right)=\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}k\left(t\right){\mathrm{e}}^{-\mathrm{i}\omega t}\mathrm{d}t=A_{\mathrm{k}},\\ k_{\mathrm{w}}\left(\omega \right)=\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}k\left(t\right){\mathrm{e}}^{-\mathrm{i}\omega t}\mathrm{d}t+\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}{A}_{\mathrm{w}}\delta \left(t\right){\mathrm{e}}^{-\mathrm{i}\omega t}\mathrm{d}t=A_{\mathrm{k}}+A_{\mathrm{w}},\\ k_{\mathrm{p}}\left(\omega \right)=\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}k\left(t\right){\mathrm{e}}^{-\mathrm{i}\omega t}\mathrm{d}t+\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}{A}_{\mathrm{p}}\delta \left(t\right){\mathrm{e}}^{-\mathrm{i}\omega t}\mathrm{d}t=A_{\mathrm{k}}+A_{\mathrm{p}},\end{array}

where k_{\mathrm{h}}\left(\omega \right), k_{\mathrm{w}}\left(\omega \right), and k_{\mathrm{p}}\left(\omega \right) denote GMS spectrum of healthy gear, wear gear, and pitting gear, respectively, and A_{\mathrm{k}}, A_{\mathrm{w}}, and A_{\mathrm{p}} denote healthy spectrum amplitude and additional amplitude caused by wear and pitting faults. In this paper, the pitting fault is more serious, namely, A_{\mathrm{p}}>A_{\mathrm{w}}. Obviously, we can get A_{\mathrm{k}}+A_{\mathrm{p}}>A_{\mathrm{k}}+A_{\mathrm{w}}>A_{\mathrm{k}}, and thus verify that the above phenomenon is intuitive.

5.7 Discussion

Finally, we want to discuss some noteworthy points of the proposed GMSI approach.

The traditional fault diagnosis methods are mainly amplitude-based diagnosis, thus highly related to the load. GMSI provides a novel technical route for fault diagnosis, i.e., shape-based diagnosis. Shape-based diagnosis has many potential advantages over amplitude-based diagnosis, such as being robust to operating conditions (speed and load) variation, easily achieving fault classification.

Additionally, the identification of GMS is also significant on the following aspects:

(i) First, the GMS can be considered to reflect the state of the gear pair at the source since possible structural transfer effects are eliminated. This is not only beneficial for fault diagnosis of the gearbox with complicated structure but also avoids selecting sensors’ location. One can install sensors at the locations that are convenient for measurement instead of close to bearing support because the GMS can be identified by housing vibration signals measured at any location.

(ii) Second, the GMS identified by the GMSI approach is highly similar to that modeled in the classical dynamics theory, which means that the dynamic model can effectively describe the real mechanical system. The fault simulation based on the dynamic model is expected to save huge manpower and material resources.

(iii) Finally, our experimental results also demonstrate that the spectrum noise level of the identified GMS can indicate gear fault severity, which is expected to further enable quantitative diagnosis.

All of the experimental results in Section 5 demonstrate the above statements. Note that the identified GMS can be regarded as another type of vibration signal. Thus, the research fields associated with the vibration signal are also beneficial to meshing stiffness, e.g., denoising algorithm, similarity index. To our best knowledge, the research on sparse measures in the field of classical fault diagnosis has been very mature, which plays an important role in improving the effect of fault diagnosis. Similarly, we believe that research on similarity indexes can also improve the effect of fault diagnosis using GMS.

6 Conclusions

This paper establishes an analytical mathematical description of the whole transfer process from gear meshing excitation to housing vibration. Then, we propose GMSI approach using housing vibration signals based on the mathematical description. Numerical simulation and experimental results demonstrate the proposed approach can realize gear fault diagnosis better than the original housing vibration signal. With the identified GMS and simple similarity indexes, the proposed approach also has the potential to generalize to other speeds and loads. Some interesting properties are discovered in the identified GMS spectra. For example, the frequency components not belonging to the meshing frequency and sideband frequencies are eliminated, and the spectrum noise level matches with the fault severity exactly. Since the frequency components of the identified GMS are very similar to theoretical modeling, the research results also validate the rationality of using time-varying stiffness to describe the actual gear meshing process. Finally, the identified GMS has a clear physical meaning, which is beneficial for fault diagnosis and classification of complicated equipment.

7 Nomenclature

Abbreviations
DOF	Degree of freedom
FEM	Finite element method
FT	Fourier Transform
FRF	Frequency response function
GMS	Gear meshing stiffness
GMSI	Gear meshing stiffness identification
IMF	Intrinsic mode function
IFT	Inverse Fourier Transform
Variables
\mathit{D}	System damping matrix
f_{ix\mathrm{R}},, f_{ix\mathrm{I}}	Real and imaginary parts of the x component of bearing dynamic force at bearing i, respectively
f_{\mathrm{m}}	Meshing frequency
{\mathit{f}}_{\mathrm{d}}	Intermediate variable of the decoupled bearing dynamic force excitation spectra vector
\mathit{F}	External excitation
{\mathit{F}}_{12}	External load
{\mathit{F}}_{\mathrm{b}}	Bearing dynamic force excitation spectra vector
{\mathit{F}}_{\mathrm{B}}	Bearing force
{\mathit{F}}_{\mathrm{d}}	Decoupled bearing dynamic force excitation spectra vector
{\mathit{F}}_{\mathrm{e}}	Load caused by no-load transfer error of the gear pair
H_{ix\mathrm{R}jx\mathrm{R}},H_{ix\mathrm{I}jx\mathrm{I}}	Frequency response functions
\mathit{H}	Transfer function matrix
{\mathit{H}}_{\mathrm{d}}	Decoupled transfer function matrix
{\mathit{H}}_{\mathrm{T}}	Transition matrix
I1, I2	Transverse moments of inertia of gears 1 and 2, respectively
J1, J2	Polar moments of inertia of gears 1 and 2, respectively
k_{12}	Meshing stiffness of the gear pair
k_{\mathrm{m}i}	Meshing stiffness of gear pair i
k_{xx}, k_{yy}, k_{\mathit{\text{zz}}}	Bearing stiffnesses of translational DOF in x-, y-, and z-direction, respectively
k_{{\theta }_x{\theta }_x}, k_{{\theta }_y{\theta }_y}	Bearing stiffnesses of rotational DOF in x and y directions, respectively
{\mathit{k}}_{\mathrm{b}}	Bearing flexibility vector
{\mathit{k}}_{\mathrm{B}}	Bearing stiffness vector
{\tilde{\mathit{k}}}_{\mathrm{m}}	Temporary variable containing meshing stiffness of all gear pairs
\mathit{K}	System stiffness matrix
{\mathit{K}}_{12}	Meshing stiffness matrix of the gear pair
{\mathit{K}}_i^{\text{s}}	Unit stiffness matrix of the system
{\tilde{\mathit{K}}}_{\mathrm{m}}	Estimated whole meshing stiffness matrix
{\mathit{K}}_{\mathrm{m}i}	Meshing stiffness matrix of gear pair i
m	Number of bearing dynamic force excitation channels
m_{\mathrm{g}}	Mass of gear
m_1, m_2	Masses of gears 1 and 2, respectively
\mathit{M}	System mass matrix
Mg	Mass matrix of a lumped mass model
{\mathit{M}}_{12}	Mass matrix of the gear pair
n	Number of response channels
n_i	Node number
r_{\mathrm{b}1},r_{\mathrm{b}2}	Radii of base circle
S2	2-norm similarity index
Se	Energy ratio similarity index
T1, T2	Torques applied on the input shaft and output shaft, respectively
\mathit{u}	Nodes displacement vector
{\mathit{u}}^{\mathrm{e}}	DOF of the beam model
{\mathit{u}}_{\mathrm{f}}	Vibration response of all nodes in the frequency domain
w_i	Element of the weighted matrix
W	Weighted matrix
x	Translational DOF of beam in x direction
x_1,x_2	Translational DOFs of gears 1 and 2 in x direction, respectively
x_{\mathrm{b}}	Translational DOF of bearing in x direction
x_{jx\mathrm{R}},x_{jx\mathrm{I}}	Real and imaginary parts of the x component of housing vibration response at location j, respectively
\mathit{x}	Response channel spectra vector
{\tilde{\mathit{x}}}_{\mathrm{b}}	Vibration response of all bearing nodes
{\mathit{x}}_{\mathrm{B}}	Displacement vector of bearing
{\mathit{x}}_{\mathrm{d}}	Decoupled response channel spectra vector
{\mathit{X}}_{12}	Generalized coordinate of the gear pair
y	Translational DOF of beam in y direction
y_1, y_2	Translational DOFs of gears 1 and 2 in y direction, respectively
y_{\mathrm{b}}	Translational DOF of bearing in y direction
\mathit{\text{z}}	Translational DOF of beam in z direction
{\mathit{\text{z}}}_1, {\mathit{\text{z}}}_2	Translational DOFs of gears 1 and 2 in z direction, respectively
{\mathit{\text{z}}}_{\mathrm{b}}	Translational DOF of bearing in z direction
\alpha	Weighted factor
{\mathit{\alpha }}_{12}	Projection vector of the gear pair
{\psi }_{12}	Angle between the positive y axis and the meshing surface
Ψ	Projection matrix of all gear nodes
\delta	Numerical error
\omega	Angular frequency
{\theta }_x, {\theta }_y, {\theta }_{\mathit{\text{z}}}	Rotational DOFs of beam in x, y, and z directions, respectively
{\theta }_{x1}, {\theta }_{x2}	Rotational DOFs of gears 1 and 2 in x direction, respectively
{\theta }_{y1}, {\theta }_{y2}	Rotational DOFs of gears 1 and 2 in y direction, respectively
{\theta }_{\mathit{\text{z}}1}, {\theta }_{\mathit{\text{z}}2}	Rotational DOFs of gears 1 and 2 in z direction
{\theta }_{x\mathrm{b}}, {\theta }_{y\mathrm{b}}, {\theta }_{\mathit{\text{z}}\mathrm{b}}	Rotational DOFs of bearing in x, y, and z directions, respectively
{\mathbf{\Theta }}_i	Projection matrix

Acknowledgements

This work was supported by the Basic Research Foundation, China (Grant No. MKF20210013). The authors would like to appreciate the anonymous reviewers for the valuable comments and suggestions to further improve the paper.

References

1 Ding F F , Tian Z G . Integrated prognosis for wind turbine gearbox condition-based maintenance considering time-varying load and crack initiation time uncertainty. International Journal of Reliability Quality and Safety Engineering, 2021, 28(4): 2150024

2 Higgs P A , Parkin R , Jackson M , Al-Habaibeh A , Zorriassatine F , Coy J . A survey on condition monitoring systems in industry. In: Proceedings of the 7th Biennial Conference on Engineering Systems Design and Analysis. Manchester: ASME, 2004, 163–178

3 Meng Z , Shi G X , Wang F L . Vibration response and fault characteristics analysis of gear based on time-varying mesh stiffness. Mechanism and Machine Theory, 2020, 148: 103786

4 Wang D , Tsui K L . Two novel mixed effects models for prognostics of rolling element bearings. Mechanical Systems and Signal Processing, 2018, 99: 1–13

5 Liu K B , Gebraeel N Z , Shi J J . A data-level fusion model for developing composite health indices for degradation modeling and prognostic analysis. IEEE Transactions on Automation Science and Engineering, 2013, 10(3): 652–664

6 Saunders C, Stitson M O, Weston J, Bottou L, Schӧlkopf B, Smola A. Support Vector Machine-Reference Manual. Technical Report CSD-TR-98-03, 2002

7 Huang N E , Shen Z , Long S R , Wu M C , Shih H H , Zheng Q A , Yen N C , Tung C C , Liu H H . The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 1998, 454(1971): 903–995

8 Gilles J . Empirical wavelet transform. IEEE Transactions on Signal Processing, 2013, 61(16): 3999–4010

9 Li N , Huang W G , Guo W J , Gao G Q , Zhu Z K . Multiple enhanced sparse decomposition for gearbox compound fault diagnosis. IEEE Transactions on Instrumentation and Measurement, 2020, 69(3): 770–781

10 Chen S Q , Dong X J , Peng Z K , Zhang W M , Meng G . Nonlinear chirp mode decomposition: a variational method. IEEE Transactions on Signal Processing, 2017, 65(22): 6024–6037

11 Dragomiretskiy K , Zosso D . Variational mode decomposition. IEEE Transactions on Signal Processing, 2014, 62(3): 531–544

12 Daubechies I , Lu J F , Wu H T . Synchrosqueezed wavelet transforms: an empirical mode decomposition-like tool. Applied and Computational Harmonic Analysis, 2011, 30(2): 243–261

13 He D Q , Liu C Y , Jin Z Z , Ma R , Chen Y J , Shan S . Fault diagnosis of flywheel bearing based on parameter optimization variational mode decomposition energy entropy and deep learning. Energy, 2022, 239: 122108

14 Pichika S V V S N , Yadav R , Geetha Rajasekharan S , Praveen H M , Inturi V . Optimal sensor placement for identifying multi-component failures in a wind turbine gearbox using integrated condition monitoring scheme. Applied Acoustics, 2022, 187: 108505

15 Shen H, Zhao F, Mao J G, Zhang C, Hu W. Fault diagnosis of abnormal fuel injection of small aviation piston engine. Journal of Aerospace Power, 2021, 36(4): 861–873 (in Chinese)

16 Tamilselvan P , Wang P F . Failure diagnosis using deep belief learning based health state classification. Reliability Engineering & System Safety, 2013, 115: 124–135

17 Lei Y G , Jia F , Lin J , Xing S B , Ding S X . An intelligent fault diagnosis method using unsupervised feature learning towards mechanical Big Data. IEEE Transactions on Industrial Electronics, 2016, 63(5): 3137–3147

18 Shao S Y , Mcaleer S , Yan R Q , Baldi P . Highly accurate machine fault diagnosis using deep transfer learning. IEEE Transactions on Industrial Informatics, 2019, 15(4): 2446–2455

19 Lu C , Wang Z Y , Qin W L , Ma J . Fault diagnosis of rotary machinery components using a stacked denoising autoencoder-based health state identification. Signal Processing, 2017, 130: 377–388

20 Hoang D T , Tran X T , Van M , Kang H J . A deep neural network-based feature fusion for bearing fault diagnosis. Sensors, 2021, 21(1): 244

21 Wang X , Gu H Y , Wang T Y , Zhang W , Li A H , Chu F L . Deep convolutional tree-inspired network: a decision-tree-structured neural network for hierarchical fault diagnosis of bearings. Frontiers of Mechanical Engineering, 2021, 16(4): 814–828

22 Liu J , Zhou K B , Yang C Y , Lu G L . Imbalanced fault diagnosis of rotating machinery using autoencoder-based SuperGraph feature learning. Frontiers of Mechanical Engineering, 2021, 16(4): 829–839

23 Hong T Z , Wang Z , Luo X , Zhang W N . State-of-the-art on research and applications of machine learning in the building life cycle. Energy and Building, 2020, 212: 109831

24 Yin H , Li Z Z , Zuo J K , Liu H D , Yang K , Li F . Wasserstein generative adversarial network and convolutional neural network (WG-CNN) for bearing fault diagnosis. Mathematical Problems in Engineering, 2020, 2020: 2604191

25 Gao K L , Guo W Y , Yu X C , Liu B , Yu A Z , Wei X P . Deep induction network for small samples classification of hyperspectral images. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, 2020, 13: 3462–3477

26 Hu Z X , Wu D Y , Nie F P , Wang R . Generalization bottleneck in deep metric learning. Information Sciences, 2021, 581: 249–261

27 Wang S N , Xieeryazidan A , Zhang X F , Zhou J X . An improved computational method for vibration response and radiation noise analysis of two-stage gearbox. IEEE Access, 2020, 8: 85973–85988

28 Li Y X, Sun L M. Structural deformation reconstruction by the Penrose−Moore pseudo-inverse and singular value decomposition-estimated equivalent force. Structural Health Monitoring, 2020 (in press)

29 Chen Z , Wang Z , Wang Z H , Chan T H T . Comparative studies on the criteria for regularization parameter selection based on moving force identification. Inverse Problems in Science and Engineering, 2021, 29(2): 153–173

30 Aucejo M . Structural source identification using a generalized Tikhonov regularization. Journal of Sound and Vibration, 2014, 333(22): 5693–5707

31 Zhang Y F , Zhu J S . Damage identification for bridge structures based on correlation of the bridge dynamic responses under vehicle load. Structures, 2021, 33: 68–76

32 Hwang J S , Kareem A , Kim H . Wind load identification using wind tunnel test data by inverse analysis. Journal of Wind Engineering and Industrial Aerodynamics, 2011, 99(1): 18–26

33 Chao M , Hongxing H , Feng X . The identification of external forces for a nonlinear vibration system in frequency domain. Proceedings of the Institution of Mechanical Engineers, Part C: Journal of Mechanical Engineering Science, 2014, 228(9): 1531–1539

34 Sun Q. Aerodynamic load identification of helicopter rotor blade based on hybrid force-modal analysis. Thesis for the Master’s Degree. Nanjing: Nanjing University of Aeronautics and Astronautics, 2017 (in Chinese)

35 Dai H , Long X H , Chen F , Bian J . Experimental investigation of the ring-planet gear meshing forces identification. Journal of Sound and Vibration, 2021, 493: 115844

36 He D A , Chen F , Chao X , Long X H . Numerical calculation and experimental measurement for gear mesh force of planetary gear transmissions. Mechanical Systems and Signal Processing, 2022, 162: 108085

37 Yu X L , Li Z W , He Q B , Yang Y , Du M G , Peng Z K . Gearbox fault diagnosis based on bearing dynamic force identification. Journal of Sound and Vibration, 2021, 511: 116360

38 Yamamoto N , Tsuji T , Murase M , Urakabe K , Yoshimura T . Decomposition of frequency response functions into paths. Mechanical Systems and Signal Processing, 2022, 166: 108457

39 Drača S. Finite element model of a double-stage helical gear reduction. Thesis for the Master’s Degree. Windsor: University of Windsor, 2006

40 Zhang Y M , Wang Q B , Ma H , Huang J , Zhao C Y . Dynamic analysis of three-dimensional helical geared rotor system with geometric eccentricity. Journal of Mechanical Science and Technology, 2013, 27(11): 3231–3242

41 Plunt J . Finding and fixing vehicle NVH problems with transfer path analysis. Sound and Vibration, 2005, 39(11): 12–16

42 Jia Y , Yang Z C , Guo N , Wang L . Random dynamic load identification based on error analysis and weighted total least squares method. Journal of Sound and Vibration, 2015, 358: 111–123

43 Calvetti D, Morigi S, Reichel L, Sgallari F. Tikhonov regularization and the L-curve for large discrete ill-posed problems. Journal of Computational and Applied Mathematics, 2000, 123(1–2): 423–446

44 Choi H G, Thite A N, Thompson D J. Comparison of methods for parameter selection in Tikhonov regularization with application to inverse force determination. Journal of Sound and Vibration, 2007, 304(3–5): 894–917

45 Saad Y , Schultz M H . GMRES: a generalized minimal residual algorithm for solving nonsymmetric linear systems. SIAM Journal on Scientific and Statistical Computing, 1986, 7(3): 856–869

46 Huangfu Y F , Chen K K , Ma H , Li X , Han H Z , Zhao Z F . Meshing and dynamic characteristics analysis of spalled gear systems: a theoretical and experimental study. Mechanical Systems and Signal Processing, 2020, 139: 106640

47 Chen K K , Huangfu Y F , Ma H , Xu Z T , Li X , Wen B C . Calculation of mesh stiffness of spur gears considering complex foundation types and crack propagation paths. Mechanical Systems and Signal Processing, 2019, 130: 273–292