Analysis of Time-Varying Mesh Stiffness and Dynamic Response of Gear Transmission System with Pitting and Cracking Coupling Faults

Abstract

The gear transmission system is an important part of the mechanical system, so it is essential to judge its running state accurately. To solve the difficult problem of identifying the components of coupling faults, this paper derives the calculation method of gear time-varying mesh stiffness for coupling faults of pitting and cracking based on the energy method and considering the coupling between teeth, establishes the dynamics model of two-stage gear transmission system with coupling faults and studies the influence of coupling faults on gear time-varying mesh stiffness and dynamic characteristics. The accuracy of the proposed method is verified by experiments. The results show that both pitting and cracking can lead to a reduction in mesh stiffness. The stiffness of pitting will fluctuate irregularly due to the influence of pitting on the tooth surface, while the stiffness of cracked teeth is relatively smooth. The coupling fault stiffness is dominated by more serious faults. By analyzing the periodic impact components in time domain and the sideband components around the harmonics in frequency domain the faulty gears in the transmission system can be distinguished. It provides an effective reference for the diagnosis of faulty gears.

1. Introduction

Gears are important components of transmission systems and have been widely used in various items of large mechanical equipment, industrial production and military industries due to their precise transmission ratios, long service life and high reliability. However, they may fail over time due to the inevitable fatigue that occurs in gear systems under heavy loads and harsh working conditions [1]. Inadequate lubrication and overloading can cause deformation, which can lead to wear, cracking, pitting, spalling and even tooth breakage, directly affecting the safe operation of equipment. Particularly in some large-scale engineering fields, gear failures often cause huge economic losses and threaten the safety of people and property [2]. Therefore, it is particularly important to study the mechanism of gear failure [3,4]. The more common forms of gear failure are cracking and pitting [5]. Cracking occurs mainly at the tooth root during the meshing process, while pitting occurs on the tooth surface. However, whether pitting or cracking occurs, it will damage the surface material, reduce the gear teeth’s effective cross-sectional area, change the gears’ meshing stiffness and ultimately affect the dynamic characteristics of the gear system [6].

Time-varying mesh stiffness of gears is one of the main internal excitations in gear dynamics [7]. It is a time-varying phenomenon caused by changes in the number and location of tooth contacts [8]. In gear fault diagnosis, the square-wave method, the potential energy method and other methods are commonly used to calculate the time-varying mesh stiffness of gears. In the case of gear cracks, the crack path is usually simplified as a straight line [9] or a slightly curved line starting from the tooth root [10,11]. Most studies assume that tooth cracks extend across the entire width of the tooth and have a constant crack depth. For healthy gears running continuously, gear mesh stiffness is a periodic function. F. Chaari et al. [12,13] used a square wave to approximate the time-varying mesh stiffness of the gear. The period of a square wave is called the mesh period, which is equal to the duration of one rotation divided by the number of teeth. The square wave can reflect changes in the number of tooth contacts, but ignores changes in the tooth contact position, and the size of the square wave is essentially estimated based on personal experience. Hence, the accuracy of the gear time-varying stiffness given by this method is not high. Therefore, many researchers have improved the analytical models for stiffness calculation. Weber [14] derived a comprehensive deformation calculation method using energy integration, including tooth bending, shear and compression. Cornell [15] further proposed a numerical integral method including root fillet and elastic deformation of tooth foundation based on Weber’s research. Yang and Lin [16] established a potential energy method based on the potential energy principle. They considered the total energy stored in a pair of gears as the sum of Hertz contact energy, bending energy and axial compression energy, which correspond to the Hertz contact stiffness, bending stiffness and axial compression stiffness, respectively. Later, Tian [17] added another energy term called shear energy, corresponding to the shear stiffness, and then used this method to derive the mesh stiffness of gears with cracks. This method calculates the mesh stiffness of gears with faults more efficiently and has been widely used. Since previous scholars did not consider the deformation of the gear body when calculating the stiffness, i.e., the gear body was assumed to be rigid, Sainsot et al. [18] proposed an empirical formula for tooth deformation caused by fillet-foundation deflection, which has been widely used to study the effects of various factors such as wear, cracking, stripping and pitting on mesh stiffness. Ankur Saxena [19] studied the effects of various tooth cracking conditions on the system’s modal characteristics and frequency response function. The modal and frequency response characteristics induced by gear cracking were compared. Kramberger et al. [20] indicated that cracks mostly initiated at the point of the maximum principal stress in the tensile side of a gear tooth. Wan et al. [21] proposed a modified model that considers the relative position of the root and base circles, and optimized and corrected the stiffness calculation. The above model ignores the gear structural coupling deformation caused by the gear body when calculating the double-tooth contact area. Xie et al. [22] proposed a time-varying mesh stiffness calculation method that considers the structural coupling effect between teeth and verified its effectiveness by finite element analysis. This further improved the accuracy of the energy method. Regarding pitting failure, scholars’ research mainly focuses on improving the similarity between the pitting analysis model and the actual pitting failure model. Ankur Saxena [23] proposed an analytical method to calculate the TVMS of the spur gear for different spall shapes, size and location considering sliding friction. Bilal El Yousfi [24] double-discretized the tooth surface to consider the variation of defect depth in the width and length direction of the gear teeth, and proposed a method to calculate meshing stiffness of spur gears based on potential energy. Cheng et al. [25] proposed a rectangular pitting model and put forward a method to evaluate the degree of pitting damage. Chen et al. [26] regarded the pitting pit as cylindrical. They considered the influence of crack and pitting pit deterioration on the time-varying mesh stiffness of gears when both existed on the same gear, which was verified by finite element analysis. Luo et al. [27,28] proposed a pitting pit model based on ellipsoidal geometry and verified the accuracy of the model using the finite element method. Lei et al. [29] proposed a pitting pit distribution model based on one-dimensional Gaussian distribution and simulated the variation of gear mesh stiffness under different degrees of failure. Chen et al. [30] derived the comprehensive mesh stiffness affected by pitting using a two-dimensional Gaussian distribution model and further analyzed the conditions of multiple tooth surface pitting.

In summary, the above research has focused mainly on the mechanism and detection of single faults in a single gear pair. However, in important mechanical and industrial fields, various coupling faults may occur in the gear transmission system under heavy loads and insufficient lubrication, among which pitting and cracking are most likely to occur. To address this problem, this paper proposes a stiffness calculation model for coupling, pitting and cracking faults based on the energy method considering the coupling effect between teeth, analyzes the influence of different degrees and types of faults on the time-varying mesh stiffness and establishes a two-stage gear transmission system dynamic model containing single faults and coupling faults based on the finite element method. The effects of single and coupling faults on the dynamic characteristics of the transmission system are analyzed, and the accuracy of the model is verified through experiments.

2. Calculation of Gear Time-Varying Mesh Stiffness

The time-varying meshing stiffness, as an internal excitation, is particularly important in the calculation of the dynamic characteristics of the gear. Therefore, accurate calculation of the meshing stiffness of the gear is conducive to accurate dynamic simulation and dynamic characteristics of the gear system.

This chapter calculates the time-varying mesh stiffness of healthy and faulty gears, and the parameters of gears are shown in

Table 1. Gear parameter.

Parameters	Drive Gear	Driven Gear
Number of teeth	36	90
Pressure angle (°)	20	20
Face width (mm)	12	12
Modulus (mm)	1.5	1.5
Rotary inertia (kg·m2)	0.000380	0.003492
Mass (kg)	0.1836	1.3114

.

2.1. Calculation of Meshing Stiffness of Normal Gear

As a mathematical method, the potential energy method is combined with the mechanics of materials to divide the object and integrate each part to obtain accurate results. Therefore, this method is also used to calculate the stiffness of gears. According to numerous studies, it is believed that the energy in teeth can be divided into four parts. The shear, bending, axial and Hertzian energies are replaced by U_s, U_b, U_a and U_h, respectively. The above four kinds of energy formula calculation are described as follows [26]:

$$\left\{\begin{array}{l}{U}_h=\frac{F^2}{2k_h}\\ U_b=\frac{F^2}{2k_b}={{\displaystyle \int }}_0^d\frac{{\left[F_b\left(d-x\right)-F_{a}h\right]}^2}{2E{I}_x}dx\\ U_s=\frac{F^2}{2k_s}={{\displaystyle \int }}_0^d\frac{1.2F_b^2}{2G{A}_x}dx\\ U_a=\frac{F^2}{2k_a}={{\displaystyle \int }}_0^d\frac{F_a^2}{2E{A}_x}dx\end{array}\right.$$

(1)

$$I_x=\frac{1}{12}{\left(2h_x\right)}^{3}L$$

(2)

$$A_x=\left(2h_x\right)L$$

(3)

where F represents the acting force by the meshing tooth in the contact point. F_a, F_b are radial and tangential forces, G, E, L, represent shear modulus, Young’s modulus and the tooth width, respectively. I_x, A_x are the moment of inertia and cross-sectional area. the other parameters are shown in Figure 1.

According to the principle of gear, the relative positions of the dedendum circle and the base circle are not fixed, which means the integral interval needs to be handled in different cases [21]. The radii of the base circle and dedendum circle can be expressed as follows:

$$r_b=\frac{mz}{2}\cos \left(\theta \right), r_f=\frac{mz}{2}-\left(h_a^{*}+c^{*}\right)m$$

(4)

where m, z and $\theta$ represent module, number of teeth and pressure angle, respectively. $c^{*}$ and $h_a^{*}$ are tip clearance and addendum coefficients. Taking the standard gear as an example, $h_a^{*}=1$, $c^{*}=0.25$, $\theta =20°$, when $r_b=r_f$, the number of teeth z ≈ 42.

If the number of teeth is less than 42, that is r_b > r_f, as shown in Figure 1a, the bending energy of teeth can be written as:

$$\begin{array}{l}{k}_b=1/({{\displaystyle \int }}_{-{\alpha }_1}^{{\alpha }_2}\{3({\alpha }_2-\alpha )\cos \alpha [1+\cos {\alpha }_1(({\alpha }_2-\alpha )\sin \alpha -\\ \cos \alpha ){]}^2/2EL{H}^3\}\mathrm{d}x+{{\displaystyle \int }}_0^{r_b-r_f}\{[\cos {\alpha }_1(d-x_1)-\\ \sin {\alpha }_{1}h{]}^2/E{I}_{x_1}\}\mathrm{d}{x}_1)\end{array}$$

(5)

Similarly, the shear stiffness k_s and axial compression stiffness k_a can be obtained as follows:

$${\mathrm{k}}_s=1/({\displaystyle {\int }_{-{\alpha }_2}^{{\alpha }_1}\frac{(1+\nu )({\alpha }_2-\alpha )\cos \alpha {\cos }^2{\alpha }_1}{ELH}}\mathrm{d}\alpha +{\displaystyle {\int }_0^{r_b-r_f}\frac{1.2{\cos }^2{\alpha }_1}{G{A}_x}d{x}_1})$$

(6)

$${\mathrm{k}}_a=1/({\displaystyle {\int }_{-{\alpha }_2}^{{\alpha }_1}\frac{\left({\alpha }_2-\alpha \right)\cos \alpha {\sin }^2{\alpha }_1}{2ELH}}\mathrm{d}\alpha +{\displaystyle {\int }_0^{r_b-r_f}\frac{{\sin }^2{\alpha }_1}{E{A}_x}d{x}_1})$$

(7)

where I_x denotes the moment of inertia of the section at a distance x₁ from the base circle. Combined with the geometric relation of involute, d, x and h_x can be expressed as:

$$\left\{\begin{array}{l}d=r_b\left[\cos {\alpha }_1+\left({\alpha }_1+{\alpha }_2\right)\sin {\alpha }_1-\cos {\alpha }_2\right]\\ x=r_b\left[\cos \alpha +\left(\alpha +{\alpha }_2\right)\sin \alpha -\cos {\alpha }_2\right]\\ h=r_b\sin {\alpha }_2\\ H=\sin \alpha +({\alpha }_2-\alpha )\cos \alpha \end{array}\right.$$

(8)

If the dedendum circle is greater than base circle, that is $r_b<r_f$, as shown in Figure 1b, the bending energy of the tooth can be written as:

$$k_b=1/({{\displaystyle \int }}_{-{\alpha }_1}^{{\alpha }_2}\{3({\alpha }_2-\alpha )\cos \alpha {[1+\cos {\alpha }_1(({\alpha }_2-\alpha )\sin \alpha -\cos \alpha )]}^2/2EL{H}^3\}\mathrm{d}x)$$

(9)

$${\mathrm{k}}_s=1/({\displaystyle {\int }_{-{\alpha }_2}^{{\alpha }_1}\frac{(1+\nu )({\alpha }_2-\alpha )\cos \alpha {\cos }^2{\alpha }_1}{ELH}}\mathrm{d}\alpha )$$

(10)

$${\mathrm{k}}_a=1/({\displaystyle {\int }_{-{\alpha }_2}^{{\alpha }_1}\frac{\left({\alpha }_2-\alpha \right)\cos \alpha {\sin }^2{\alpha }_1}{2ELH}}\mathrm{d}\alpha )$$

(11)

The Hertz contact stiffness is expressed as [29]:

$$k_h=\frac{EL\pi }{4(1-{\mu }^2)}$$

(12)

According to the theory of Muskhelishvili, in addition to the tooth deformation, the deflection of the fillet-foundation has a great influence. In this case, the tooth is regarded as a rigid body and the fillet-foundation as an elastic body. The flexible deformation stiffness of the fillet-foundation can be calculated by:

$$\frac{1}{k_f}=\frac{{\cos }^2{\alpha }_1}{EL}\{L^{*}{\left(\frac{u}{S_f}\right)}^2+M^{*}\left(\frac{\mathrm{u}}{S_f}\right)+P^{*}(1+Q^{*}{\tan }^2{\alpha }_1)\}$$

(13)

where the coefficients L*, M*, P*, Q* are the same as in [18]. $u$ and $S_f$ are shown in Figure 2.

In most previous calculation models, the stiffness of healthy teeth in the double teeth-meshing region was commonly calculated by directly summing the corresponding single teeth-meshing stiffness. Such a method ignores the influence of the two engaged teeth sharing the same body, which makes the stiffness of the double teeth-meshing region significantly larger compared to the results of the finite element method (FEM) [31]. This will directly affect the results of the dynamic solution. The structural coupling between the teeth in the double teeth-meshing region is schematically shown in Figure 2.

Referring to [22], the structural coupling stiffnesses k_f21 and k_f12 can be expressed, respectively, as:

$$\begin{array}{l}\frac{1}{k_{f21}}=\frac{cos{\alpha }_{1}cos{\alpha }_2}{EL}\{L_2\left(\frac{u_1{u}_2}{S_f^2}\right)+[M_{2}tan{\alpha }_2+P_2]\left(\frac{u_1}{S_f}\right)+\left[Q_{2}tan{\alpha }_1+R_2\right]\left(\frac{u_2}{S_f}\right)+\\ \left[S_{2}tan{\alpha }_1+T_2\right]tan{\alpha }_2+U_{2}tan{\alpha }_1+V_2\}\end{array}$$

(14)

$$\begin{array}{l}\frac{1}{k_{f12}}=\frac{cos{\alpha }_{1}cos{\alpha }_2}{EL}\{L_1\left(\frac{u_1{u}_2}{S_f^2}\right)+[M_{1}tan{\alpha }_1+P_1]\left(\frac{u_2}{S_f}\right)+\left[Q_{1}tan{\alpha }_2+R_1\right]\left(\frac{u_1}{S_f}\right)+\\ \left[S_{1}tan{\alpha }_2+T_1\right]tan{\alpha }_1+U_{1}tan{\alpha }_2+V_1\}\end{array}$$

(15)

where 1/k_f21 denotes the displacement of meshing point 1 under the condition that F₁ = 0 and F₂ is the unit force. 1/k_f12 is the displacement of meshing point 2 when F₂ = 0 and F₁ is the unit force. Parameters such as L_i, M_i, P_i, Q_i, R_i, S_i, T_i, U_i and V_i (i = 1,2) are the same as those in [22], u₁, u₂ and s_f are the same as mentioned above.

Therefore, the comprehensive mesh stiffness can be obtained from Equation (16):

$$k=\left\{\begin{array}{l}\frac{1}{\frac{1}{k_{pf1}}+\frac{1}{k_{pt1}}+\frac{1}{k_{h1}}+\frac{1}{k_{gt1}}+\frac{1}{k_{gf1}},} \mathrm{Single} \mathrm{tooth} \mathrm{meshing} \mathrm{area} \\ {\displaystyle \sum _{i=1}^2{k}_i}, \mathrm{Double} \mathrm{tooth} \mathrm{meshing} \mathrm{area}\end{array}\right.$$

(16)

where $k_i=\frac{1}{\frac{F_j}{F_i{k}_{pfij}}+\frac{1}{k_{pfi}}+\frac{1}{k_{pti}}+\frac{1}{k_{hi}}+\frac{1}{k_{gti}}+\frac{1}{k_{gfi}}+\frac{F_j}{F_i{k}_{gfij}}}\left(i\ne j\right)$, $k_{hi}\left(i=1,2\right)$ denotes the Hertz contact stiffness. $k_{nti}\left(i=1,2;n=p,g\right)$ is the stiffness of the tooth part. It includes bending stiffness, shear stiffness and axial compression stiffness. $k_{nfi}\left(i=1,2;j=p,g\right)$ is the flexible deformation stiffness of the fillet-foundation. $k_{nfij}(i,j=1,2,i\ne j;n=p,g)$ is the coupling stiffness of the structure between teeth.

2.2. Calculation of Time-Varying Mesh Stiffness of Gear Pair with Fault

Due to the concurrent and coupling effects of faults, cracking and pitting are most likely to occur simultaneously under heavy load and poor lubrication. In [26] the authors study pitting and cracking on the same tooth, but the gear failure has the driving wheel and driven wheel fault coupling. Therefore, the influence of pitting and cracking on the time-varying meshing stiffness of the driving and the driven wheel in the transmission system is studied in this paper. The fault distribution is shown in Figure 3.

2.2.1. Calculation of Time-Varying Meshing Stiffness of Gear with Root Crack

It is shown in [32] that the crack propagation path is smooth, continuous and, in most cases, fairly straight, with only a slight curvature. Therefore, in this study, it is assumed that the cracked gear is a cantilever beam model, the crack exists at the dedendum along the tooth width direction through the crack, the intersection angle between the crack and the gear center line is a constant ν and the tooth profile curve remains intact [21]; the details are shown in Figure 4.

The number of teeth of the gear where the crack is located is 90 teeth, so the case that the dedendum circle is larger than the base circle is considered. Because the crack does not change the length of the gear contact line, the Hertz contact stiffness does not change. In addition, the gear can still bear axial compression force after the root cracks, so it can be considered that the axial compression stiffness will not change. In this paper, it is assumed that there are two stages of crack growth as shown in Figure 4. Figure 4a shows the initial stage of crack growth, and Figure 4b shows the more serious stage after crack growth [21].

As shown in Figure 4a, the integration area can be divided into normal and fault areas. The moment of inertia and the cross-sectional area of the gear in the region from the root to g_c away from the root are changed due to the crack in the fault region. The expression is as follows [21]:

$$A_{\mathrm{x}}=\left\{\begin{array}{l}\left(h_c+h_{\mathrm{x}}\right)L \left(x\le g_{\mathrm{c}}\right)\\ 2h_{\mathrm{x}}L \left(x>g_{\mathrm{c}}\right)\end{array}\right.$$

(17)

$$I_{\mathrm{x}}=\left\{\begin{array}{l}{\left(h_c+h_{\mathrm{x}}\right)}^{3}L/12 \left(x\le g_{\mathrm{c}}\right)\\ 2h_{\mathrm{x}}{}^{3}L/3 \left(x>g_{\mathrm{c}}\right)\end{array}\right.$$

(18)

In this case, the axial compression stiffness is calculated in the same way as Equation (11), while the bending stiffness and shear stiffness can be expressed as [21]:

$$\begin{array}{l}{k}_{\mathrm{b}}=1/({\displaystyle {\int }_{-{\alpha }_{\mathrm{g}}}^{{\alpha }_2}\{12{\{1+\cos {\alpha }_1\left[\right({\alpha }_2-\alpha )\sin \alpha -\cos \alpha \left]\right\}}^2({\alpha }_2-}\alpha )\cos \alpha /\{EL[\sin {\alpha }_2-\left(q/r_{\mathrm{b}}\right)\sin \nu +\sin \alpha +\\ \left({\alpha }_2-\alpha \right)\cos \alpha {]}^3\left\}\right\}\mathrm{d}a+{\displaystyle {\int }_{-{\alpha }_1}^{-{\alpha }_{\mathrm{g}}}\left\{3\right\{1+\cos {\alpha }_1\left[\right({\alpha }_2-}\alpha )\sin \alpha -\cos \alpha {\left]\right\}}^2({\alpha }_2-\alpha )\cos \alpha /\{2EL[\sin \alpha +\\ ({\alpha }_2-\alpha )\cos \alpha {]}^3\left\}\right\}\mathrm{d}\alpha )\\ \end{array}$$

(19)

$$\begin{array}{l}{k}_{\mathrm{s}}=1/({\displaystyle {\int }_{-{\alpha }_{\mathrm{g}}}^{{\alpha }_{{}_2}}\{2.4(1+\nu \left)\right({\alpha }_2-\alpha )\cos \alpha {\cos }^{{}^2}{\alpha }_1/} (EL(\sin {\alpha }_2-(q/r_{\mathrm{b}})\sin \nu +\sin \alpha +({\alpha }_2-\\ \alpha )\cos \alpha \left)\right)\}\mathrm{d}a+{\displaystyle {\int }_{-{\alpha }_{{}_1}}^{-{\alpha }_{\mathrm{g}}}\left\{2.4\right(1+\nu \left)\right({\alpha }_2-}\alpha )\cos \alpha {\cos }^{{}^2}{\alpha }_1/\{EL[\sin \alpha +({\alpha }_2-\alpha )\cos \alpha \left]\right\}\}\mathrm{d}a)\end{array}$$

(20)

As shown in Figure 4b, the cross-sectional area and moment of inertia of the gear can be expressed in the following form:

$$A_{\mathrm{x}}=\left(h_c+h_{\mathrm{x}}\right)L$$

(21)

$$I_{\mathrm{x}}={\left(h_c+h_{\mathrm{x}}\right)}^{3}L/22$$

(22)

In this case, bending stiffness and shear stiffness can be expressed as:

$$\begin{array}{l}{k}_{\mathrm{b}}=1/({\displaystyle {\int }_{{-\mathsf{\alpha }}_{\mathrm{g}}}^{{\mathsf{\alpha }}_2}\{12{\{1+\cos {\alpha }_1\left[\right({\alpha }_2-\alpha )\sin \alpha -\cos \alpha \left]\right\}}^2({\alpha }_2-} \alpha )\cos \alpha /\{EL[\sin {\alpha }_2-(q/r_{\mathrm{b}})\sin \nu +\sin \alpha +\\ ({\alpha }_2-\alpha )\cos \alpha {]}^3\left\}\right\}\mathrm{d}\alpha )\end{array}$$

(23)

$$k_{\mathrm{s}}=1/({\displaystyle {\int }_{{-\mathsf{\alpha }}_{\mathrm{g}}}^{{\mathsf{\alpha }}_2}\left\{2.4\right(1+\nu \left)\right({\alpha }_2-\alpha )\cos \alpha {\cos }^2{\alpha }_1/}\{EL[\sin {\alpha }_2-(q/r_{\mathrm{b}})\sin \nu +\sin \alpha +({\alpha }_2-\alpha )\cos \alpha \left]\right\}\}\mathrm{d}a$$

(24)

The calculation results of the time-varying meshing stiffness of the gear with crack fault are shown in Figure 5. The results show that the meshing stiffness of the gear gradually decreases with the increase in the crack degree. The degree increases gradually, because the crack causes the change of the inertia moment and cross-sectional area of the gear.

On the other hand, it can be observed from the figure that the stiffness greatly decreases at the initial position and decreases less with the subsequent meshing. This is because the location of the faulty gear is the driven wheel, and its meshing process is from the addendum to the dedendum. This indicates that the first double-tooth region when the faulty tooth engages in meshing has the greatest impact on the vibration of the system after the root crack fault occurs in the driven wheel. With the meshing of the driven wheel, the bending moment generated by the meshing force in the vertical direction gradually decreases, and the bending potential energy gradually decreases. The stiffness reduction caused by the reduction in cross-sectional area and moment of inertia caused by the crack is relatively reduced.

2.2.2. Pitting Gear Modeling

Due to too thin oil conditions, friction between material surfaces leads to surface defects or cracks, which transform into pitting. With the increased gear running time, the pitting degree further deteriorates from a healthy condition to a severe degree. In addition, when the contact point coincides with the gear pitch line, pure rolling will occur, and it is difficult to form lubricating oil films. If the lubrication is insufficient or the load is too large, the gear pitch line will naturally produce fatigue particles [33]. Therefore, pitting due to fatigue usually occurs in the banded area below the pitch line, which indicates that pits are distributed below the pitch line along the height of the tooth surface. The actual shape of pitting is usually irregular, but to facilitate the establishment of the model, this paper assumes that the shape of pitting is a regular circle [7,34]. Many previous models were designed to simulate the distribution of pits and consider the propagation of pits. As the surface cracks spread, a series of small pits are formed. After this, the pitted gear continues to mesh with other healthy gears. During the rest of the run, sliding during gear meshing squeezes oil or small particles into the surface crack, accelerating the crack propagation or pitting. As a result, these previously created pits continue to spread. Therefore, to simulate the real situation, this study simulates the gradually increasing size of these pits instead of assuming fixed pits. Moreover, in formerly healthy areas, new potholes appear over time. This phenomenon is also considered in this paper.

In summary, in this paper, the pitting distribution on gears is modeled as a two-dimensional random variable and x_i and y_i are the coordinate values of the i-th pitting pit in the direction of the width and height of the tooth surface, respectively. The relationship is shown below:

$$x_i~N(\mu ,{\sigma }^2)$$

(25)

$$y_i~U(0,L_i)$$

(26)

where L_i is the length of pitting spread along the tooth width. Besides, for the normal distribution function in Equation (25), the 3σ criterion is used to describe the distribution areas of pits in this paper. Therefore, μ and σ can be calculated as follows:

$$\mu =x_{\mathrm{p}}-\delta ,\sigma =\frac{x_{\mathrm{p}}-x_{\min }-\delta }{3}$$

(27)

where x_p is the coordinate value of the pitch circle in the x direction, x_min is determined by the boundary of meshing area and d indicates that the distribution of pits usually has a concentration below the pitch line, illustrated in Figure 6.

In this paper, according to Figure 7, pitting is further analyzed according to its surface integral in the tooth surface area into three damage degrees: slight pitting, moderate pitting and severe pitting, as shown in

Table 2. Pitting degree parameters.

Degree	Number of Pits with Different Sizes			Pitting Area	Pit Depth ti (mm)
	dp1 = 0.2 mm	dp2 = 0.3 mm	dp3 = 0.4 mm
1	20			6.3%	0.1
2	84	20		27.8%	0.15
3	204	84	20	41.7%	0.2

.

$$n_i=\left\{\begin{array}{l}\frac{HL{P}_i}{\pi d_{\mathrm{p}1}^2/4},i=1;\\ \frac{HL{P}_i}{\pi d_{\mathrm{p}1}^2/4}-{\displaystyle \sum _{j=1}^i{n}_{i-j+1}{(\frac{d_{\mathrm{p} j}}{d_{\mathrm{p}1}})}^2,i\ge 2.}\end{array}\right.$$

(28)

where P_i is the percentage of the pitting area; i = 1, 2, 3, where 1 represents slight pitting, 2 represents moderate pitting and 3 represents severe pitting; d_pi is the pitting diameter of the i-th pitting degree; and H is the height of the pitting area.

2.2.3. Calculation of Time-Varying Meshing Stiffness of Pitting Gear

From the analysis of the potential energy method, the contact line length during gear meshing greatly influences the calculation of the potential energy and the time-varying meshing stiffness of the gear. The contact length during gear meshing is equal to the tooth width L for healthy gears. However, for pitted gears, the effective contact length is decreased. Given the decreased tooth contact length represented by ΔL, the effective contact length would be L−L. After modeling pits in Section 2.2.2, the decreased tooth contact length ΔL can be calculated as follows.

In previous studies, the segmentation method usually calculates the meshing stiffness of gears with tooth surface pitting. This paper uses the image processing method to obtain the mesh line length in combination with [29]. It is assumed that the tooth contact line consists of N pixels, as shown by the golden line in Figure 6. For each pixel, a special color or gray indicates whether the point is normal or pitted. More specifically, if the point is in the normal region, it is set to white with a grayscale of 255; if the point is in the pitting area, it is set to black, and the grayscale is 0. In this way, all the pixels on the tooth contact line are distinguished by their colors. Finally, if the black number is ΔN, the reduced tooth contact length is ΔL = ΔNL/N.

Therefore, the Hertzian contact stiffness, bending stiffness, shear stiffness and axial compression stiffness of pitted teeth were calculated as follows [29]:

$$k_{\mathrm{h}}=\frac{E(L-\Delta L_1)\pi }{4(1-{\mu }^2)}$$

(29)

$$\begin{array}{l}\frac{1}{k_{\mathrm{b}}}={\displaystyle {\int }_0^d\frac{{[\cos {\alpha }_1(d-x)-\sin {\alpha }_{1}h]}^2}{E{I}_x}\mathrm{d}x-}{\displaystyle {\int }_0^d\frac{{[\cos {\alpha }_1(d-x)-\sin {\alpha }_{1}h]}^2}{E{I}_x}\frac{\Delta L}{L}\mathrm{d}x+{\displaystyle {\int }_0^d\frac{{[\cos {\alpha }_1(d-x)-\sin {\alpha }_{1}h]}^2}{E{I}^{\prime }{}_x}\frac{\Delta L}{L}\mathrm{d}x}}\\ +{\displaystyle {\int }_0^{r_b-r_f}\frac{{[\cos {\alpha }_1(d+x_1)-\sin {\alpha }_{1}h]}^2}{E{I}_{x1}}\mathrm{d}{x}_1}\end{array}$$

(30)

$$\frac{1}{k_{\mathrm{a}}}={\displaystyle {\int }_0^d\frac{{\sin }^2{\alpha }_1}{E{A}_x}\mathrm{d}x-}{\displaystyle {\int }_0^d\frac{{\sin }^2{\alpha }_1}{E{A}_x}\frac{\Delta L}{L}\mathrm{d}x+{\displaystyle {\int }_0^d\frac{{\sin }^2{\alpha }_1}{E{A}_x^{\prime }}\frac{\Delta L}{L}\mathrm{d}x+{\displaystyle {\int }_0^{r_b-r_f}\frac{{\sin }^2{\alpha }_1}{E{A}_{x1}}\mathrm{d}{x}_1}}}$$

(31)

$$\frac{1}{k_{\mathrm{s}}}={\displaystyle {\int }_0^d\frac{1.2{\cos }^2{\alpha }_1}{G{A}_x}\mathrm{d}x-}{\displaystyle {\int }_0^d\frac{1.2{\cos }^2{\alpha }_1}{G{A}_x}\frac{\Delta L}{L}\mathrm{d}x+{\displaystyle {\int }_0^d\frac{1.2{\cos }^2{\alpha }_1}{G{A}_x^{\prime }}\frac{\Delta L}{L}\mathrm{d}x+{\displaystyle {\int }_0^{r_b-r_f}\frac{1.2{\cos }^2{\alpha }_1}{G{A}_{x1}}\mathrm{d}{x}_1}}}$$

(32)

where ΔL₁ is the contact length at the meshing position, as shown by the red solid line in Figure 6, ΔL is the reduction in the contact length from the base circle x, as shown by the yellow line in Figure 6, and I_x and A_x are the effective moment of inertia and the effective cross-sectional area of the distance x from the base circle, respectively, when pitting corrosion occurs.

$$I_x^{\prime }=\frac{1}{12}{\left(2h_x-t_i\right)}^{3}L$$

(33)

$$A_x^{\prime }=\left(2h_x-t_i\right)L$$

(34)

where t_i is the depth of pitting. By substituting Equations (33) and (34) into Equations (29)–(32), the following formulae can be obtained [29]:

$$\begin{array}{l}\frac{1}{k_{\mathrm{b}}}={{\displaystyle \int }}_{-{\alpha }_1}^{{\alpha }_2}\frac{3\left({\alpha }_2-\alpha \right)\cos \alpha {\{1+\cos {\alpha }_1[({\alpha }_2-\alpha )\sin \alpha -\cos \alpha ]\}}^2}{2EL{[\sin \alpha +({\alpha }_2-\alpha )\cos \alpha ]}^3}\mathrm{d}\alpha -{{\displaystyle \int }}_{-{\alpha }_1}^{{\alpha }_2}\frac{3\Delta L\left({\alpha }_2-\alpha \right)\cos \alpha {\{1+\cos {\alpha }_1[({\alpha }_2-\alpha )\sin \alpha -\cos \alpha ]\}}^2}{2E{L}^2{[\sin \alpha +({\alpha }_2-\alpha )\cos \alpha ]}^3}\mathrm{d}\alpha \\ +{{\displaystyle \int }}_{-{\alpha }_1}^{{\alpha }_2}\frac{3\Delta L\left({\alpha }_2-\alpha \right)\cos \alpha {\{1+\cos {\alpha }_1[({\alpha }_2-\alpha )\sin \alpha -\cos \alpha ]\}}^2}{2E{L}^2{[\sin \alpha +({\alpha }_2-\alpha )\cos \alpha -\frac{t_i}{2R_b}]}^3}\mathrm{d}\alpha +{{\displaystyle \int }}_0^{r_{\mathrm{b}}-r_f}\frac{3{[R_b-R_b\cos {\alpha }_1\cos {\alpha }_2+x_1\cos {\alpha }_1]}^2}{2EL{h}_{x_1}^3}\mathrm{d}{x}_1\end{array}$$

(35)

$$\begin{array}{l}\frac{1}{{\mathrm{k}}_{\mathrm{s}}}={\displaystyle {\int }_{-{\alpha }_1}^{{\alpha }_2}\frac{1.2(1+\nu )({\alpha }_2-\alpha )\cos \alpha {\cos }^2{\alpha }_1}{EL[\sin \alpha +({\alpha }_2-\alpha )\cos \alpha ]}}\mathrm{d}\alpha -{\displaystyle {\int }_{-{\alpha }_1}^{\alpha 2}\frac{1.2\Delta L(1+\nu )({\alpha }_2-\alpha )\cos \alpha {\cos }^2{\alpha }_1}{E{L}^2[\sin \alpha +({\alpha }_2-\alpha )\cos \alpha ]}}\mathrm{d}\alpha \\ +{\displaystyle {\int }_{-{\alpha }_1}^{{\alpha }_2}\frac{1.2\Delta L(1+\nu )({\alpha }_2-\alpha )\cos \alpha {\cos }^2{\alpha }_1}{E{L}^2[\sin \alpha +({\alpha }_2-\alpha )\cos \alpha -\frac{t_i}{2R_b}]}}\mathrm{d}\alpha +{\displaystyle {\int }_0^{r_{\mathrm{b}}-r_f}\frac{1.2(1+\nu ){\cos }^2{\alpha }_1}{EL{h}_{x1}}\mathrm{d}{x}_1}\end{array}$$

(36)

$$\begin{array}{l}\frac{1}{{\mathrm{k}}_{\mathrm{a}}}={\displaystyle {\int }_{-{\alpha }_1}^{\alpha 2}\frac{\left({\alpha }_2-\alpha \right)\cos \alpha {\sin }^2{\alpha }_1}{2EL[\sin \alpha +({\alpha }_2-\alpha )\cos \alpha ]}}\mathrm{d}\alpha -{\displaystyle {\int }_{-{\alpha }_1}^{\alpha 2}\frac{\Delta L\left({\alpha }_2-\alpha \right)\cos \alpha {\sin }^2{\alpha }_1}{2E{L}^2[\sin \alpha +({\alpha }_2-\alpha )\cos \alpha ]}}\mathrm{d}\alpha \\ +{\displaystyle {\int }_{-{\alpha }_1}^{\alpha 2}\frac{\Delta L\left({\alpha }_2-\alpha \right)\cos \alpha {\sin }^2{\alpha }_1}{2E{L}^2[\sin \alpha +({\alpha }_2-\alpha )\cos \alpha -\frac{t_i}{2R_b}]}}\mathrm{d}\alpha +{\displaystyle {\int }_0^{r_b-r_f}\frac{{\sin }^2{\alpha }_1}{2EL{h}_{x1}}\mathrm{d}{x}_1}\end{array}$$

(37)

The time-varying meshing stiffness of gears with different degrees of pitting can be evaluated using the stiffness calculation formula mentioned above. The calculation results of the time-varying meshing stiffness of gears with pitting are shown in Figure 8.

From Figure 8, it can be observed that as the degree of pitting increases, the meshing stiffness gradually decreases. As the pitting spreads further, there is also a significant decrease in the stiffness of the initial engagement of the gears. At the same time, it can be observed that the influence of pitting on stiffness is different from that of cracks. The stiffness of a gear pair with cracks shows a smooth curve, whereas pitting shows more pronounced irregular fluctuations. This is because pitting affects the effective contact length of the gear teeth, which in turn affects the stiffness calculation. Some scholars have found that pitting primarily affects the time-varying meshing stiffness of gears by influencing the Hertzian contact stiffness [30].

2.2.4. Stiffness Calculation of Pitting Corrosion-Crack Composite Failure

In this study, it is assumed that the driving gear has pitting. In contrast, the driven gear has cracks, resulting in two types of gear fault situations during the meshing process, as shown in Figure 9. When two faulty gears mesh, their comprehensive meshing stiffness is calculated using Equation (16). The calculated results of the composite fault stiffness are shown in Figure 9, Figure 10, Figure 11 and Figure 12. In this paper, the crack degree is expressed by its propagation depth, q = 0.5 mm, 0.75 mm, 1 mm, 1.25 mm, 1.5 mm, 2 mm. The degree of pitting is expressed as slight (SL), moderate (M) and severe (Se).

Figure 10a–c shows the variation of the time-varying mesh stiffness with the variation of the gear crack depth, while keeping the pitting severity constant. As can be observed from the graph, with an increase in crack depth, the time-varying mesh stiffness gradually decreases, which is evident from the significant drop in the first double-tooth region during the faulty gear meshing and a gradual decrease in the drop magnitude in the second double-tooth region. In comparison, the time-varying mesh curve of Figure 5, which describes only the crack fault, is smooth. However, for the composite fault stiffness described in Figure 10, the time-varying mesh stiffness shows irregular fluctuations. As the pitting severity increases, this type of fluctuation becomes more pronounced.

Figure 11a–c describes the mesh stiffness of gear coupling cracks of three depths (1 mm, 1.5 mm and 2 mm) with slight pitting. The coupling mesh stiffness is lower than the single-fault mesh stiffness and approaches the crack stiffness. After the crack depth exceeds 1.5 mm, the stiffness decreases significantly. From these figures, it can be concluded that during the slight pitting stage the coupling fault behaves similarly to a single crack fault. This is because during the slight pitting stage the material peeling on the tooth surface is relatively small, and its effect on the comprehensive mesh stiffness is smaller than that of the crack. However, the features of pitting stiffness can still be observed in the coupling stiffness, which exhibits irregular fluctuations.

The number of pits at the severe pitting stage is twice that of the moderate pitting stage, and the pitting craters have spread over almost the entire tooth surface, as shown in Figure 6d. Due to the further expansion of the failed area caused by pitting, the stiffness reduction in the double-tooth and single-tooth areas is extremely significant compared to slight and moderate pitting, as shown in Figure 13. The effect of pitting on stiffness exceeds the effect of cracking at all stages. The coupling stiffness is close to the pitting stiffness.

In addition, the accuracy of the stiffness calculation method mentioned above has been demonstrated in [26,29,30]. The details of the finite element method are as follows. Finite element software is used for the finite element calculation in the paper. The finite element model uses the element type SOLID 187. The teeth are mapped with tetrahedral elements. In the finite element model, the linear material is assumed to be linear elastic. Surface to surface contact (Conta174 and Targe170) is adopted in the finite element modeling process. Gear pair solution setting and stiffness calculation are per [35]. The finite element mesh model is shown in Figure 14.

The comparison between the stiffness obtained from the proposed method and the results from finite element simulations is shown in Figure 15. The comparison shows that both the proposed method and the finite element method produce consistent results in the double-tooth meshing area for both normal and faulty gears. However, there are some differences between the two methods in the single-tooth meshing region. This is because, in the actual meshing process, there is a phenomenon of tooth engagement advance and tooth disengagement lag due to the deformation of the gear teeth, resulting in a decrease in the carrying time of the single-tooth area. Nevertheless, the results indicate that the proposed method and the finite element method are consistent in terms of numerical values and trends, demonstrating the proposed method’s effectiveness in calculating the stiffness.

3. Simulation and Experimental

3.1. Modeling of Two-Stage Gear Transmission System

The gearbox’s physical picture and schematic diagram is shown in Figure 16. The gear transmission system is divided into three parts: gear, bearing and shaft, where 1a, 1b, 2a, 2b, 3a and 3b are bearings; p1 and g1 form the first pair of spur gears; and p2 and g2 form the second pair of spur gears. Gear, shaft and bearing parameters are given in

Table 3. Gear parameters of system.

Parameters	P1	g1	P2	g2
Number of teeth	29	95	36	90
Pressure angle (°)	20	20	20	20
Face width (mm)	12	12	12	12
Modulus (mm)	1.5	1.5	1.5	1.5
Rotary inertia (kg·m2)	0.000205	0.001810	0.000380	0.003492
Mass (kg)	0.1085	1.4648	0.1836	1.3114
The angle of inclination of the tooth line (°)	0	0	0	0
The coefficient of addendum	1	1	1	1
The coefficient of bottom clearance	0.25	0.25	0.25	0.25
The modification coefficient	0	0	0	0

,

Table 4. Shaft parameters of system.

Parameters	Shaft 1	Shaft 1	Shaft 1
Length (mm)	240	160	180
Radius (mm)	20	20	20
Shear modulus (Pa)	8 × 1010	8 × 1010	8 × 1010
Elastic modulus (Pa)	2.1 × 1011	2.1 × 1011	2.1 × 1011
Density (kg/m3)	7850	7850	7850

and

Table 5. Bearing parameters of system.

Parameters	Value
Number of rolling elements	8
Pitch diameter (mm)	37.65
Diameter of inner raceway (mm)	28.7
Diameter of outer raceway (mm)	46.6
Roller diameter (mm)	8.7
Radial clearance (mm)	0.5
Curvature radius (mm)	4.5
Inner raceway curvature sum (1/mm)	0.3078
Outer raceway curvature sum (1/mm)	0.1956
Goodness of fit	0.5172

. The gear transmission system is an elastic system with infinite degrees of freedom, but the finite element method provides us with the possibility of solving it [36,37]. In this paper, according to the structure of the transmission system, the finite element method is used to divide it into three parts: shaft segment unit, gear unit and bearing unit. The finite element model is shown in Figure 17.

3.2. Axial Segment Element Modeling

Considering the flexibility of the shaft, this paper adopts Timoshenko beam elements to establish shaft segment units. As the model in this paper is a parallel-axis spur gear transmission system, the transmission shaft mainly bears the effects of torsion and bending, with a relatively small force acting along the axial direction. Therefore, the axial excitation is ignored. The first node is connected to the input point and is constrained by the torsion degree of freedom. Thus, only two translation degrees of freedom exist at the first node. The other nodes are considered to have three degrees of freedom: two translational degrees of freedom, x and y, and one rotational degree of freedom θ. The stiffness matrix and the mass matrix of the shaft components are calculated as follows:

$$k_e=\left[\begin{array}{cccccc}\frac{GA}{Kl}& 0& 0& -\frac{GA}{Kl}& 0& 0\\ 0& \frac{GA}{Kl}& 0& 0& -\frac{GA}{Kl}& 0\\ 0& 0& \frac{GA}{l}& 0& 0& -\frac{GA}{l}\\ -\frac{GA}{Kl}& 0& 0& \frac{GA}{Kl}& 0& 0\\ 0& -\frac{GA}{Kl}& 0& 0& \frac{GA}{Kl}& 0\\ 0& 0& -\frac{GA}{l}& 0& 0& \frac{GA}{l}\end{array}\right]$$

(38)

$$M_e=\frac{\pi r^{3}l}{6}\left[\begin{array}{cccccc}2& 0& 0& 1& 0& 0\\ 0& 2& 0& 0& 1& 0\\ 0& 0& \frac{2J}{A}& 0& 0& \frac{J}{A}\\ 1& 0& 0& 2& 0& 0\\ 0& 1& 0& 0& 2& 0\\ 0& 0& \frac{J}{A}& 0& 0& \frac{2J}{A}\end{array}\right]$$

(39)

In the equation, l and A represent the length and cross-sectional area of the shaft element, respectively. G, J and K represent the shear modulus, polar moment of inertia and cross-sectional shape factor, respectively. The shaft segment elements are shown in Figure 18.

Assuming that the displacement column vector of the shaft element nodes in the local coordinate system for the i-th beam element is $q_e=[x_i,y_i,{\theta }_{zi},x_{i+1},y_{i+1},{\theta }_{zi+1}]$.

The motion differential equation for the shaft element is expressed as:

$$M_e{\stackrel{..}{q}}_e+C_e{\stackrel{.}{q}}_e+K_e{q}_e=0$$

(40)

where $M_e$ is the consistent mass matrix of the i-th shaft element, $K_e$ is the stiffness matrix of the shaft element and $C_e$ is the damping matrix of the shaft element, calculated using the Rayleigh damping method, with the calculation formula as follows:

$$C_e=a_0{M}_e+a_1{K}_e$$

(41)

where a₀ and a₁ are the proportional coefficients of the Rayleigh damping for the mass matrix and stiffness matrix, respectively.

3.3. Modeling of Gear Meshing Unit

$q_s=[x_{\mathrm{p}},y_{\mathrm{p}},{\theta }_{z\mathrm{p}},x_{\mathrm{g}},y_{\mathrm{g}},{\theta }_{z\mathrm{g}}]$ is selected as the node displacement column vector of gear meshing unit, and the vibration displacement of each gear is projected to the direction of the meshing line. The relative total deformation of gear meshing unit along the direction of the meshing line can be expressed as:

$$\delta =V{q}_s-e_s$$

(42)

where e_s is the integrated error of gear meshing, $V$ is the projection vector of the upward displacement of each gear pair along the meshing line, which can be expressed in the following form.

$$V=[\sin \alpha ,\cos \alpha ,-r_{\mathrm{p}},-\sin \alpha ,-\cos \alpha ,-r_{\mathrm{g}}]$$

(43)

where r_p and r_g are the base circle radius of the driving and driven wheels, respectively, and α is the pressure angle. The meshing unit is shown in Figure 19.

According to Newton’s second law, the differential equation of motion of the spur gear meshing element can be expressed as:

$$\left\{\begin{array}{l}{m}_{\mathrm{p}}\stackrel{..}{x_{\mathrm{p}}}+c_{\mathrm{s}}\stackrel{.}{\delta }\sin \alpha +k_{\mathrm{s}}(t)\delta \sin \alpha =0\\ m_{\mathrm{p}}\stackrel{..}{y_{\mathrm{p}}}+c_{\mathrm{s}}\stackrel{.}{\delta }\cos \alpha +k_{\mathrm{s}}(t)\delta \cos \alpha =0\\ I_{\mathrm{p}}\stackrel{..}{{\theta }_{\mathrm{zp}}}-c_{\mathrm{s}}\stackrel{.}{\delta }{r}_{\mathrm{p}}-k_{\mathrm{s}}(t)\delta r_{\mathrm{p}}=0\\ m_{\mathrm{g}}\stackrel{..}{x_{\mathrm{g}}}-c_{\mathrm{s}}\stackrel{.}{\delta }\sin \alpha -k_{\mathrm{s}}(t)\delta \sin \alpha =0\\ m_{\mathrm{g}}\stackrel{..}{y_{\mathrm{g}}}-c_{\mathrm{s}}\stackrel{.}{\delta }\cos \alpha -k_{\mathrm{s}}(t)\delta \cos \alpha =0\\ I_{\mathrm{g}}\stackrel{..}{{\theta }_{\mathrm{zg}}}-c_{\mathrm{s}}\stackrel{.}{\delta }{r}_{\mathrm{g}}-k_{\mathrm{s}}(t)\delta r_{\mathrm{g}}=0\end{array}\right.$$

(44)

where m_p, m_g is the mass of the driving and driven wheels, respectively, I_zp, I_zg are, respectively, the moment of inertia of the driving and driven wheel around the z axis, C_s is the damping of the meshing element and k_s(t) is the time-varying meshing stiffness of gears.

The motion differential equation matrix of the gear meshing element can be expressed in the following form:

$$M_{\mathrm{s}}{\stackrel{..}{q}}_{\mathrm{s}}+C_{\mathrm{s}}({\stackrel{.}{q}}_{\mathrm{s}}-\stackrel{.}{e})+K_e(q_e-e)=0$$

(45)

where M_s is the mass matrix of the meshing element, The specific form is M_s = diag[m_p, m_p, I_zp, m_g, m_g, I_zg], K_s is the stiffness matrix of the gear mesh element that can be written as K_s = k_s(t)V^TV, C_s is the damping matrix of the gear mesh element that can be written as C_s = c_sV^TV and e is the mesh error vector of the gear mesh element.

3.4. Modeling of Bearing Unit

The bearing supports the shaft system and transmits vibration from the gears to the housing. The bearing element is modeled using springs and dampers as shown in Figure 20. The time-varying bearing stiffness is expressed in Equation (46):

$$k_j(t)=k_{\mathrm{a}}+k_0\sin (2\pi f_{b}t+{\beta }_0)$$

(46)

where j is the radial direction of the bearing, j = x, y, k_a = 8.5 × 10⁸ N/m is the static stiffness of the bearing, k₀ is the fluctuation amplitude of bearing stiffness, f_b is the bearing passing frequency and β₀ is the bearing phase angle.

3.5. Overall Dynamics Model of Two-Stage Gear Transmission System

The differential equation of the overall motion of the system is obtained by integrating each unit equation, which is expressed as follows:

$$M\stackrel{..}{X}(t)+C\stackrel{.}{X}(t)+KX(t)=F$$

(47)

where X(t) is the overall node displacement column vector, M, C, K is the total mass matrix, total damping matrix and total stiffness matrix of the system and F is the external load column vector. The total stiffness matrix of the system is shown in Figure 21.

3.6. Simulation Results

Section 2 obtains the time-varying meshing stiffness of normal gear, cracked gear, pitted gear and gear coupling pitted gear and cracked gear. The Newmark-β method is used to solve the vibration signal results of the two-stage gear system under the conditions of health and failure. The Newmark direct integration method is a generalization of the integral form of linear acceleration method and a simplified algorithm for linear systems with multiple degrees of freedom. When the selected control parameters satisfy a certain relation, the method is unconditionally stable, and the time step size does not affect the stability of the solution. Therefore, this method is used to solve the multi-degree-of-freedom system dynamics in this paper. The basic parameters of the simulation model are as follows: input shaft rotation frequency f_input = 28 Hz, intermediate shaft rotation frequency f_middle = 8.547 Hz and output shaft rotation frequency f_output = 3.419 Hz. The meshing frequencies for the two-stage gear drive are 812 Hz (f_m1) and 307.705 Hz (f_m2), respectively. The bearing ball passing frequencies are 86 Hz (f_b1), 26.2891 Hz (f_b2) and 10.516 Hz (f_b3), respectively. Two more severe failures are selected as components of the subsequent simulation system to observe the phenomena caused by the failures. The two defects are severe pitting and 2 mm cracks, respectively. Pitting failure frequency is f_pit = f_middle = 8.547 Hz. The period is T_pit = 0.117 s, the crack failure frequency is f_crack = _foutput = 3.419 Hz. The period is T_crack = 0.292 s. The load moment is 16.7 N.

This paper further analyzes the gear fault using the vibration acceleration signal. Figure 22 shows the acceleration signals of a normal gear system, a gear system with cracks, a gear system with pitting and a gear system with both pitting and cracks. Obviously, the time-domain signal of the normal gear system is stable without any abnormality, while the vibration acceleration signal of the system containing only cracks or pitting generates a large periodic impact pulse. The impact period of the cracked gear system is 0.292 s, which is the reciprocal of the rotation frequency of the shaft where the cracked gear is located in the gear pair. The impact period of the gear system with pitting is 0.117 s, which is the reciprocal of the rotation frequency of the shaft where the pitting gear is located in the gear pair. The vibration acceleration signals of the gear system containing cracks and pitting simultaneously generate three kinds of periodic impact pulses; one is the meshing period of the cracked gear meshing with normal gear (0.292 s), the second is the meshing period of the pitted gear meshing with the normal gear (0.117 s) and the third is the meshing period of the cracked gear meshing with the pitted gear (0.584 s). The reason for the periodic effect is that the gear is pitted or cracked, so when the gear pair meshes with the faulty gear, the time-varying meshing stiffness of the gear pair decreases, which is manifested as an impact in the vibration response. Suppose both the driving and the driven wheels have faults during meshing. In that case, the time-varying meshing stiffness of the gear pair will drop sharply, manifesting itself as a shock in the vibration response, and the impact amplitude will be greater than that of the single fault. When the gear meshes with the faulty gear, the normal gear participates in the meshing, the time-varying meshing stiffness returns to the normal value, and the impact disappears. When the gear continues to rotate, the impact signal representing the fault will appear periodically as the rotation progresses. Compared with the whole amplitudes of the normal and the faulty gearboxes, pitting and cracking faults are local faults that cause only a local increase in the vibration signal amplitudes and do not have an overall effect on the increase or decrease in the vibration signal.

The time domain curve is transformed into a frequency domain curve by Fast Fourier Transform to observe more detailed fault characteristics. Figure 23 shows the frequency domain signals of a normal gear system, a gear system with cracks, a gear system with pitting and a gear system with both pitting and cracks. When a gear fault occurs, many sidebands appear in the spectrum diagram in addition to the mesh frequency. The sidebands occur mainly around the mesh frequency, and the interval between the sidebands is the rotation frequency of the axis where the faulty gear is located. Therefore, the sideband interval of the pitting spectrum signal is f_pit = 8.547 Hz. The sideband interval of the crack fault spectrum signal side corresponds to f_crack = 3.419 Hz. The system with simultaneous pitting and cracking shows a more complex sideband, and the sideband signal contains f_pit = 8.547 Hz and f_crack = 3.419 Hz. Therefore, the fault location in the transmission system can be found by analyzing the sideband components, which provides a theoretical basis for fault diagnosis.

3.7. Experimental Results

To verify the proposed composite failure model, a two-stage gear drive system test rig has been set up. The two-stage gear drive test system consists of a drive motor, an acceleration sensor, a controller, a gearbox, a magnetic powder brake and a data acquisition system (DT9837, frequency: 8000 Hz). Gear accuracy grade is 8. The cracked gear and pitted gears are shown in Figure 24. The vibration acceleration signal of the gear system can be collected from the accelerometer located at the bearing end plate, and the sampling frequency is 10.24 kHz. The experimental input rotation frequency is 28 Hz. Due to slip in the transmission process, the actual input is about 27.345 Hz, the intermediate shaft rotation frequency is 8.347 Hz, the output shaft rotation frequency is 3.339 Hz, first stage drive frequency is 793.005 Hz, the second stage drive frequency is 300.492 Hz and the load is 4 V (about 16.7 N).

Experimental vibration acceleration results are shown in Figure 25. It can be found from the time domain of the normal gear system that the state is stable and no abnormal state appears in the running process, while the time domain response of the gear system with cracks and pitting has obvious periodic impact, and the impact period is consistent with the simulation results. The spectrum diagram of the experimental signal is shown in Figure 26. The sideband component in the frequency domain of the experimental signal is relatively complex; even in the normal state, the sideband component will appear due to noise. Although there is noise in the fault state, the sideband component generated by cracks and pitting faults is also prominent, and the interval between the edge frequencies is also represented as the fault frequency, which is consistent with the simulation results. However, due to the attenuation of vibration energy in the process of experimental measurement and the influence of environmental noise, the amplitude of the simulated signal is different from that of the experimental signal.

The accelerometers are mounted on the housings in the experiment, and the simulation result is the bearing acceleration. Energy decay is inevitable in the transmission path of the vibration signal. The interface between gear and shaft, inner race and outer race and outer race and housing will cause significant energy loss. The simulation does not consider the mass eccentricity and assembly error of the faulty gear. Therefore, the amplitude of the experimental results is lower than that of the simulated results and the excitation of the test is more complex.

4. Conclusions

In this paper, a stiffness calculation model of pitted and cracked gear composite faults considering structural coupling was proposed based on the energy method. A dynamic model of two-stage gear transmission system with pitted and cracked faults was established. A two-stage gear box test platform verified the accuracy of the proposed model. The main conclusions drawn are as follows:

(1)

The presence of cracks and pitting reduces the meshing stiffness of the gear, but the stiffness curve after cracks is smooth, while the stiffness curve after pitting has irregular fluctuations. The composite failure stiffness will approach the failure stiffness of a fault type with greater fault degree.

(2)

In the time domain, pitting and cracking as local faults will produce periodic vibration and impact with the operation of gears. Coupling faults will produce three kinds of vibration shocks; one is pitting impact, one is cracking impact and one is coupling faults impact. The amplitude of coupling faults is obviously larger than that of single fault impact.

(3)

In the frequency domain, the occurrence of pitting and cracking faults will lead to the occurrence of sidebands near each order of harmonics in the spectrum, and the interval between sidebands is mainly the rotation frequency of the shaft where the faulty gear is located. The coupling faults will result in a more complex sideband. The sideband composed of two fault frequencies will appear near each harmonic.

(4)

Although the experimental signal is in general agreement with the simulation signal, the gap between the simulation signal and the experimental signal is inevitable due to a series of factors such as the extraction position of the experimental speed signal, attenuation of vibration energy in the transmission process, friction in the transmission process and various errors.

Vibration analysis of a faulty gear system is an important source of information for fault diagnosis, and it is of great significance to establish a dynamic model of the fault transmission system. The work in this paper will provide some theoretical support for follow-up research coupling fault dynamic modeling and coupling fault diagnosis. In future work, we will pay more attention to the effects of different types of coupling faults on the meshing and dynamic characteristics of multi-stage gear systems.