*Article* **Dynamic Analysis of a High-Contact-Ratio Spur Gear System with Localized Spalling and Experimental Validation**

**Zhenbang Cheng 1,2, Kang Huang 1,3, Yangshou Xiong 4,\* and Meng Sang 1,3**


**Abstract:** The dynamic characteristics and tooth spalling fault features are studied for the highcontact-ratio spur gear bearing system. The bending torsional dynamic model is proposed in this study for the gear bearing system with an ellipsoid spalling fault. This model also considers timevarying meshing stiffness, tooth friction, fractal gear backlash, and comprehensive transmission error. The meshing stiffness of the system is evaluated using the potential energy method. The bifurcation diagram, time-domain waveform, Poincaré map, phase map, frequency spectrum, and related threedimensional map are used as tools to analyze the system's dynamic response qualitatively. The results reveal that the system's motion with ellipsoid tooth spalling defect exhibits rich dynamic behavior. The response of the proposed dynamic model is consistent with experimental results in the frequency domain. Therefore, the developed dynamic model can predict the system's vibration behavior with localized spalling fault. Hence, it could also provide a theoretical foundation for future spall defect diagnosis of the gear transmission system.

**Keywords:** high-contact-ratio gear; tooth spalling; meshing stiffness; nonlinear dynamic behavior

#### **1. Introduction**

The contact ratio is an important index to indicate the smoothness and the uniformity of the gear system's transmission load [1]. By increasing the addendum coefficient, enhancing the number of teeth, reducing the pressure angle, and adjusting the modification coefficient, an involute spur gear with a contact ratio greater than two can be obtained, called high-contact-ratio spur gear [2]. Compared with standard gears, high-contact-ratio spur gear has more pairs of teeth meshing at the same time. As a result, high-contact-ratio gear has the advantages of high load carrying capacity [3]. Thus, it is widely used in automobiles, power tools, and other fields. As a common failure form of gear, tooth spalling is frequently encountered under excessive load, which initially occurs at the local position of the tooth surface. Localized tooth spalling induces serious damage to the gear transmission system gradually. Three main aspects are considered while characterizing tooth spalling fault: meshing stiffness, dynamic characteristics, and experimental validation. The model of damaged gear pairs is analyzed to estimate the mesh stiffness. The mesh stiffness is imported into the gear dynamic model to evaluate the non-linear dynamic characteristics. The dynamic response obtained from the simulation is compared with the experimental results to verify the dynamic model.

Tooth spalling affects the tooth profile and results in a change of meshing stiffness. Several researchers have developed many methods to estimate the meshing stiffness of the gears. Wang [4] and Zhan [5] have developed the finite element methods. The finite

**Citation:** Cheng, Z.; Huang, K.; Xiong, Y.; Sang, M. Dynamic Analysis of a High-Contact-Ratio Spur Gear System with Localized Spalling and Experimental Validation. *Machines* **2022**, *10*, 154. https://doi.org/ 10.3390/machines10020154

Academic Editors: Sven Matthiesen and Thomas Gwosch

Received: 12 January 2022 Accepted: 14 February 2022 Published: 18 February 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

element method can accurately simulate the contact state of the system, but it is timeconsuming. The analytical method is very efficient in terms of computation [6]. For example, Sun [7] studied the mesh stiffness of spur gear pairs with tooth modifications based on the thin slice assumption. For the gear with rectangular tooth spalling, Chaari researched the mesh stiffness of spur-gear pair [8], Jiang [9] and Han [10] studied the helical gear, while Luo studied the planetary gear [11]. In references [12,13], the meshing stiffness of gears with rectangular spalls of various widths, lengths, and positions were compared. Saxena et al. [14] estimated gear mesh stiffness for rectangular, circular, and V-shaped spalling. It is assumed that the bottom of the tooth spalling is flat in the literature. In practice, the shape of the tooth spall usually has a curvilinear bottom with a gradually changing dent depth rather than a suddenly decreased tooth thickness. Compared to rectangular solid, the ellipsoid is close to the actual geometry of the tooth spalling. This geometric model makes the calculation of mesh stiffness more accurate. The mesh stiffness of the gear with ellipsoid tooth spall is analyzed in detail in this work.

The estimation of mesh stiffness lays a foundation for gear dynamics. Dynamic analysis is valuable for operating status monitoring and gear fault diagnosis [15,16]. Using statistical methods to evaluate spall severity was suitable under low velocity and low excitation [17]. Dadon et al. [18] researched the effect of different gear imperfections on fault detection. Ma and Chen [19] studied the differences in vibration signals of tooth crack and tooth spall. Modification coefficients were introduced to research the impact of the spalls on gear dynamics [20]. Chen et al. [21] compared the mesh characteristics of helical gears with spalling faults using analytical and finite-element methods. Based on the theoretical and experimental study, Huangfu et al. [22] investigated the meshing and dynamic characteristics of a spalled gear system. Luo et al. [23] demonstrated that tooth spall and sliding friction have an evident effect on kinetic performance. Shi et al. [24] discussed the dynamic characteristics of the gear with double-teeth spalling fault. In terms of global dynamics, Ma [25] employed the singularity theory to evaluate the bifurcation characteristics of the gear system with tooth spalling.

Two drawbacks in these previous models are found. Mesh stiffness of the high-contactratio gear with ellipsoid spalling defects is not calculated. Secondly, motion state analysis of the gear with localized spalling defects under excitation frequency and gear backlash is not considered. The above deficiencies are the main contributions of this study. Highlights of this paper are listed as follows.


The rest of the article is arranged as follows. Section 2 describes the modeling of the mesh stiffness through the precise tooth profile equation of the high-contact-ratio gear and the ellipsoid equation. Section 3 illustrates the proposed dynamic model of the gearbearing system. Section 4 discusses the numerical results. Section 5 is about the vibration experiments. The conclusions are presented in Section 6.

#### **2. Mesh Stiffness Computation**

The alternating meshing process of two and three gear pairs causes the variation of mesh stiffness. It also plays the role of internal excitation of the gear system. It is of great significance to develop an analytical model to calculate the mesh stiffness.

#### *2.1. Accurate Tooth Profile Equation*

The tooth profile of high contact ratio gear comprises tooth tip arc line, involute tooth profile, and tooth root transition curve. The motion of a rack-type cutting tool is equivalent to the meshing of the rack and pinion. In the process of machining, the cutter's machining pitch line is always tangent to the gear's machining pitch circle [26]. The profile of the racktype cutting tool is shown in Figure 1. As presented in Figure 2, during gear machining, the involute part of the tooth profile is cut directly by the straight part of the cutter (BC), and the fillet part of the cutter (AB) cuts the transition part. The transition part is the isometric curve of the extended involute. This extended involute is depicted by the center of the tool's rounded corner. An *x–y* coordinate system with the center of the gear as the origin is shown. The gear teeth are usually considered as a cantilever beam of a variable cross-section. The coordinates of any contact point: *i*, on the involute part are expressed as follows.

$$\begin{cases} \ x\_i = r\_b[(a\_i + \theta\_b)\sin\alpha\_i + \cos\alpha\_i] \\\ y\_i = r\_b[(a\_i + \theta\_b)\cos\alpha\_i - \sin\alpha\_i] \end{cases} \tag{1}$$

*rb* denotes the radius of the gear base circle. *α<sup>i</sup>* is the working load angle acting on the contact point: *i*. *θ<sup>b</sup>* is half of the tooth base arc angle (*θ<sup>b</sup> =* (*π*/*2 + 2X*·tan*α*0)/*N + invα*0). *x* is the displacement coefficient. *N* is the number of teeth. *inv*(·) represents the involute function of the pressure angle. Coordinates of point: *j*, on the transition curve are defined as follows.

$$\begin{cases} \begin{aligned} x\_{\parallel} &= r \cos \varrho - (a\_1/\sin \gamma + r\_\rho) \sin(\gamma - \varrho) \\ y\_{\parallel} &= r \sin \varrho - (a\_1/\sin \gamma + r\_\rho) \cos(\gamma - \varrho) \\ &\quad (a\_0 \le \gamma \le \pi/2) \end{aligned} \tag{2}$$

*r* and *r<sup>ρ</sup>* are the radii of the gear pitch and the top corner of the tool, respectively. *a* is the distance from the center of the top corner to the center line of the tool (*a =* (*ha\* + c\**)*m* − *rρ, a*<sup>1</sup> *= a* − *xm, ϕ =* (*a*1/tan*γ + b*)/*r, b = πm*/4 *+ ha*\*·*m*·tan*α*<sup>0</sup> *+ rρ*·cos*α*0).

**Figure 1.** Rack-type cutting tool profile.

**Figure 2.** Tooth profile curve diagram.

#### *2.2. Analytical Model of Meshing Stiffness*

Many tooth surface spalls have progressively varying depth and a curved base surface in practice as shown in Figure 3. The ellipsoid tooth spall is formed by removing the intersection of the gear and the ellipsoid. As shown in Figure 4, the ellipsoid tooth spall's maximum width, length, and center depth are *ws*, *ls*, and *hs*, respectively. The starting and ending position is denoted by *xstart and xend*, respectively. *θ<sup>s</sup>* is the angle between the tangent of tooth spall and involute curve. In this work, the position (*xstart*, *θs*) and the severity (*ws, ls, hs*) are fixed.

**Figure 3.** The ellipsoid tooth spall on the pinion.

**Figure 4.** A cross-sectional view of the ellipsoid tooth spall.

Table 1 displays the comparison of key geometry parameters of the ellipsoid tooth spall for healthy and spall gear.

**Table 1.** The ellipsoid tooth spall data.


An ellipsoid's parametric function is defined as follows.

$$\frac{x^2}{a^2} + \frac{y^2}{b^2} + \frac{z^2}{c^2} = 1\tag{3}$$

*a*, *b*, and *c* are the radii of the ellipsoid along *x, y,* and *z* axis. Equations (4) and (5) are derived assuming *a* and *b* as equal.

$$a = b = r\_{s\\_max} = \frac{(0.5w\_s)^2 + h\_s^2}{2h\_s} \tag{4}$$

e

$$\mathcal{L} = \frac{1}{2} \sqrt{\frac{r\_{s\\_max}^2 l\_s^2}{2r\_{s\\_max}h\_s - h\_s^2}}\tag{5}$$

In the *x–y* plane, the coordinates of the ellipsoid center are derived as follows.

$$\begin{cases} \begin{aligned} \mathbf{x\_{05}} &= \mathbf{x\_{start}} + r\_{s\\_\text{max}} \cos\left(\frac{\pi}{2} - \theta\_\text{s} - \arcsin\left(\frac{0.5 w\_\text{s}}{r\_{s\\_\text{max}}}\right)\right) \\\ y\_{os} &= -\tan(\theta\_\text{s})(\mathbf{x\_{os}} - \mathbf{x\_{start}}) + \mathbf{y}(\mathbf{x\_{start}}) + \frac{r\_{s\\_\text{max}} - h\_\text{s}}{\cos\theta\_\text{s}}, \quad (h\_\text{s} < r\_{s\\_\text{max}}) \end{aligned} \end{cases} \tag{6}$$

Figure 4b shows the cross-section of the ellipsoid tooth spall at a distance *x*, The ellipse equation is as follows.

$$\frac{y\_s^2}{b\_x^{\prime 2}} + \frac{z\_s^2}{c\_x^{\prime 2}} = 1\tag{7}$$

$$b\_{\mathbf{x}} = \sqrt{b^2 \left(1 - \frac{\left(\mathbf{x}\_{\text{os}} - \mathbf{x}\_t\right)^2}{a^2}\right)}\tag{8}$$

$$c\_x = \sqrt{c^2 \left(1 - \frac{\left(\mathbf{x}\_{os} - \mathbf{x}\_t\right)^2}{a^2}\right)}\tag{9}$$

The length *lxs* is calculated is expressed as follows.

$$d\_{\rm xs} = \sqrt{c\_x^2 \left(1 - (y\_{\rm os} - h\_x)^2 / b\_x^2\right)}\tag{10}$$

As shown in Figure 4b, *hx* is the half-height of the gear tooth cross-section at a distance *x*. The gear tooth contact length *Le* is given as follows.

$$L\_{\mathfrak{e}} = L - \mathfrak{L}\_{\mathfrak{x}s} \tag{11}$$

*L* denotes the gear tooth width. The maximum depth of the ellipsoid tooth spall at a distance *x* is determined as follows.

$$h\_{\rm xs} = b\_{\rm x} - (y\_{\rm os} - h\_{\rm x}) \tag{12}$$

The corresponding area of the portion of the ellipse is deduced as follows.

$$A\_{xs} = 2\int\_{-b\_x}^{-\left(y\_{\infty} - h\_x\right)} \sqrt{\left(1 - \frac{y\_s}{b\_x^2}\right)} c\_x^2 dy\_s \tag{13}$$

The cross-section area *Ax* of the ellipsoid tooth spall at a distance *x* is given below.

$$A\_{\mathbf{x}} = 2h\_{\mathbf{x}}L - A\_{\mathbf{x}s} \tag{14}$$

*Ax* causes a shift in the cross-neutral section's axis. The corresponding displacement *δxs* between these two central axes is calculated as follows.

$$\delta\_{\rm xs} = \frac{A\_{\rm xs} h\_{\rm oxs}}{A\_{\rm x}} \tag{15}$$

The area moment of inertia of the ellipse segments is expressed as follows.

$$I\_{zs} = \int\_{-b\_{\rm x}}^{-(y\_{\rm as} - h\_{\rm x})} 2y\_s^2 \sqrt{\left(1 - \frac{y\_s^2}{b\_{\rm x}^2}\right)} c\_{\rm x}^2 dy\_s - A\_{\rm xs} (y\_{\rm os} - h\_{\rm oxs})^2, \quad (h\_{\rm x} \le y\_{\rm os}) \tag{16}$$

Under the tooth spall conditions, the area moment of inertia of the gear tooth crosssection with respect to the new axis is modified as follows.

$$I\_z = \frac{2Lh\_x^3}{3} + 2Lh\_x \delta\_{xs}^2 - \left(I\_{zs} + A\_{xs} (h\_{axs} + \delta\_{xs})^2\right) \tag{17}$$

Figure 5 shows the gear profile. Segments *AB* and *BC* represent the transition curve and involute curve, respectively. *P* is the meshing point, and the corresponding pressure angle is *αp*. *F* denotes the meshing force that is decomposed into *Fx* and *Fy*, respectively. Based on the potential energy method [27], comprehensive mesh stiffness *k*(*t*) is deduced.

$$k(t) = \sum\_{i=1}^{3} \frac{1}{\frac{1}{k\_{hi}} + \frac{1}{k\_{hi}^T} + \frac{1}{k\_{si}^T} + \frac{1}{k\_{fi}^T} + \frac{1}{k\_{fi}^T} + \frac{1}{k\_{hi}^T} + \frac{1}{k\_{si}^T} + \frac{1}{k\_{ai}^T} + \frac{1}{k\_{fi}^T}} \tag{18}$$

**Figure 5.** Gear profile.

Subscript *i* (*i =* 1, 2, 3) denotes the pair of the engaging gear tooth. Based on the tooth profile equations, each of the stiffness is expressed as follows.

1 *ka* = ⎧ ⎪⎪⎪⎪⎪⎪⎨ ⎪⎪⎪⎪⎪⎪⎩ *xB xA* sin2 *α<sup>p</sup> EA*<sup>1</sup> *dx*<sup>1</sup> <sup>+</sup> *xp xB* sin2 *α<sup>p</sup> EA*<sup>2</sup> *dx*<sup>2</sup> (*xB* < *<sup>x</sup>* ≤ *xstart*) *xB xA* sin2 *α<sup>p</sup> EA*<sup>1</sup> *dx*<sup>1</sup> <sup>+</sup> *xs xB* sin<sup>2</sup> *α<sup>p</sup> EA*<sup>2</sup> *dx*<sup>2</sup> <sup>+</sup> *xp xs* sin2 *α<sup>p</sup> EA*<sup>3</sup> *dx*<sup>3</sup> (*xstart* < *<sup>x</sup>* ≤ *xend*) *xB xA* sin2 *α<sup>p</sup> EA*<sup>1</sup> *dx*<sup>1</sup> <sup>+</sup> *xs xB* sin<sup>2</sup> *α<sup>p</sup> EA*<sup>2</sup> *dx*<sup>2</sup> <sup>+</sup> *xe xs* sin2 *α<sup>p</sup> EA*<sup>3</sup> *dx*<sup>3</sup> <sup>+</sup> *xp xe* sin<sup>2</sup> *α<sup>p</sup> EA*<sup>4</sup> *dx*<sup>4</sup> (*xend* < *<sup>x</sup>*) (19) 1 *kb* = ⎧ ⎪⎪⎪⎪⎪⎪⎨ ⎪⎪⎪⎪⎪⎪⎩ *xB xA M*<sup>2</sup> 1 *EI*<sup>1</sup> *dx*<sup>1</sup> + *xP xB M*<sup>2</sup> 2 *EI*<sup>2</sup> *dx*<sup>2</sup> (*xB* < *x* ≤ *xstart*) *xB xA M*<sup>2</sup> 1 *EI*<sup>1</sup> *dx*<sup>1</sup> + *xs xB M*<sup>2</sup> 2 *EI*<sup>2</sup> *dx*<sup>2</sup> + *xp xs M*<sup>2</sup> 3 *EI*<sup>3</sup> *dx*<sup>3</sup> (*xstart* < *x* ≤ *xend*) *xB xA M*<sup>2</sup> 1 *EI*<sup>1</sup> *dx*<sup>1</sup> + *xs xB M*<sup>2</sup> 2 *EI*<sup>2</sup> *dx*<sup>2</sup> + *xe xs M*<sup>2</sup> 3 *EI*<sup>3</sup> *dx*<sup>3</sup> + *xp xe M*<sup>2</sup> 4 *EI*<sup>4</sup> *dx*<sup>4</sup> (*xend* < *x*) (20)

$$\frac{1}{k\_{s}} = \begin{cases} \int\_{x\_{A}}^{x\_{B}} \frac{1.2 \cos^{2} a\_{F}}{GA\_{1}} dx\_{1} + \int\_{x\_{B}}^{x\_{F}} \frac{1.2 \cos^{2} a\_{F}}{GA\_{2}} dx\_{2} & (\mathbf{x}\_{B} < \mathbf{x} \le \mathbf{x}\_{\text{start}})\\\\ \int\_{x\_{A}}^{x\_{B}} \frac{1.2 \cos^{2} a\_{F}}{GA\_{1}} dx\_{1} + \int\_{x\_{B}}^{x\_{B}} \frac{1.2 \cos^{2} a\_{F}}{GA\_{2}} dx\_{2} + \int\_{x\_{A}}^{x\_{F}} \frac{1.2 \cos^{2} a\_{F}}{GA\_{3}} dx\_{3} & (\mathbf{x}\_{\text{start}} < \mathbf{x} \le \mathbf{x}\_{\text{end}})\\\\ \int\_{x\_{A}}^{x\_{B}} \frac{2.2 \cos^{2} a\_{F}}{GA\_{1}} dx\_{1} + \int\_{x\_{B}}^{x\_{A}} \frac{1.2 \cos^{2} a\_{F}}{GA\_{2}} dx\_{2} + \int\_{x\_{A}}^{x\_{A}} \frac{1.2 \cos^{2} a\_{F}}{GA\_{3}} dx\_{3} + \int\_{x\_{C}}^{x\_{F}} \frac{1.2 \cos^{2} a\_{F}}{GA\_{4}} dx\_{4} & (x\_{cd1} < x) \end{cases} \tag{21}$$

$$\frac{1}{k\_{f}} = \frac{\cos^{2} a\_{P}}{EL} \left\{ L^{\*} \left(\frac{u\_{f}}{s\_{f}}\right)^{2} + M^{\*} \left(\frac{u\_{f}}{s\_{f}}\right) + P^{\*} \left(1 + Q^{\*} \tan^{2} a\_{P}\right) \right\} \tag{22}$$

$$\frac{1}{k\_{\rm li}} = \frac{\pi EL\_{\rm c}}{4(1 - \nu^2)}\tag{23}$$

*E*, *G*, and *υ* stand for Young's modulus, shear modulus, and Poisson's ratio, respectively. *A*<sup>1</sup> and *I*<sup>1</sup> denote the cross-sectional area and corresponding inertia within the transition curve, respectively. *A*2, *A*<sup>4</sup> and *I*2, *I*<sup>4</sup> are the cross-sectional areas of the healthy part of the involute curve and the related area moments of inertia. For the spall part of the gear profile, *A*<sup>3</sup> and *I*<sup>3</sup> are evaluated by Equations (14) and (17), respectively. *uf*, *sf*, *L\**, *M\**, *P,\** and *Q\** are given in reference [28].

Figure 6 shows the curve of mesh stiffness. The incidence of tooth spall reduces the mesh stiffness. A localized spall fault is considered. Hence, it occurs once per revolution. Thus, compared with the healthy gear, the reduction in mesh stiffness occurs mainly in the double-tooth and the triple-tooth zones.

**Figure 6.** Time-varying mesh stiffness.

#### **3. Dynamic Model of System**

The support of the rolling bearing is assumed to be rigid. The gear-bearing translationtorsion dynamic lumped parameter model is established, as is shown in Figure 7, considering the influence of gear backlash, damping, comprehensive transmission error, friction force, and time-varying meshing stiffness.

**Figure 7.** Gear dynamic model.

#### *3.1. Gear-Bearing System*

*m, I, T,* and *θ* represent the mass, moment of inertia, torque, and torsional angular displacement, respectively. Subscripts *p* and *g* indicate quantities associated with pinion and gear, respectively. *rb*, *kx,* and *ky* denote the base circle radius, vertical radial support stiffness, and horizontal radial support stiffness. The time-varying meshing stiffness, meshing

damping, and backlash of the gear pair are expressed by *km*, *cm*, and *2b*, respectively. *e*(*t*) is the comprehensive static error along the tangent direction of the gear base circle.

As shown in Figure 7, the dynamic model has six degrees of freedom, including two rotational degrees of freedom and four translational degrees of freedom along with the horizontal and vertical directions. The generalized coordinate array is expressed as follows.

$$q = \begin{bmatrix} \theta\_{\mathcal{P}}, \theta\_{\mathcal{S}'}, \mathbf{x}\_{\mathcal{P}'}, \mathbf{x}\_{\mathcal{S}'} y\_{\mathcal{P}'} y\_{\mathcal{S}} \end{bmatrix}^T \tag{24}$$

#### *3.2. Dynamic Meshing Force and Frictional Force*

The relative displacement between pinion and gear along the line of action is expressed as follows.

$$\delta(t) = r\_{bp}\theta\_{\mathcal{P}} - r\_{b\mathcal{g}}\theta\_{\mathcal{S}} + (\mathbf{x}\_{\mathcal{P}} - \mathbf{x}\_{\mathcal{S}})\cos\varphi\_{\mathcal{P}\mathcal{S}} + (\mathcal{y}\_{\mathcal{P}} - \mathcal{y}\_{\mathcal{S}})\sin\varphi\_{\mathcal{P}\mathcal{S}} - c(t) \tag{25}$$

The dynamic meshing force between gears consists of elastic meshing forces caused by time-varying stiffness and viscous meshing forces caused by meshing damping, denoted as follows. .

$$F\_d = c\_m \delta(t) + k\_m f(\delta(t), b(t)) \tag{26}$$

The gap function formula can be calculated as follows.

$$f(\delta(t), b) = \begin{cases} \delta(t) - b(t) & (\delta(t) > b(t)) \\ 0 & (|\delta(t)| \le b(t)) \\ \delta(t) + b(t) & (\delta(t) < -b(t)) \end{cases} \tag{27}$$

*b*(*t*) represents the gear backlash. Based on fractal theory, it is expressed as follows [29].

$$b(t) = b\_0 - \frac{R\_{d1}}{R\_{\rm ac}(D\_1)} \sum\_{k=0}^{+\infty} \lambda^{(D\_1 - 2)k} \sin(\lambda^k t) - \frac{R\_{d2}}{R\_{\rm ac}(D\_2)} \sum\_{k=0}^{+\infty} \lambda^{(D\_2 - 2)k} \sin(\lambda^k t) \tag{28}$$

*b0*, *λ, Ra*1/*Ra*2, and *D*1/*D*<sup>2</sup> denote initial gear backlash, characteristic scale coefficient, actual surface roughness, and fractal dimension, respectively. *Rac*(*D*) is the function to get the corresponding *Ra* with a particular fractal dimension, which can be obtained from the reference [30]. The friction force between tooth surfaces during gear meshing is deduced as follows.

$$F\_f = \eta \mu F\_\mathrm{d} \tag{29}$$

The direction of friction is variable, and the direction coefficient depends on the following formula.

*η* = *sign*(*ω*p*KN*<sup>1</sup> − *ωgKN*2) (30)

Friction force arms (Figure 8) are deduced as follows.

$$\begin{aligned} \text{KN}\_1 &= \sqrt{\left(r\_1 + r\_2\right)^2 - \left(r\_{b1} + r\_{b2}\right)^2} - \sqrt{r\_{d2}^2 - r\_{b2}^2} + r\_{b1}\omega\_1 t\\ \text{KN}\_2 &= \sqrt{r\_{d2}^2 - r\_{b2}^2} - r\_{b1}\omega\_1 t \end{aligned} \tag{31}$$

The friction coefficient is related to relative sliding velocity, tooth surface roughness, contact pressure, lubrication situation, etc. Hence, it is particularly difficult to predict the friction coefficient value. The tooth friction model mainly includes the Coulomb friction model, Buckingham empirical formula, Benedict Kelly model, and the friction model based on Elastohydrodynamic lubrication (EHL) theory. Among these, the friction model proposed by Xu [31] reasonably agrees with the measured data, so it is usually adopted to predict the friction coefficient in this paper. It is based on non-Newtonian, thermal EHL theory, and multiple linear regression analysis. The friction coefficient calculation results are shown in Figure 9.

**Figure 8.** Geometrical relationship of gear pair.

**Figure 9.** Time−varying friction coefficient.

The friction torque on the gear can be expressed as follows.

$$\begin{aligned} \mathcal{M}\_p &= \eta \mu F\_d \text{KN}\_1\\ \mathcal{M}\_\% &= \eta \mu F\_d \text{KN}\_2 \end{aligned} \tag{32}$$

#### *3.3. Differential Equations of Motion*

Based on the time-varying meshing stiffness, transmission error, the fractal backlash, the frictional force, and Newton's law of motion, the differential equation of vibration motion of the gear-bearing coupling system is deduced as follows.

$$\begin{cases} \begin{aligned} \label{10} I\_p \ddot{\theta}\_p + F\_d r\_{bp} - M\_p &= T\_d \\ I\_\xi \ddot{\theta}\_\mathcal{g} - F\_d r\_{b\mathcal{g}} + M\_\mathcal{g} &= -T\_b \\ m\_p \ddot{\mathbf{x}}\_p + c\_{px} \dot{\mathbf{x}}\_p + k\_{px} \mathbf{x}\_p + F\_d \cos \boldsymbol{\varrho}\_{p\mathcal{g}} + F\_f \sin \boldsymbol{\varrho}\_{p\mathcal{g}} &= 0 \\ m\_\mathcal{g} \ddot{\mathbf{x}}\_\mathcal{g} + c\_{\mathcal{g}x} \dot{\mathbf{x}}\_\mathcal{g} + k\_{\mathcal{g}x} \mathbf{x}\_\mathcal{g} - F\_d \cos \boldsymbol{\varrho}\_{p\mathcal{g}} - F\_f \sin \boldsymbol{\varrho}\_{p\mathcal{g}} &= 0 \\ m\_p \ddot{\mathbf{y}}\_p + c\_{p\mathcal{y}} \dot{\mathbf{y}}\_p + k\_{p\mathcal{y}} \mathbf{y}\_p + F\_d \sin \boldsymbol{\varrho}\_{p\mathcal{g}} - F\_f \cos \boldsymbol{\varrho}\_{p\mathcal{g}} &= 0 \\ m\_\mathcal{g} \ddot{\mathbf{y}}\_\mathcal{g} + c\_{\mathcal{g}\mathcal{g}} \dot{\mathbf{y}}\_\mathcal{g} + k\_{\mathcal{g}\mathcal{g}} \mathbf{y}\_\mathcal{g} - F\_d \sin \boldsymbol{\varrho}\_{p\mathcal{g}} + F\_f \cos \boldsymbol{\varrho}\_{p\mathcal{g}} &= 0 \end{aligned} \tag{33}$$

Setting *x*<sup>1</sup> = *θp*, *x*2=*x*<sup>1</sup> , *x*<sup>3</sup> = *θg*, *x*<sup>4</sup> = *x*<sup>3</sup> , *x5* = *xp*, *x*<sup>6</sup> = *x*<sup>5</sup> , *x*<sup>7</sup> = *xg*, *x*<sup>8</sup> = *x*<sup>7</sup> , *x*<sup>9</sup> = *yp*, *x*<sup>10</sup> = *x*<sup>9</sup> , *x*<sup>11</sup> = *yg*, and *x*<sup>12</sup> = *x*<sup>11</sup> , the above equation is transformed into the following form.

$$\begin{cases} \begin{aligned} \mathbf{x}\_1' &= \mathbf{x}\_2 \\ \mathbf{x}\_2' &= \frac{T\_d}{l\_p} + \frac{M\_{PL}}{l\_p} - \frac{F\_d r\_{p\_2}}{l\_p} \\ \mathbf{x}\_3' &= \mathbf{x}\_4 \\ \mathbf{x}\_4' &= \frac{-T\_b}{l\_g} - \frac{M\_{Pl}}{l\_g} + \frac{F\_d r\_{b1}}{l\_g} \\ \mathbf{x}\_5' &= \mathbf{x}\_6 \\ \mathbf{x}\_6' &= -\frac{-F\_d \cos\varphi\_{PL}}{m\_p} - \frac{F\_f \sin\varphi\_{PL}}{m\_p} - \frac{\varepsilon\_{pr}\dot{x}\_p}{m\_p} - \frac{k\_{p\times p}}{m\_p} \\ \mathbf{x}\_7' &= \mathbf{x}\_8 \\ \mathbf{x}\_8' &= \frac{F\_d \cos\varphi\_{PL}}{m\_p} + \frac{F\_f \sin\varphi\_{PL}}{m\_g} - \frac{\varepsilon\_{rr}\dot{x}\_r}{m\_g} - \frac{k\_{r\times r}\mathbf{x}\_r}{m\_g} \\ \mathbf{x}\_9' &= \mathbf{x}\_{10} \\ \mathbf{x}\_{10'} &= \frac{-F\_d \sin\varphi\_{PL}}{m\_p} + \frac{F\_f \cos\varphi\_{PL}}{m\_p} - \frac{\varepsilon\_{rr}\dot{y}\_p}{m\_p} - \frac{k\_{p\times p}}{m\_p} \\ \mathbf{x}\_{11'} &= \mathbf{x}\_{12} \end{aligned} \tag{34}$$

#### **4. Numerical Simulation and Discussion**

The high-contact-ratio gear system is non-linear. Table 2 lists its main parameters. The parameters in blue boxes are gear basic parameters. Gear basic parameters were designed via the KISSsoft software. Mesh damping ratio was obtained from the reference [32]. Additionally, the other parameters were calculated from the gear system via the Solidworks software. The fourth-order Runge–Kutta numerical integration method is used to solve the above differential equations of the system through the MATLAB software. The simulated data are processed for generating a bifurcation diagram, three-dimensional frequency spectrum, Poincaré map, phase map, etc.



#### *4.1. Effect of Excitation Frequency*

Gear excitation frequency often changes with working conditions. It is one of the key parameters affecting the dynamic characteristics of the gear system and is often used as a variable parameter to compare the systematic dynamical performance. Here, surface roughness *Ra* is taken equal to 0.8 μm, and fractal dimensions *D*<sup>1</sup> and *D*<sup>2</sup> are equivalent to 1.1. Figure 10 shows the bifurcation characteristics of the lateral displacement (*xp*) of the pinion changing with the excitation frequency Ω. Its three-dimensional frequency map is illustrated in Figure 11. Both figures reach an agreement on systematic bifurcation behaviors.

**Figure 10.** Bifurcation diagram under different excitation frequency.

**Figure 11.** Three−dimensional frequency spectrum under different excitation.

The system is in periodic motion at low excitation frequency, as shown in Figure 12a,b. When Ω equals 0.6218, resonance occurs in the system. The system enters chaotic regions A1, A2, and A3 sequentially with the increasing excitation frequency. The system follows quasi-periodic motion between the chaotic regions, as shown in Figure 12e to Figure 12h. How the systems enter chaos is different. The system enters chaotic region A1 through quasi-periodic motion, as shown in Figure 12c,d. However, the system enters chaotic regions A2 and A3 by way of crisis, as shown in Figure 12i–l.

**Figure 12.** Evolution process of system response for Ω ∈ [0.5, 1.3].

As shown in Figure 11, the main components of the spectrum of gear system are composed of the meshing frequency and its harmonics. In the chaotic region, the frequency spectrum is a continuous line. The value range of the chaotic regions is also inconsistent with Figure 10.

Chaos means the uncontrollability and unpredictability of the system movement, which aggravates the vibration and noise of the system. In practice, such movement is to be avoided by appropriate measures. Through the above analysis of bifurcation characteristics and frequency analysis, the frequency regions and critical bifurcation values of chaotic motions and motions of different periods are obtained. Hence, the desired motion state can be obtained artificially.

#### *4.2. Effect of Gear Backlash*

Gear backlash is one of the main nonlinearities of the system. It generally changes the wear and the deformation of components of the system. The excitation frequency was set as a constant value (Ω *=* 1.2), and the other parameters remained unchanged in the study of bifurcation characteristics of the system.

The bifurcation diagram and three-dimensional frequency spectrum are shown in Figures 13 and 14, respectively. Both cases have the same trend in motion state transitions. Several frequency jumps exist in the bifurcation diagram, for example, at *b* equals to 11.3. Various motion patterns such as single periodic motion, multi-periodic motion, quasiperiodic motion, and chaotic motions are seen within the rotational speed range. There are mainly four chaotic regions in the bifurcation diagram, namely B1, B2, B3, and B4. The way the system enters a chaotic region is different. The system enters B1, B2, and B4 through quasi-periodic motion, as shown in Figure 15b,c,e,f,j,k. However, the system enters B3 through the period-doubling route, as shown in Figure 15h,i.

**Figure 13.** Bifurcation diagram under different gear backlash.

**Figure 14.** Three-dimensional frequency spectrum under different gear backlash.

**Figure 15.** Evolution process of system response at the range of b∈ [17, 75].

#### **5. Experimental Validation**

A single-stage high-contact-ratio gear test rig is designed and developed to validate the dynamic model. Dynamic characteristics under different speed conditions to imitate real-life applications are studied. The test rig is set up at the Anhui Digital Design and manufacturing laboratory, Hefei University of Technology, China. It can measure the acceleration, acting load, and rotational speed of the system. The gearbox is coupled to the servo motor, which has a maximum speed of 3000 rpm. The prime mover is Delta make with a power capacity of 0.75 kW. ZHY-6001 piezoelectric accelerometers measure the vibration signals of the system. Two DYN-200 torque transducers track the rotational speed and torque of the gear system. A data acquisition system ZHKJ-1001 with six channels is used. The sampling rate for the experimental trials is set at 10 kHz. The experimental setup and test gears are shown in Figure 16. High-contact-ratio gear is a non-standard part. Thus, the test gears are machined by slow wire cutting. Tooth spalling defects in the gears are generated with the electric mill.

**Figure 16.** Experimental setup.

The comparison between the experimental and simulated signals in the time-domain and frequency-domain is presented in Figure 17. The gear vibration signal contains a lot of interference noise caused by the operation of mechanical equipment. The white noise of the ambient background appears in its full frequency band. In this paper, the moving average method is used to denoise the collected experimental signals. The rotational frequency of pinion and gear is *fp* (*fp = np*/60 = 10 Hz) and *fg* (*fg* = 8.7 Hz), respectively. The pinion's meshing frequency and meshing are *fm* (*fm = zp* ∗ *fp =* 270 Hz) and *Tp* (*Tp =* 1/*fm =* 0.0037 s). The characteristic fault frequency is *fs*. It is the reciprocal of the failure period and can expressed as follows.

$$fsf = \frac{1}{\chi\_{p2} - \chi\_{p1}} \text{ or } fs = \frac{1}{\chi\_{p3} - \chi\_{p2}} \tag{35}$$

As shown in Table 3, there is little difference between the fault frequency calculated from the abscissa data of three points (*P*1, *P*2, *P*3) in Figure 17a and the data of the three points (*P*1, *P*2, *P*3) in Figure 17b.

**Figure 17.** Comparison between the experimental and simulated signals in horizontal direction. (**a**) Time−domain of experimental signal, (**b**) Time−domain of simulated signal, (**c**) Frequency spectrum of the experimental signal, (**d**) Frequency spectrum of the simulated signal.


Figure 17a,b show the time-domain experimental and simulated horizontal displacement signal. It can be observed that there is a certain difference in amplitude between the two cases, which is caused by experimental errors. Compared to experimental results, the simulated time-domain signal has a clear periodic impact to *Zp* ∗ *Tp*. The frequency response characteristic of the system under the two conditions is presented in Figure 17c,d. The frequency spectrum is mainly composed of gear mesh frequency and its harmonics. Peaks appear in the frequency spectrum at *n* ∗ *fm*, where *n* is a positive integer. However, the sideband structures can be seen clearly in both cases, as indicated by the arrow. The vibration spectrum of the spalled high-contact-ratio gear system is primarily characterized by the gear mesh frequency and its harmonics, and the sidebands induced by the modulation phenomenon. The results of simulated signal agree with that of experimental signals, which shows the reliability of the dynamic model proposed in this paper.

#### **6. Conclusions**

A novel approach for calculating the meshing stiffness of gear with tooth spalling defect is proposed. Simultaneously, a gear dynamic model is developed incorporating the time-varying mesh stiffness, fractal backlash and time-varying friction. The bifurcation characteristic of the gear system is acquired through the gear dynamic model. The vibration signatures of the gear system obtained experimentally are used to validate the dynamic model. Significant contributions of the study are summarized as follows:


**Author Contributions:** Z.C. (Writing—original draft: Lead); K.H. (Project administration: Lead); Y.X., (Data processing: Lead); M.S. (Software: Lead). All authors have read and agreed to the published version of the manuscript.

**Funding:** This study was supported by the National Natural Science Foundation of China (51775156), Natural Science Foundation of Anhui Province of China (1908085QE228), Key Research and Development Project of Anhui Province (202004h07020013), the Fundamental Research Funds for the Central Universities of China (PA2020GDSK0091), and the Horizontal Cooperation Project (0045022007).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The data that support the findings of this study are available from the corresponding author upon reasonable request.

**Conflicts of Interest:** The authors state that they have no conflicting interest.

#### **References**

