Numerical Modelling of Vibration Responses of Helical Gears under Progressive Tooth Wear for Condition Monitoring

Abstract

Gear wear is a common fault that occurs in a gear transmission system that degrades the operating efficiency and may cause other catastrophic failures such as tooth breakage and fatigue. The progressive wear of a helical gear and its influences on vibration responses are rarely investigated due to the combined effects of the complicated lubrication state and the time-varying characteristic. To fill this gap, a numerical study was put forward to investigate the interactions between gear wear and dynamic response. In this study, an Archard’s wear model with elastohydradynamic lubrication (EHL) effect is adopted to simulate the helical gear wear, which is incorporated with an eight-degree of freedom dynamic model for understanding the gear dynamic at different wear degrees. The wear model shows that the gear wear mainly happens at the gear root due to the relative high slide-to-roll ratio. The dynamic modelling results demonstrate that the wear causes a reduction in time-varying gear mesh stiffness further leads to more vibration. Besides, the simulated vibration responses and experimental validation show that the wear cause increases in the amplitudes of the gear mesh frequency and its harmonics, which can reflect the evolution of progressive gear wear and can be used as monitoring features of gear wear.

1. Introduction

The helical gearbox is the most common apparatus used in various applications such as helicopters, automobiles, wind turbines and marine power trains for transmitting power and speed and altering the direction of motion. Normally, these applications work under high load conditions and variant lubrication states, and suffer from extreme environmental conditions. Thus, gears will inevitably get worn and even fail with some critical faults like tooth breakage and tooth plastic deformations. Gear tooth surface wear is a common early fault that takes place over a long period of service time. The gear wear can cause a deviation of tooth profile from its ideal profile, thus causing a reduction in gear mesh stiffness and further altering the vibration responses of a gear train. However, the effects of the progressive helical gear wear on the dynamic responses are not comprehensively investigated yet, which increases the difficulties in condition monitoring of gear wear.

In the past few decades, many numerical studies have been carried out to simulate the wear process in mechanical contact, among which the commonly used model is the Archard’s wear equation [1,2] proposed based on the theory of asperity contact. For instance, Stanković et al. [3] tried to experimentally determine the Archard’s wear coefficient for sliding bearings. Osman [4] studied the dynamic influences of tooth wear on gear trains. Moreover, there are some other literatures [5,6,7] addressing the modelling of tooth wear, but the wear depth was usually derived from the dry-contact Archard’s wear model, which did not take the lubrication effect into consideration. In fact, for most industrial applications, gears work under lubrication in which the mating surfaces are separated by a thin lubricant film to avoid direct metal-to-metal contact. Therefore, taking into account the lubrication effect is critical necessary to reflect the real wear state. Johnson [8] first proposed the concept of the load sharing in mixed elastohydradynamic lubrication (EHL), which illustrates that the load between the mating surfaces is shared by the lubricant film and asperity contact. Since then, researchers have paid particular attention to the wear occurred in lubrication state. Masjedi [9,10] introduced a method to rapidly calculate the friction coefficient and wear rate in mixed EHL regime taking into account the temperature rise, central film thickness and asperity load ratio. Afterwards, the method was employed on the wear prediction of spur gears by the authors. Based on Masjedi’s study, Wang et al. [11] utilised the method proposed to determine the wear depth for helical gear contact in mixed EHL. Although some studies have considered the wear effects under lubrication state, most of them are focused on the spur gears. Nevertheless, the time-varying characteristic caused by the helical angle brings many difficulties to the numerical simulation. Therefore, the wear simulation of helical gears under lubrication has not been fully studied.

Dynamic model analysis is a very popular method that can study the gear dynamic characteristics from the theoretical point of view. Ding [12] combined a dynamic model with a tooth wear model to investigate the interactions between tooth wear and dynamic behaviours. The results demonstrate that there is a two-way interaction between wear and dynamic behaviours. Liu [7] put forward a prediction model to study the nonlinear correlations between gear wear and dynamic responses of a spur gear set. Ma et al. [13] proposed an improved model for determining the gear mesh stiffness with tip relief, which gives good guidance to study the gear dynamics with wear. These relevant studies demonstrated that the surface wear will alter the gear involute profile, and thus cause a reduction in gear mesh stiffness, changes in contact load and also variation in the transmission errors. Finally, it leads to abnormal vibration and dynamic behaviours. The wear can cause many random and interval variables that increases the uncertainties in dynamics as well [14]. However, at the incipient stage, the wear-related signals are quite weak, which brings many difficulties to monitor the wear state of gears. Researchers are dedicated to realizing the condition monitoring of gear wear in recent years. Some promising works have been carried out to fulfil this goal. Hu [15] developed a state vector extracted from the Time Synchronous Averaging (TSA) vibration signal to evaluate the effect of wear on gear performance. The indicator consists of the sideband information which can reflect the deterioration in gearbox vibration caused by wear. Brethee [16] studied helical gear wear both numerically and experimentally, and the results show that the gear meshing frequency shows a distinguishable increasing trend with the increase in wear severity. Loutas [17] combined the vibration analysis, acoustic emission analysis and oil debris analysis to conduct condition monitoring of rotating machinery. The integration of these three methods increases the diagnostic capacity and reliability. Zhang [18] investigated the modulation signal of gear vibration and used a sideband estimator to monitor the wear process, and a promising result had been demonstrated in this work. Based on the research work presented above, vibration analysis is capable of monitoring the gear wear process, and the meshing frequency and its sideband information contain valuable information for fault diagnosis.

Aiming at the issues addressed above, this study investigates the dynamic responses of helical gear induced by the progressive gear tooth wear through an eight-degree of freedom dynamic model incorporated with an Archard’s wear model under mixed EHL regime. Moreover, the changes in gear mesh stiffness induced by gear wear are comprehensively studied. The simulated vibration responses are validated by an experimental study, which is very beneficial for the condition monitoring of progressive gear wear. This paper is organised as follows. Section 2 introduces the wear depth modelling process based on the Archard’s wear model with EHL effect. Section 3 presents the modelling process of dynamic responses of a helical gear set through an eight-degree of freedom dynamic model, in which the changes in gear mesh stiffness caused by tooth wear were investigated in detail. Section 4 illustrates the simulated dynamic responses at different running hours, which are further verified by an experimental study. Last, several conclusions are drawn in Section 5.

2. Gear Tooth Wear Modelling

2.1. Mesh Characteristics of Helical Gear

For an ideally mounted spur gear pair, when the gear tooth enters the mesh area, the contact line length of the mesh tooth pair is constant and equal to the tooth width. However, for helical gears, the contact line length of one mesh tooth pair is time-varying due to the existence of the helical angle. As shown in Figure 1, the time-varying contact line length of one mesh tooth pair of a helical gear set can be classified into two categories:

$$L_i\left(t\right)=\{\begin{array}{c}{L}_{CD}\csc \beta \gamma \left(t\right), {\epsilon }_{\alpha }\le {\epsilon }_{\beta }\\ b\sec \beta \gamma \left(t\right), {\epsilon }_{\alpha }>{\epsilon }_{\beta }\end{array}$$

(1)

where $b$ is the tooth width, $\beta$ is the helical angle and $\gamma \left(t\right)$ is the periodic piecewise linear function which can be expressed as

$$\gamma \left(t\right)=\{\begin{array}{c} \frac{t^{\prime }}{{\epsilon }_o}, 0\le t^{\prime }\le {\epsilon }_o\\ 1, {\epsilon }_o\le t^{\prime }\le {\epsilon }_t\\ 1-\frac{\left(t^{\prime }-{\epsilon }_t\right)}{{\epsilon }_o}, {\epsilon }_t\le t^{\prime }\le \epsilon \\ 0, \epsilon \le t^{\prime }\le n\end{array}$$

(2)

where ${\epsilon }_0$ is the overlap ratio and ${\epsilon }_t$ is the transverse contact ratio, the total contact ratio $\epsilon ={\epsilon }_0+{\epsilon }_t$. $t^{\prime }=mod\left(t,n{t}_c\right)/t_c$.$n=ceil\left(\epsilon \right)$, $t_c$ is the time duration of the normal pitch which can be calculated as $t_c=\frac{P_b}{{\omega }_g{r}_{bg}}$. ${\omega }_g$, $r_{bg}$ represent the angular velocity and the radius of the base circle of the gear respectively. $P_b$ is the base pitch.

Figure 2 depicts the time-varying contact line length of a pair of helical gears which was selected from the experimental study, for detailed specification see Section 4.4.1. It can be observed that when tooth pair 1 starts to enter the meshing area, its contact line length increases gradually, then remains constant for a period of time and decreases gradually until the mesh point exits the meshing area. As the gear rotates, the gear teeth enter the meshing area in a certain sequence, and the total contact line length can be determined by adding the contact line lengths of all the meshing tooth pairs in the contact region simultaneously.

Compared with spur gears, due to the helical angle, the gear mesh impact generated due to the mesh-in and mesh-out actions can be effectively reduced. However, also due to this time-varying characteristic of the helical gear, the numerical modelling of progressive gear wear is very challenging because there are many time-varying variables, and the wear is also unevenly distributed on the tooth surface due to the loading sharing characteristic. Therefore, to overcome the difficulty induced by the time-varying characteristic, the tooth surface of the helical gear was divided into a grid with i × j segments as shown in Figure 3. By determining the wear depth of each local meshing point P, the wear depth of the helical gear can be simulated. The following section will introduce a wear model to estimate the wear of helical gear in the mixed EHL condition.

2.2. Modelling of Gear Wear under Mixed EHL Regime

In order to make the wear modelling comparable with the experimental validation, the same gear parameter as the gearbox 2 from a real test study was selected for modelling. For detailed parameters, please refer to Section 5. The simulated gearbox comprises of two stages: the first stage works at low speed and the second stage works at high speed. First, the lubrication state of the two stages needs to be determined in order to accurately simulating the wear depth. According to the equation proposed by Dowson and Higginson [19], the lubrication regime can be distinguished by the film thickness ratio which equals to the thickness of the minimum oil film $h_{min}$ divided by the composite surface roughness $R_e$:

$$\lambda =\frac{h_{min}}{R_e}=\frac{1.6{\alpha }_{pv}^{0.6}{\left({\eta }_0{u}_r\right)}^{0.7}{\left(E^{\prime }\right)}^{0.03}{\rho }^{0.43}{w}^{-0.13}}{R_e}$$

(3)

where the definitions of the variables in Equation (3) can be found in

Table 1. Values of the variables for calculating film thickness ratio.

hmin	Minimum film thickness, µm	Re	composite surface roughness, 1.1 × 10−6 m
αpv	Pressure-viscosity factor, 2.03 × 10−8 m2/N	η0	Lubricant viscosity, 0.2883 Pa·s
ur	Rolling velocity, m/s	E	Elastic modulus, 200 GPa
ρ	Equivalent radius, m	E′	$2/E^{\prime }=\left(1-{\nu }_1^2\right)/E_1+\left(1-{\nu }_2^2\right)/E_2$
w	Load per unit length	ν	Poisson ratio, 0.28
kc	Thermal conductivity, W/m·K	k	Wear coefficient, 1 × 10−16 m2/N

, and $u_r, \rho , w$ can be calculated using the gear specification illustrated in

Table 2. Gearbox specifications.

Gearbox	GB1		GB2
	Stage 1	Stage 2	Stage 1	Stage 2
Teeth number, Z	Zr1/Zr2 = 49/55	Zr3/Zr4 = 13/59	Zi1/Zi2 = 59/13	Zi3/Zi4 = 47/58
Helix angle, β (°)	27	13	13	27
Centre distance, a (mm)	74	74	74	74
Normal pressure angle, αn (°)	20	20	20	20
Transverse pressure angle (°)	22.22	20.48	20.48	22.22
Gear tooth width, b (mm)	25	36	36	25
Normal Module, mn (mm)	1.25	2	2	1.25
Transverse contact ratio, εt	1.6514	1.5927	1.5927	1.6518
Overlap ratio, εo	2.8902	1.2889	1.2889	2.8902
Young’s modulus, E (GPa)	200	200	200	200
Poisson ratio, ν	0.28	0.28	0.28	0.28
Surface roughness, σ (µm)	0.8	0.8	0.8	0.8
Vickers hardness, Hd (Pa)	6.865 × 109	6.865 × 109	6.865 × 109	6.865 × 109

and the corresponding working conditions for each stage, respectively.

It can be noticed in Figure 4 that, for the low-speed stage, λ is around 1.48, and for the high-speed stage, $\lambda$ is around 5. According to the Stribeck curve [20], the first stage is working in mixed EHL regime and the second stage is working in full EHL regime. This implies that the first stage is more prone to wear due to the asperity contact in the mixed EHL regime. Thus, the following content will introduce a wear model in mixed EHL regime to determine the wear depth.

According to the publications produced by Archard [1], the traditional Archard’s wear equation which describes the wear of dry contact can be written as

$$V^{\prime }=K\frac{F_n}{H}S$$

(4)

where $V^{\prime }$ is the volume of wear. $S$ is the relative sliding distance between the two meshing points. $K$ is the non-dimensional wear index. $F_n$ is the normal force applied on the contact region. $H$ is the hardness of material. Dividing both sides of Equation (4) by time t yields

$$V_t=K\frac{F_n}{H}{u}_s$$

(5)

where $V_t$ stands for the wear volume per unit time and $u_s$ is the sliding velocity between two mating surfaces. However, the mating gears usually work in the mixed EHL regime, in which the load between the mating surfaces is shared by the asperities and the lubricating film simultaneously. According to research by Masjedi and Khonsari [21], the wear volume per unit time in the mixed EHL regime can be rewritten as

$$V_{\mathrm{lub}}=K\Psi \left(\frac{L_a}{100}\right)\frac{F_n}{H}{u}_s$$

(6)

where $V_{\mathrm{lub}}$ is the wear volume per unit of time in the mixed lubrication, and $L_a$ and $\psi$ are the asperity load ratio of the two contact surfaces and the fractional film defect, respectively, which can be expressed as

$$\Psi =1-\exp \left[-\frac{X}{u_s{t}_0}\exp \left(-\frac{E_a}{R_g{T}_s}\right)\right]$$

(7)

$$\begin{gathered} L_a=0.005W^{-0.408}{U}^{-0.088}{G}^{0.103} \\ \hspace{1em} \times \ln \left(1+4470{\overline{\sigma }}^{6.015}{V}^{1.168}{W}^{0.485}{U}^{-3.741}{G}^{-2.898}\right) \end{gathered}$$

(8)

where ${\rm X}$, $t_0$ $E_a$ are parameters relating to lubrication which can be determined by experiments. $R_a$ is the molar gas constant and $T_s$ is the contact temperature. The parameters of $W, U, G, V$ and $\overline{\sigma }$ are functions of the dimensionless load, speed, modulus, surface hardness and surface roughness, respectively, the expression of which can be found in [21].

Dividing both sides of Equation (6) by $u_s$ yields

$$\frac{V_{\mathrm{lub}}}{u_s}=K\Psi \left(\frac{L_a}{100}\right)\frac{F_n}{H}$$

(9)

Noticed that $V_{\mathrm{lub}}=Ah$ and $p=F_n/A$ ($A$ is the contact area, $h$ is the wear depth per unit of time, $p$ is the average pressure value of all points in a local contact area), Equation (9) can be transformed to

$$\frac{h}{u_s}=K\Psi \left(\frac{L_a}{100}\right)\frac{p}{H}$$

(10)

Therefore, the wear depth $h_{ij}$ of each local contact point $P_{ij}$ on the gear surface is

$$h_{ij}={{\displaystyle \int }}_0^{S}k{\Psi }_{ij}\left(\frac{L_a}{100}\right)p_{ij}dS$$

(11)

where $k$ is the wear coefficient and can be written as $k=K/H$. According to Equation (11), the wear depth $h_{ij}$ at point $ij$ on the gear tooth surface in mixed EHL is determined by the contact pressure $p_{ij}$, the fractional film defect ${\psi }_{ij}$ and the sliding distance $S$.

To determine the contact pressure, the contact between the two mating teeth was regarded as the contact of two cylinders with parallel axes. As stated by the Hertzian contact theory, the half width of the contact area is

$$a_h=\sqrt{\frac{4F_t\rho }{\pi bE}}$$

(12)

where $F_t$ is the tangential force, and $\rho$ and $E$ can be found in Table 1.

Therefore, the contact pressure between the mating point is

$$p=\frac{2F_t}{\pi a_{h}b}$$

(13)

Figure 5 shows the changing process of a pair of mating points of the gear and pinion in the contact area, respectively. As the two contact points have different peripheral velocities, the corresponding sliding distances of the gear and pinion can be, respectively, obtained by the following formulae.

$$S_p=2a_h\frac{u_p-u_g}{u_p}$$

(14)

$$S_g=2a_h\frac{u_g-u_p}{u_g}$$

(15)

where the subscript g, p stands for the pinion and gear respectively. $u_g$ and $u_p$ represent the peripheral velocities of the two mating points on the gears that can be estimated by

$$u_g={\omega }_g\left(\frac{d_g}{2}\sin \alpha +y\right)$$

(16)

$$u_p={\omega }_p\left(\frac{d_p}{2}\sin \alpha -y\right)$$

(17)

where ${\omega }_g$, ${\omega }_p$ denotes the rotational velocity. $d_g, d_p$ denotes the gear pitch diameter. $\alpha$ represents the gear pressure angle. $y$ represents the distance between the pitch point and the mating point.

As stated above, the contact pressure and the relative sliding distance of the mating point can be estimated by Equations (13)–(15), then the wear depth of the mating gears can be obtained by substituting them into Equation (11). The accumulated gear tooth wear after $n$ revolutions can be estimated as

$$h_{ij,\left(n\right)}=h_{ij,\left(n-1\right)}+{{\displaystyle \int }}_0^{S}k{\psi }_{ij}\left(\frac{L_a}{100}\right)p_{ij}d{S}_{ij}$$

(18)

where $h_{ij,\left(n\right)}$ stands for the wear depth at the meshing point $P_{ij}$ on the tooth surface after n revolution, $h_{ij,\left(n-1\right)}$ stands for the wear depth after $n-1$ revolution, and $p_{ij}$ stands for the corresponding pressure at meshing point $P_{ij}$.

According to the wear model, a simulation of helical gear wear was carried out based on the parameters and working conditions of the first stage of the gearbox in the laboratory as depicted in Figure 15, the detailed parameters for the simulation can be found in Table 2. Figure 6a presents the wear depth simulated in the Archard’s wear model in EHL regime. In the direction of the line of action (LOA), when the meshing point moves from the tooth root to the tooth pitch, the wear depth decreases gradually, when it moves from the tooth pitch to the tooth tip, the wear depth increases gradually. Besides, due to the relative high sliding-to-roll ratio, the maximum wear took place close to the tooth root. In the direction of the tooth width, due to changes in the load sharing ratio and meshing tooth pairs, the wear depth shows uneven distribution on the tooth flank. Figure 6b shows the wear depth derived from the dry contact Archard’s wear model; the magnitude of the wear depth is much greater than that derived from the Archard’s wear model in EHL regime. This is because, in the absence of lubricant, the mating surfaces are in direct contact, which will generate greater friction and thus accelerate the wear process. However, in the mixed EHL regime, the load between the mating surfaces is shared by the lubricant film and the asperities, which can effectively reduce the direct asperity contact and significantly reduce the wear rate.

2.3. Wear Depth after Different Running Hours

To study the changes in wear with operating hours, four different times were selected with an interval of 180 h according to the experimental study in Section 4.4.2. Figure 7 shows the maximum tooth wear depth under different running hours. It can be observed that as time passes, the accumulated wear depth increases gradually. In the direction of LOA, the maximum wear depth is around 10 µm, which occurred on the tooth root due to the high relative sliding-to-roll ratio. In the direction of tooth width, the wear distributes the same as the sharing load allocated on the tooth surface. It implies that the high load results in more wear on the gear. According to the wear depth obtained, the following section will introduce a dynamic model to study the influence of the wear on the gear dynamic behaviour.

3. Dynamic Modelling of Helical Gear

3.1. Potential Energy Method for TVMS Calculation

The potential energy method is a well-accepted method for calculating gear mesh stiffness. The method treats the gear as a nonuniform cantilever beam as shown in Figure 8. According to the work in [22], the time-varying mesh stiffness (TVMS) of a spur gear pair includes the Hertzian mesh stiffness $k_h$, bending mesh stiffness $k_b$, shear mesh stiffness $k_s$, axial compressive stiffness $k_a$ and fillet-foundation stiffness $k_f$, which can be determined by

$$U_h=\frac{F^2}{2k_h}, \mathrm{where} k_h=\frac{\pi Eb}{4\left(1-{\nu }^2\right)}$$

(19)

$$U_b=\frac{F^2}{2k_b}={{\displaystyle \int }}_0^d\frac{{\left[F_b\left(d-x\right)-F_a{y}_m\right]}^2}{2E{I}_x}dx$$

(20)

$$U_s=\frac{F^2}{2k_s}={{\displaystyle \int }}_0^d\frac{1.2F_b^2}{2G{A}_x}dx$$

(21)

$$U_a=\frac{F^2}{2k_a}={{\displaystyle \int }}_0^d\frac{F_a^2}{2E{A}_x}dx$$

(22)

$$U_f=\frac{F}{k_f}=\frac{F{\cos }^2{\alpha }_m}{Eb}\left\{L^{*}{\left(\frac{{\mu }_f}{S_f}\right)}^2+M^{*}\left(\frac{{\mu }_f}{S_f}\right)+P^{*}\left(1+Q^{*}{\tan }^2{\alpha }_m\right)\right\}$$

(23)

where $F$ stands for the contact force. $F_a, F_b$ are the radial and tangential forces. ${\alpha }_m$ is the operating pressure angle. $G$ is the shear modulus. $S_n$ is the shear factor.$u_f, S_f, h, d$ and $x$ are depicted in Figure 8. $I_x$ and $A_x$ are the inertia moment and area of the shaded part, respectively, as illustrated in the figure, the calculation of which can refer to the publication [23]. The coefficients $L^{*}, M^{*}, P^{*}, Q^{*}$ can be obtained by the polynomial function proposed by Sainsot [24].

3.2. TVMS Calculation of Helical Gear with Tooth Wear

For helical gear, the gear tooth can be regarded as a combination of multiple spur gears with different angles as shown in Figure 9, thus the total TVMS can be determined by summing up the stiffnesses of all the spur gear slices at the meshing area at every time instant.

When the wear occurred on the gear tooth surface, the gear tooth profile will be altered as shown in Figure 10, thus the tooth thickness of each spur gear slice will be reduced. Correspondingly, the coordinate of the meshing point will be changed as well. As a result, the gear mesh stiffness will be affected. In order to more accurately assess the gear mesh stiffness of worn gears, the involute profile of each spur gear slice was updated with its corresponding wear depth:

$$\{\begin{array}{c}{x}_p=r_{bg}\left(\sin \alpha -\alpha \cos \alpha \right)-h_{p,n}\sin {\alpha }_m\\ y_p=r_{bg}\left(\cos \alpha +\alpha \sin \alpha \right)-h_{p,n}\cos {\alpha }_m\end{array}$$

(24)

where $\alpha$ is the transverse pressure angle at pitch point. The mesh stiffness $k_i^t$ of each spur gear slice at time $t$ can be estimated by the following equation:

$$k_i^t=\frac{1}{k_{hi}^t+k_{bgi}^t+k_{bpi}^t+k_{sgi}^t+k_{spi}^t+k_{agi}^t+k_{api}^t+k_{fgi}^t+k_{fpi}^t}$$

(25)

The total TVMS $k_m^t$ at the time instant $t$ equals to the summation of the stiffnesses of all the spur gear slices at the meshing region:

$$k_m^t={{\displaystyle \sum }}_{i=1}^n{k}_i^t$$

(26)

The damping coefficient of the meshing gears can be calculated as

$${\mathrm{c}}_m^t=2\xi \sqrt{k_m^t{m}_e}$$

(27)

where $m_e$ stands for the equivalent mass of the gears and $\xi$ denotes the viscous damping ratio.

3.3. An Eight-Degree of Freedom Dynamic Model

To study the dynamic responses of the helical gear caused by progressive tooth wear, a dynamic model which has eight degrees of freedom was put forward, as demonstrated in Figure 11. The gearbox used in the experimental study was a two-stage gearbox. However, the study of lubrication state illustrates that the low-speed stage is more prone wearing due to the poor lubrication state, so the dynamic model will only focus on the low-speed stage for simplifying the computation and eliminating the irrelevant factors. In this model, the shaft, gear and pinion are considered rigid, and the ball bearings were treated as linear spring–damping unit. In total, this model includes two torsional motions and six translational motions along the direction of the line of action (LOA), radial direction and axial direction. The governing equations can be written as

$$\{\begin{array}{c}{m}_g{\ddot{y}}_g+c_{gy}{\dot{y}}_g+k_{gy}{y}_g=F_y\\ m_g{\ddot{z}}_g+c_{gz}{\dot{z}}_g+k_{gz}{z}_g=-F_z\\ m_g{\ddot{x}}_g+c_{gx}{\dot{x}}_g+k_{gx}{x}_g=-F_f\\ I_g{\ddot{\theta }}_g=T_d-F_y{r}_{bg}+F_f\left(r_{bg}\tan {\alpha }_m-H\right)\\ m_p{\ddot{y}}_p+c_{py}{\dot{y}}_p+k_{py}{y}_p=-F_y\\ m_p{\ddot{z}}_p+c_{pz}{\dot{z}}_p+k_{pz}{z}_p=F_z\\ m_p{\ddot{x}}_p+c_{px}{\dot{x}}_p+k_{px}{x}_p=F_f\\ I_p{\ddot{\theta }}_p=-T_L+F_y{r}_{bp}-F_f\left(r_{bp}\tan {\alpha }_m+H\right)\end{array}$$

(28)

where I_i represents the inertia moment. m_i stands for the mass of the gear. ${\theta }_i$ represents the angle of rotation. $x_i, y_i, z_i$ denote the translational displacement in x, y, z directions, respectively. c_i and k_i, respectively, denote the damping constant and stiffness constant of the bearing. T_d, T_L are the loads applied on the driving gear and driven pinion. r_bg, r_bp are the radius of the base circle of gear and pinion. H is the distance from the mesh point to the pitch point. F_m, F_f are the dynamic meshing force and friction force that can be estimated by

$$\{\begin{array}{c}{F}_y=k_m^t\cos \beta \left({\tilde{y}}_g-{\tilde{y}}_p\right)+c_m^t\sin \beta \left({\dot{\tilde{y}}}_g-{\dot{\tilde{y}}}_p\right)\\ F_z=k_m^t\sin \beta \left({\tilde{z}}_g-{\tilde{z}}_p\right)+c_m^t\cos \beta \left({\dot{\tilde{z}}}_g-{\dot{\tilde{z}}}_p\right)\\ F_f=\mu \left(t\right)\sqrt{F_y^2+F_z^2}\end{array}$$

(29)

where $k_m^t$, $c_m^t$ can be derived from the Equations (26) and (27). $\mu \left(t\right)$ is the time-varying friction coefficient found in [25]. ${\dot{\tilde{y}}}_g$, ${\dot{\tilde{y}}}_p$, ${\dot{\tilde{z}}}_g$, ${\dot{\tilde{z}}}_p$ are the derivatives of ${\tilde{y}}_g$, ${\tilde{y}}_p$, ${\tilde{z}}_g$, ${\tilde{z}}_p$, which can be calculated as

$$\{\begin{array}{c}{\tilde{y}}_g=y_g+r_{bg}{\theta }_g, {\tilde{y}}_p=y_p-r_{bp}{\theta }_p\\ {\tilde{z}}_g=z_g+{\tilde{y}}_g\tan \beta , {\tilde{z}}_p=z_p-{\tilde{y}}_p\tan \beta \end{array}$$

4. Modelling Results and Experimental Study

4.1. Changes in Gear Meshing Stiffness and Meshing Force

As aforementioned, gear mesh stiffness can be derived using the potential energy method in which the tooth thickness plays an important role to determine the stiffness. Figure 12a demonstrates the total time-varying mesh stiffness which fluctuates periodically due to the changes in total meshing tooth pairs. Besides, note that the mesh stiffness shows a decreasing trend with the growth of the wear depth, which mainly occurs due to the reduction in tooth thickness. Figure 12b–d presents the dynamic meshing force along the y, z, x directions. As can be observed in the zoom-in figures, the meshing force increases gradually with the increase in running hours. This is because the gear wear reduces the mesh stiffness, causing more deformation in contact area, which then leads to an increase in meshing force.

4.2. Raw Spectrum Comparison between Simulated Vibration and Experimental Vibration

The simulated vibration responses are in three dimensions, which will finally transmit to the body of the gearbox through the shaft and bearing eventually, so the vibration of the gearbox body was the superposition of the vibration in these three directions. Hereby, the root mean square (RMS) value was adopted to evaluate the integral vibration f_int of the gearbox body:

$$f_{\mathrm{int}}=\sqrt{\frac{f_{\mathrm{xv}}{}^2+f_{\mathrm{yv}}{}^2+f_{\mathrm{zv}}{}^2}{3}}$$

(30)

where f_xv, f_yv, f_zv stand for the vibration signals in x, y, z directions in the frequency domain. f_int stands for the integral vibration in the frequency domain. Figure 13a exhibits the spectrum of the vibration signal acquired experimentally. The test rig description is detailed in Section 4.4.1. As shown in the figure, the gear meshing frequency and its corresponding harmonics can be clearly observed. Figure 13b shows the simulated vibration spectrum derived from the dynamic model. The two spectra have similar vibration responses that the fundamental meshing frequency is dominated in both of the two figures. Moreover, the harmonics show relatively low amplitudes but an increasing trend with the increase of harmonic order. For the experimental vibration spectrum, there are some other frequency components as it includes the background noise from the bearings, motor and the second stage of meshing gear pair, etc.

4.3. Simulation Results under Different Operating Hours

To get a distinguishable indication of the changes in the vibrations, the amplitudes of the first four harmonics of the meshing frequency were extracted from the spectrum of the integral vibration under the different operating hours. The variation in the amplitude was depicted in Figure 14. As the figure demonstrates, the amplitudes of the meshing frequency and its corresponding harmonics show an increasing trend with the increase of operating time. This is because the wear depth accumulated with the increase of operating time, then leads to the increase in the meshing force in both x, y, z dimensions as previously mentioned, further cause the increase in the gear vibrations. The incremental trend can be an indicator for reflecting the wear evolution process.

4.4. Experimental Results under Different Operating Hours

4.4.1. Description of Experimental Bench

To simulate the process of the progressive wear more precisely and naturally in the laboratory environments, a run-to-fatigue test that lasted more than 800 h was performed on a helical gearbox test rig under variant working conditions, which makes the simulation close to the scenarios of many real-world gearbox applications like wind turbines and helicopters. Figure 15 details the basic structure of the test rig. Two industrial gearboxes were installed back-to-back, which were, respectively, connected with an AC drive motor (15 kW) and a DC load generator (17.5 kW). The geometric parameters of the two helical gearboxes used in this study are shown in Figure 15. For the first stage of gearbox1 (GB1), the gear acts as a driving gear while the pinion is driven. The operating conditions were manipulated by a Variable Speed Drive and Load Controller (VSD). Two accelerometers were mounted on the block surfaces of GB1 and GB2 respectively for collecting the vibration signals. A sixteen-channel data acquisition system was used to collect the vibration, acoustic, oil temperature and encoder signals in 30 s under each working condition.

The test was run for 868 h and then stopped due to the noticeably high peak of GB2 vibration that occurred. During this process, the gearbox was never dissembled to ensure the installation state of the gears and the rig is the same throughout the whole test. After opening the gearbox, the wear was observed to happen in GB2 at the low-speed stage. As presented in Figure 16, the visible wear was generated naturally close to the gear tooth root due to the relative high sliding-to-roll ratio. This is in line with the theoretical study.

4.4.2. Feature Extracted from Experimental Vibration Spectrum

The data collected were divided into two sets: the first set is the data collected before 324 h, during which the gears were working at the break-in stage, which leads to more vibration and noise. However, as time passes, the gear tooth surfaces will be gradually ground and become smooth. The following stage is the progressive wear developing period under mixed EHL. In this period, the vibration and noise vary with the changes in surface topography of the mating tooth. With the development of the wear severity, the vibro-acoustic characteristics of GB2 deteriorate gradually. Therefore, to get more meaningful results, this study will focus on the data collected after 324 h at the operating condition of 70% speed (1039 rpm) and full load.

In order to indicate the development of wear, four sets of data with an equal time interval of 180 h were selected from the data for further analysis. Due to the superiority of time-synchronized averaging (TSA) in terms of noise reduction and signal enhancement, TSA was used to process the vibration signals for suppressing the background noise from bearing, motor and the other gear stage. Figure 17 demonstrates the TSA spectrum of the GB2 vibration signal under different operating hours. The high peaks in the frequency domain represent the meshing frequency f_m3 and its corresponding harmonics. A clear increase in amplitude can be observed from the TSA spectrum, which reflects the deterioration of the gear dynamic due to wear.

Gear mesh impact is the main source of gearbox vibration and noise. Many gear faults can be identified by analysing the characteristic frequencies of a gearbox. Figure 18 depicts the evolution trends of the amplitudes of f_m3 and its harmonics as shown in Figure 18a–d. A noticeable increase in the amplitudes at the frequency of $1\times f_{\mathrm{m}3}$,$2\times f_{\mathrm{m}3}$,$3\times f_{\mathrm{m}3}$ and $4\times f_{\mathrm{m}3}$ can be observed. However, the changes in $f_{\mathrm{m}4}$ and its harmonics did not show any regular variation related to gear wear. First, the reason for the change between the two meshing frequencies $f_{\mathrm{m}3}$ and $f_{\mathrm{m}4}$ is that the wear mainly occurred at the low-speed stage because the relatively low slide-to-roll ratio cannot form a sufficient lubrication film, leading to inadequate lubrication and rapid wear. Moreover, the changes in the amplitudes of $f_{\mathrm{m}3}$ and its harmonics demonstrate that the wear can significantly enlarge the gear mesh impacts due to the reduction in gear contact ratio and the contact pressure angle, thereby further dramatically deteriorates the gear dynamic characteristics. Through the experimental study, the evolution trends of $f_{\mathrm{m}3}$ and its harmonics can effectively indicate the development of progressive gear wear. In addition, the experimental results have a good agreement with the numerical results indicating that the numerical model can give a reliable prediction of gear dynamics, and also providing a reliable theoretical basis for accurate fault diagnosis of progressive gear wear.

5. Conclusions

This study investigated the evolution process of the dynamic responses of a helical gear set induced by progressive gear wear both numerically and experimentally. First, the wear depth of the helical gear set was derived from an Archard’s wear model in mixed EHL condition. Then, the gear dynamic responses were derived from an eight-degree of freedom dynamic model which comprehensively considers the changes in gear meshing position due to wear. Besides, a gear run-to-fatigue test was conducted to validate the numerical results. The results obtained in this study can effectively reflect the progressive wear process of helical gears, and, at the same time, provide a theoretical basis for the condition monitoring of gear trains. Based on this research, the following conclusions can be drawn:

(1) The wear depth simulated from the Archard’s wear model with EHL effect is more reasonable than that derived from the traditional Archard’s wear model, as the lubricant can significantly reduce the friction effect between the mating surfaces. Besides, the wear mainly happens at the region close to the gear root due to the relative high sliding-to-roll ratio. The wear distribution along gear tooth width direction was also identified by load sharing.

(2) The simulation results show that the wear will cause a reduction in gear mesh stiffness and an increase in contact force, which further deteriorates the vibration responses. Correspondingly, the gear mesh frequency gives more valuable information to reflect gear wear process as the amplitudes of gear meshing frequency and its harmonics have increasing trends with the increase of wear depth.

(3) The experimental results show a good consistency with the numerical results. The wear significantly changes the vibration amplitudes of the meshing frequency and its harmonics. The increasing trends of the amplitudes can be effective indicators for condition monitoring the progressive wear in gear train. By employing appropriate advanced signal processing method, the gear tooth progressive wear can be identified using the changes in the amplitudes of the gear meshing frequency, which further can be used for diagnosing the gear wear.