Dynamic Characteristics Analysis of Gear-Bearing System Considering Dynamic Wear with Flash Temperature

Copyright: © 2021 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/). 1 School of Mechanical Engineering and Automation, Northeastern University, Shenyang 110819, China; 1910093@stu.neu.edu.cn (J.X.); 1810123@stu.neu.edu.cn (R.C.); 1910094@stu.neu.edu.cn (Z.Y.); 1810124@stu.neu.edu.cn (H.Y.) 2 MICAS, Department of Electrical Engineering ESAT, Katholieke Universiteit Leuven (KU LEUVEN), 3001 Leuven, Belgium; linlin.wang@student.kuleuven.be * Correspondence: xpli@me.neu.edu.cn


Introduction
With the development of science and technology, clean energy has become an integral part of society. Wind energy is a kind of renewable clean energy, and wind power generation is an important part of non-fossil energy power generation. Furthermore, wind power generation plays a key role in reducing the use of fossil energy and alleviating environmental pollution. The gear transmission system, which is in the gearbox of the wind power generation, has an important impact on the stability of the system, and many scholars have carried out modeling analysis on it [1,2]. Based on a 1.5 MW wind turbine, Chen et al. [3] established a planetary gear system considering random wind speed excitation, and they studied the change of system response under the random component of comprehensive error. Zhao et al. [4] established a nonlinear dynamic model of the two-stage planetary gear and one-stage parallel axial gear transmission. Sun et al. [5] established an 8-DOF dynamic differential equation, and the dynamic responses caused by progressive tooth wear was studied. Wang et al. [6] transformed the semi-custom system with a rigid displacement model into a custom system. Then, the dynamic characteristics of a one-stage planetary gear transmission system and two-stage parallel axial gear transmission system were analyzed. Xu et al. [7] divided the gearbox of the wind power system into a transmission stiffness are introduced to calculate parameters (initial tooth surface wear, dynamic tooth surface wear, and tooth surface temperature). Then, dynamic characteristics of the gearbearing system with the initial tooth surface wear are studied. In Section 3, the meshing stiffness, considering tooth surface temperature, is established based on the flash temperature theory. Then, the meshing forces of gear bearing considering the backlash, which has fractal characteristics, are calculated. Furthermore, the influence of meshing force and friction coefficient on the tooth surface temperature and meshing stiffness is analyzed. In Section 4, based on the tooth surface temperature model and the initial tooth surface wear, the dynamic wear model of the tooth surface is established. Furthermore, the effects of the tooth surface initial wear, the dynamic wear, the tooth contact temperature, and the damping ratio on the dynamic characteristics of the gear-bearing system are analyzed.

System Model
The gear-bearing system in references [14,16] is taken as the research object, and the model diagram is shown in Figure 1. The meshing stiffness, considering the tooth surface flash temperature and the backlash with tooth surface wear, are calculated. Then, the dynamic wear model of a gear is constructed, and the interaction between the surface fractal dimension, temperature, and wear is analyzed.
In the model, the gear is an involute spur gear, and the bearing is a sliding bearing. Considering the time-varying meshing stiffness, backlash, tooth surface friction, comprehensive transmission error, nonlinear oil film force, and radial vibrations of the shaft and the torsional vibration of two gears, the coordinate system xiOiyi (i = 1 or 2) is established, in which O1 and O2 are the geometric centers of the two bearings, Oj1 and Oj2 are the centers of the two transmission shafts, and Op and Og are the geometric centers of the driving and driven gears, respectively. To calculate parameters (initial tooth surface wear, dynamic tooth surface wear, and tooth surface temperature), the gear-bearing systems considering different backlash and meshing stiffness are defined. Then, the gear-bearing system considering backlash with fractal characteristics is called System Ⅰ, through which the meshing force is calculated; the gear-bearing system considering backlash with initial tooth surface wear is called System Ⅱ, through which the system response under a wear condition is analyzed; and the gear-bearing system considering the dynamic wear model is called System Ⅲ.
In this paper, the backlash considering the initial tooth surface wear (Section 4.1) is calculated, and it is substituted into the gear-bearing system (System II). We are using the parameters μ = 0.3, ki = 1, ζh = 0.01, ζ1 = ζ2 = ζp = ζg = 0.01, ep = eg = 0.01, e = 0.1, Fm = 0.105, Fa In the model, the gear is an involute spur gear, and the bearing is a sliding bearing. Considering the time-varying meshing stiffness, backlash, tooth surface friction, comprehensive transmission error, nonlinear oil film force, and radial vibrations of the shaft and the torsional vibration of two gears, the coordinate system x i O i y i (i = 1 or 2) is established, in which O 1 and O 2 are the geometric centers of the two bearings, O j1 and O j2 are the centers of the two transmission shafts, and O p and O g are the geometric centers of the driving and driven gears, respectively.
To calculate parameters (initial tooth surface wear, dynamic tooth surface wear, and tooth surface temperature), the gear-bearing systems considering different backlash and meshing stiffness are defined. Then, the gear-bearing system considering backlash with fractal characteristics is called System I, through which the meshing force is calculated; the gear-bearing system considering backlash with initial tooth surface wear is called System II, through which the system response under a wear condition is analyzed; and the gear-bearing system considering the dynamic wear model is called System III.
In this paper, the backlash considering the initial tooth surface wear (Section 4.1) is calculated, and it is substituted into the gear-bearing system (System II). We are using the parameters µ = 0.3, k i = 1, ζ h = 0.01, ζ 1 = ζ 2 = ζ p = ζ g = 0.01, e p = e g = 0.01, e = 0.1, F m = 0.105, F a = 0.01, and D = 1.5, and the system component parameters are shown in Tables 1-3 in reference [16]. The Runge-Kutta method is used to calculate, and the dynamic equations are shown in Equation (8) in reference [16]. The bifurcation diagram, the top Lyapunov exponent diagram, and the three-dimensional spectrum diagram are obtained, respectively, as shown in Figure 2.
states, such as nT-periodic, quasi-periodic, and chaotic in the process of the relative speed changing from 0.5 to 4. In Figure 2a, region A represents the chaotic motion state and its corresponding top Lyapunov exponent is greater than zero; region B represents the nTperiodic motion state and its corresponding top Lyapunov exponent is less than zero; other regions are quasi-periodic motion states between the chaotic and periodic motion state. In region A, there is a sudden change of the top Lyapunov exponent, and the system motion states have also changed correspondingly, which also indicates that the system is unstable in this speed range. At the same time, it can be observed from Figure 2b that the system contains 1/3fm, 1/2fm, 2/3fm, fm, and 2fm harmonic responses, which can prove that the system has rich nonlinear factors.

Calculation of Tooth Surface Flash Temperature
When the gear is running under high speed and heavy load, an abundant of heat will be generated due to the friction consumption of the tooth surface. Then, it will increase the tooth surface contact temperature and affect the stability as well as the service life of the gear system. The influence of contact temperature on the stiffness and dynamic characteristics of the system was analyzed in reference [11], but the gear system is simple. Based on the multi-degree-of-freedom gear-bearing system, the paper not only analyzes the influence of temperature on stiffness and system response but also studies the dynamic interaction between temperature and wear.
The tooth contact temperature consists of two parts: the base temperature ΔM and the instantaneous flash temperature Δf of the tooth surface [11], which is given by Equation (1).
According to the flash temperature theory [11], Δf can be expressed as: From Figure 2a, it can be observed that the system has experienced a series of motion states, such as nT-periodic, quasi-periodic, and chaotic in the process of the relative speed changing from 0.5 to 4. In Figure 2a, region A represents the chaotic motion state and its corresponding top Lyapunov exponent is greater than zero; region B represents the nT-periodic motion state and its corresponding top Lyapunov exponent is less than zero; other regions are quasi-periodic motion states between the chaotic and periodic motion state. In region A, there is a sudden change of the top Lyapunov exponent, and the system motion states have also changed correspondingly, which also indicates that the system is unstable in this speed range. At the same time, it can be observed from Figure 2b that the system contains 1/3f m , 1/2f m , 2/3f m , f m , and 2f m harmonic responses, which can prove that the system has rich nonlinear factors.

Calculation of Tooth Surface Flash Temperature
When the gear is running under high speed and heavy load, an abundant of heat will be generated due to the friction consumption of the tooth surface. Then, it will increase the tooth surface contact temperature and affect the stability as well as the service life of the gear system. The influence of contact temperature on the stiffness and dynamic characteristics of the system was analyzed in reference [11], but the gear system is simple. Based on the multi-degree-of-freedom gear-bearing system, the paper not only analyzes the influence of temperature on stiffness and system response but also studies the dynamic interaction between temperature and wear.
The tooth contact temperature consists of two parts: the base temperature ∆ M and the instantaneous flash temperature ∆ f of the tooth surface [11], which is given by Equation (1).
According to the flash temperature theory [11], ∆ f can be expressed as: where u is the temperature rise coefficient, and µ m is the friction coefficient. (For the convenience of subsequent analysis, the friction coefficient µ m in this section and the friction coefficient µ in the friction force of the gear-bearing system are represented by different symbols). Here, f e is the normal load on the tooth surface per unit tooth width, v i (i = p, g) is the tangential speed of the tooth surface, which is given by Equation (3) [21], g i is the heat conduction coefficient, ρ i is the material density, B is the half width of the contact belt, and c i is the specific heat capacity.
where ω p and ω g are the rotating speeds of the driving and driven gears, respectively, d 1 and d 2 are the diameter of the graduated circle of the driving and driven gears, α is the pressure angle of the graduated circle (α = 20 • ), and y is the distance from the meshing point to the pitch.

Deformation Caused by the Tooth Surface Flash Temperature
The increasing of tooth surface temperature will change the contact state of teeth, and then, it will lead to change of the tooth surface profile [11], which is given by: where ∆ f (t) is the difference between the tooth surface contact temperature ∆ B (t) and the base temperature ∆ M , s bi (i = p, g) is the thickness of the driving and driven gears, λ f is the linear expansion coefficient, α ki (i = p, g) is the pressure angle of the addendum circle of the driving and driven gears, and invα is the involute function. u bi (i = p, g) is the thermal deformation of the base circle of the driving and driven gears when the system works stably, given by: where ∆(r 0i ) (i = p, g) is the temperature of transmission shafts of the driving and driven gears, ∆(r bi ) (i = p, g) is the temperature of the base circle of the driving and driven gears, and r 0i is the shaft radius of the driving and driven gears.

Meshing Stiffness Caused by the Tooth Surface Flash Temperature
Due to the gear transmission, an abundant heat is generated, and it leads to a change of the tooth surface contact conditions. Then, the meshing force and deformation of the tooth surface also change. Thus, the meshing stiffness fluctuates. According to the Hertz contact theory, the change of meshing stiffness caused by temperature can be calculated as: where k Detai (i = p, g) is the stiffness of the driving and driven gears caused by temperature, b is the tooth width, and P is the meshing force. In this paper, the meshing force P is calculated by Equation (7). In addition, the displacement and the velocity in Equation (7) are calculated by system I. Thus, the meshing force is related to fractal dimension D.
Due to the manufacturing and installation errors of gears, the backlash will produce certain fluctuations and uncertainties. In this paper, the fractal characteristics of backlash are used to describe the influence of errors, which is caused by the gear manufacturing and installation on backlash. The fractal characteristics of backlash are expressed by Equation (8) [17,22]. where L is the sampling length, G is the characteristic parameter, D is the fractal dimension, and γ is the spatial frequency of the contour. Since the variable in Equation (8) is time t, the sampling length L = T (Period). Then, we take parameters as µ = 0.3, k i = 1, ζ h = 0.01, ζ 1 = ζ 2 = ζ p = ζ g = 0.01, e p = e g = 0.01, e = 0.1, F m = 0.105, F a = 0.01, D = 1.43, 1.5, and 1.6. Based on System I and Equation (7), the meshing force is calculated. Since the meshing force changes periodically with time, the meshing force in one period is selected for comparative analysis (the 2716th period is adopted in this paper). Then, the backlash diagrams, the system phase diagrams, the Poincare section diagrams, and the meshing force diagrams are obtained by System I (s 3-5).
From Figure 3, it can be observed that when D = 1.43, the amplitude of the backlash with fractal characteristics is high (Figure 3a). Then, the Poincare section diagram (red region in Figure 3b) of the system is a set of three regiment points, and it indicates that the system is in the quasi-three-periodic motion state. In addition, the amplitude of meshing force ( Figure 3c) is large, and the fluctuation is obvious. In Figure 4, the fractal dimension D is 1.5. It can be observed that the amplitude of the backlash decreases and the set of three regiment points (red region in Figure 4b) begins to contract, which means that the system begins to change from a quasi-periodic motion state to a stable periodic motion state. Meanwhile, the amplitude of the meshing force ( Figure 4c) reduces, and the curve of it tends to be smooth. As the fractal dimension D increases, it can be observed from Figure 5 that when D = 1.6, the amplitude of the backlash ( Figure 5a) is low, and the system is in the stable three periodic motion state (Figure 5b). At the same time, the meshing force is small, and the curve is smooth with little volatility. From Figure 6, it can be observed that with the increase of D, the amplitude of the meshing force decreases, and the tooth surface flash temperature also decreases. It shows that the meshing force of the gear system of the backlash with small fluctuation is relatively stable and then the tooth surface temperature is also relatively stable. Figure 7 shows changing of the meshing stiffness under different fractal dimensions D. However, the meshing force and deformation caused by the temperature are both changing with different fractal dimensions D. Therefore, the variation law of stiffness is not analyzed, and only the influence of temperature in the dynamic wear model (Section 4) on the system response is studied.    From Figure 6, it can be observed that with the increase of D, the amplitude of the meshing force decreases, and the tooth surface flash temperature also decreases. It shows that the meshing force of the gear system of the backlash with small fluctuation is relatively stable and then the tooth surface temperature is also relatively stable. Figure 7 shows changing of the meshing stiffness under different fractal dimensions D. However, the meshing force and deformation caused by the temperature are both changing with different fractal dimensions D. Therefore, the variation law of stiffness is not analyzed, and only the influence of temperature in the dynamic wear model (Section 4) on the system response is studied.      Above all, with increasing the fractal dimension D, the amplitude of the backlash with fractal characteristics reduces. Then, according to the change of the system response, the system changes from an unstable quasi-periodic motion state to the stable periodic motion state, and the meshing force gradually decreases. However, the change trend of the meshing force in one period is similar.
The meshing forces are substituted into Equations (2)-(6) to calculate the tooth surface flash temperature and the stiffness k Detai of the driving and driven gears. The specific parameters are the same in reference [11]. According to reference [23], the comprehensive stiffness generated by temperature is shown by Equation (9).
According to reference [16], the stiffness can be given in the form of Fourier expansion, and combined with reference [23], the meshing stiffness, considering the tooth surface flash temperature, can be calculated by Equation (10). Then, the curves of the tooth surface flash temperature and dimensionless meshing stiffness are obtained by Figures 6 and 7.  From Figure 6, it can be observed that with the increase of D, the amplitude of the meshing force decreases, and the tooth surface flash temperature also decreases. It shows that the meshing force of the gear system of the backlash with small fluctuation is relatively stable and then the tooth surface temperature is also relatively stable. Figure 7 shows changing of the meshing stiffness under different fractal dimensions D. However, the meshing force and deformation caused by the temperature are both changing with different fractal dimensions D. Therefore, the variation law of stiffness is not analyzed, and only the

Dynamic Wear Model of the Tooth Surface
Tooth surface wear is a dynamic process and the contact temperature, as well as system excitation, have an important impact on it. At the same time, wear can also change the conditions of the tooth surface contact, which has an impact on the stiffness and tooth surface temperature. Furthermore, the tooth surface wear is generally calculated by Archard theory [24], which is given by: where V is the wear volume of the material, s is the relative sliding distance, W is the positive pressure of the contact surface, H is the hardness of the material, and K is the wear coefficient. The meshing process of the gear is discretized into the movement of each meshing point and the wear amount at different meshing point is calculated, which is shown by: It can be observed that the calculation of the tooth surface wear mainly involves three parameters: wear coefficient k, sliding distance s, and contact pressure P. Then, the wear coefficient k is 5 × 10 −16 m 2 /N, the contact pressure is calculated by Equation (7), and the sliding distance is calculated by the equation in reference [21], which is given by:

Dynamic Wear Model of the Tooth Surface
Tooth surface wear is a dynamic process and the contact temperature, as well as system excitation, have an important impact on it. At the same time, wear can also change the conditions of the tooth surface contact, which has an impact on the stiffness and tooth surface temperature. Furthermore, the tooth surface wear is generally calculated by Archard theory [24], which is given by: where V is the wear volume of the material, s is the relative sliding distance, W is the positive pressure of the contact surface, H is the hardness of the material, and K is the wear coefficient.
The meshing process of the gear is discretized into the movement of each meshing point and the wear amount at different meshing point is calculated, which is shown by: It can be observed that the calculation of the tooth surface wear mainly involves three parameters: wear coefficient k, sliding distance s, and contact pressure P. Then, the wear coefficient k is 5 × 10 −16 m 2 /N, the contact pressure is calculated by Equation (7), and the sliding distance is calculated by the equation in reference [21], which is given by: where v p and v g are the sliding speeds of the driving and driven gears, respectively, which are shown by Equation (3), and a H is the half-width of the contact zone. In order to analyze and calculate the dynamic wear model more accurately, a certain initial wear amount (H initial ) is adopted in the numerical calculation of the system. The calculation flow of dynamic wear is shown in Figure 8. Then, the specific steps are as follows: (1) The fractal theory is used to describe the backlash in order to represent the error caused by gear manufacturing and installation. Then, the meshing force P(t, D) is obtained through System I. (2) Based on the Archard theory and flash temperature theory, the initial tooth surface wear and the flash temperature are calculated. (3) The backlash b h (1) with initial wear and the stiffness k Detai (t, D, 1) with flash temperature in the first round are calculated, which represent the initial conditions of the gear system. Then, they are substituted into System III. (4) According to System III, the dynamic meshing force P dynamic (t, D, n + 1) can be obtained. (5) The dynamic tooth surface wear H dynamic (t, D, n + 1) and the flash temperature in the (n + 1)th round are calculated. (6) The backlash b h (n + 1) in the (n + 1)th round is obtained by adding the dynamic wear H dynamic (t, D, n+1) to the backlash b h (n) in the (n)th round. Then, the stiffness with flash temperature k Detai (t, D) in the (n + 1)th round is obtained based on the Hertz theory. (7) The backlash b h (n + 1) and the stiffness with flash temperature k Detai (t, D) in the (n + 1)th round are substituted into System III. Then, go back to step 4 until the total time is reached.
In order to analyze and calculate the dynamic wear model more accurately, a certain initial wear amount (Hinitial) is adopted in the numerical calculation of the system. The calculation flow of dynamic wear is shown in Figure 8. Then, the specific steps are as follows: (1) The fractal theory is used to describe the backlash in order to represent the error caused by gear manufacturing and installation. Then, the meshing force P(t, D) is obtained through System Ӏ. (2) Based on the Archard theory and flash temperature theory, the initial tooth surface wear and the flash temperature are calculated. (3) The backlash bh(1) with initial wear and the stiffness kDetai (t, D, 1) with flash temperature in the first round are calculated, which represent the initial conditions of the gear system. Then, they are substituted into System Ⅲ. (4) According to System Ⅲ, the dynamic meshing force Pdynamic(t, D, n + 1) can be obtained. (5) The dynamic tooth surface wear Hdynamic(t, D, n + 1) and the flash temperature in the (n + 1)th round are calculated. (6) The backlash bh(n + 1) in the (n + 1)th round is obtained by adding the dynamic wear Hdynamic(t, D, n+1) to the backlash bh(n) in the (n)th round. Then, the stiffness with flash temperature kDetai(t, D) in the (n + 1)th round is obtained based on the Hertz theory. (7) The backlash bh(n + 1) and the stiffness with flash temperature kDetai(t, D) in the (n + 1)th round are substituted into System III. Then, go back to step 4 until the total time is reached.

Calculation of Initial Wear
Based on System Ӏ, the meshing forces with different fractal dimensions D are calculated. Then, based on Equation (12), the tooth surface wear Δh in one period T is obtained. It is assumed that the amount of wear per round is the same and the tooth surface wear of 10 6 rounds is calculated, that is, Δh × 10 6 . The accumulated wear of the tooth surface is

Calculation of Initial Wear
Based on System I, the meshing forces with different fractal dimensions D are calculated. Then, based on Equation (12), the tooth surface wear ∆h in one period T is obtained. It is assumed that the amount of wear per round is the same and the tooth surface wear of 10 6 rounds is calculated, that is, ∆h × 10 6 . The accumulated wear of the tooth surface is shown in Figure 9; furthermore, it is called the initial tooth surface wear (H initial ) in this paper. From Figure 9, it can be observed that with increase of D, the amplitude and the fluctuation of the tooth surface wear decrease.

Influence of Dynamic Wear on the System Response
Based on backlash considering the initial tooth surface wear (H initial ) and the meshing stiffness considering the tooth surface flash temperature, the gear-bearing system considering the dynamic wear model (System III) is calculated. Taking the parameters as µ = 0.3, k i = 1, ζ h = 0.01, ζ 1 = ζ 2 = ζ p = ζ g = 0.01, µ m = 0.06, e p = e g = 0.01, e = 0.1, F m = 0.105, F a = 0.01, and ω = 1.5, the dynamic characteristics of the gear system with 100 rounds are calculated and simulated. Then, the phase diagrams, the Poincare section diagrams, and the time domain diagrams of the system are obtained, as shown in Figures 10-12. Mathematics 2021, 9, x FOR PEER REVIEW 10 of 17 shown in Figure 9; furthermore, it is called the initial tooth surface wear (Hinitial) in this paper. From Figure 9, it can be observed that with increase of D, the amplitude and the fluctuation of the tooth surface wear decrease.

Influence of Dynamic Wear on the System Response
Based on backlash considering the initial tooth surface wear (Hinitial) and the meshing stiffness considering the tooth surface flash temperature, the gear-bearing system considering the dynamic wear model (System Ⅲ) is calculated. Taking the parameters as μ = 0.

Influence of Dynamic Wear on the System Response
Based on backlash considering the initial tooth surface wear (Hinitial) and the meshing stiffness considering the tooth surface flash temperature, the gear-bearing system considering the dynamic wear model (System Ⅲ) is calculated. Taking the parameters as μ = 0.

Influence of Dynamic Wear on the System Response
Based on backlash considering the initial tooth surface wear (Hinitial) and the meshing stiffness considering the tooth surface flash temperature, the gear-bearing system considering the dynamic wear model (System Ⅲ) is calculated. Taking the parameters as μ = 0.     It can be observed that due to the different initial tooth surface wear (H initial ), the system responses are slightly different after the same numbers of periods. When the fractal dimension D is 1.43, there are less blanks in the phase diagram (Figure 10a), which indicates that the phase trajectories of the system are widely distributed in the space; regions 1, 2, and 3 in the Poincare section diagram (Figure 10b) are dispersive, and the regularity of motion is weak; at the same time, in Figure 10c, the displacement in region 5 (the timedomain diagram) fluctuates greatly, and all the above phenomena indicate that the system is unstable. When the fractal dimension D is 1.5, the blank area in Figure 11a becomes larger, which indicates that the concentration region of phase trajectories is smaller than that of D = 1.43; this phenomenon can also be obtained from the Poincare section diagram (Figure 11b), where regions 1, 2, and 3 begin to contract; at the same time, the displacement of region 5 in Figure 11c also shows a periodic change with small fluctuation, which proves that the system tends to be stable. As the fractal dimension continues to increase, the initial wear amount is low when the fractal dimension D is 1.6. The difference between Figures 11a and 12a is not obvious, the set of points in Figure 12b continues to contract, and the changing of displacement in Figure 12c basically conforms to periodic motion. Furthermore, all these indicate that the motion state of the system is more stable. It can be observed that due to the different initial tooth surface wear (Hinitial), the system responses are slightly different after the same numbers of periods. When the fractal dimension D is 1.43, there are less blanks in the phase diagram (Figure 10a), which indicates that the phase trajectories of the system are widely distributed in the space; regions 1, 2, and 3 in the Poincare section diagram (Figure 10b) are dispersive, and the regularity of motion is weak; at the same time, in Figure 10c, the displacement in region 5 (the timedomain diagram) fluctuates greatly, and all the above phenomena indicate that the system is unstable. When the fractal dimension D is 1.5, the blank area in Figure 11a becomes larger, which indicates that the concentration region of phase trajectories is smaller than that of D = 1.43; this phenomenon can also be obtained from the Poincare section diagram (Figure 11b), where regions 1, 2, and 3 begin to contract; at the same time, the displacement of region 5 in Figure 11c also shows a periodic change with small fluctuation, which proves that the system tends to be stable. As the fractal dimension continues to increase, the initial wear amount is low when the fractal dimension D is 1.6. The difference between Figures 11a and 12a is not obvious, the set of points in Figure 12b continues to contract, and the changing of displacement in Figure 12c basically conforms to periodic motion. Furthermore, all these indicate that the motion state of the system is more stable.
In conclusion, the fractal characteristics of backlash represent the uncertainty and amplitude. That is, the amplitude and the uncertainty of backlash decrease with the increase of the fractal dimension D. As a result, the initial tooth surface wear (Hinitial) becomes smaller. Then, the system response (system Ⅲ) becomes stable gradually. Thus, the stable backlash can slow down the tooth surface wear.
The calculation of dynamic wear is a cumulative process. Since the tooth surface wear of the previous period will affect the wear of the next period, the wear amount and the system response of each period are different. Then, the system responses (system Ⅲ) of 150 rounds are calculated by using the above parameters. Since the front periods will be affected by the transient response, the dynamic cumulative wear (Hdynamic) of different rounds with D = 1.43 is analyzed. Then, the dynamic cumulative wear amount (Hdynamic) and the system responses are obtained, as shown in Figures 13-15.
In Figure 13, it can be observed that the dynamic cumulative wear (Hdynamic) increases with the rotation of the gear. Since the number of rounds calculated is small, the difference of tooth surface wear between rounds is not obvious. Thus, the system response (system III) is analyzed in order to study the difference of wear of different rounds. Figure 14a  In conclusion, the fractal characteristics of backlash represent the uncertainty and amplitude. That is, the amplitude and the uncertainty of backlash decrease with the increase of the fractal dimension D. As a result, the initial tooth surface wear (H initial ) becomes smaller. Then, the system response (system III) becomes stable gradually. Thus, the stable backlash can slow down the tooth surface wear.
The calculation of dynamic wear is a cumulative process. Since the tooth surface wear of the previous period will affect the wear of the next period, the wear amount and the system response of each period are different. Then, the system responses (system III) of 150 rounds are calculated by using the above parameters. Since the front periods will be affected by the transient response, the dynamic cumulative wear (H dynamic ) of different rounds with D = 1.43 is analyzed. Then, the dynamic cumulative wear amount (H dynamic ) and the system responses are obtained, as shown in Figures 13-15.
In Figure 13, it can be observed that the dynamic cumulative wear (H dynamic ) increases with the rotation of the gear. Since the number of rounds calculated is small, the difference of tooth surface wear between rounds is not obvious. Thus, the system response (system III) is analyzed in order to study the difference of wear of different rounds. Figure 14a-c represent the phase diagrams obtained by the corresponding displacement and velocities with 75-100, 100-125, and 125-150 rounds, respectively, while Figure 15a-c represent the Poincare section diagrams obtained by the corresponding displacement and velocities with 75-100, 100-125, and 125-150 rounds, respectively. Then, it can be observed from Figure 14 that the blank part in the phase diagram decreases with rotation of the gear, which indicates that the displacement and velocity of the system are widely distributed. From Figure 15, it can be observed that with the rotation of the gear, the concentration of regions 1 and 2 becomes worse, the points in Poincare diagram begin to disperse gradually, and region 3 gradually extends to the distance. All these show that the stability of the system becomes worse, and the accuracy of the dynamic wear model is also proven.
with 75-100, 100-125, and 125-150 rounds, respectively. Then, it can be observed from Figure 14 that the blank part in the phase diagram decreases with rotation of the gear, which indicates that the displacement and velocity of the system are widely distributed. From Figure 15, it can be observed that with the rotation of the gear, the concentration of regions 1 and 2 becomes worse, the points in Poincare diagram begin to disperse gradually, and region 3 gradually extends to the distance. All these show that the stability of the system becomes worse, and the accuracy of the dynamic wear model is also proven.   Figure 14 that the blank part in the phase diagram decreases with rotation of the gear, which indicates that the displacement and velocity of the system are widely distributed. From Figure 15, it can be observed that with the rotation of the gear, the concentration of regions 1 and 2 becomes worse, the points in Poincare diagram begin to disperse gradually, and region 3 gradually extends to the distance. All these show that the stability of the system becomes worse, and the accuracy of the dynamic wear model is also proven.   Figure 14 that the blank part in the phase diagram decreases with rotation of the gear, which indicates that the displacement and velocity of the system are widely distributed. From Figure 15, it can be observed that with the rotation of the gear, the concentration of regions 1 and 2 becomes worse, the points in Poincare diagram begin to disperse gradually, and region 3 gradually extends to the distance. All these show that the stability of the system becomes worse, and the accuracy of the dynamic wear model is also proven.

Influence of µ m on the System Response
In the dynamic wear model, in addition to the initial tooth surface wear (H initial ) and the number of rounds n, the meshing stiffness considering the tooth surface flash temperature also has an important influence on the system response. From Equation (2), it can be observed that the friction coefficient µ m of the tooth surface and the relative speed have important effects on the tooth surface flash temperature and then, they affect the stiffness and system response. Due to the difference of meshing forces at different relative speeds, System I is chosen to be calculated. Taking the parameters as µ = 0.3, k i = 1, e g = 0.01, ζ h = 0.01, ζ 1 = ζ 2 = ζ p = ζ g = 0.01, e p = e g = 0.01, e = 0.1, F m = 0.105, F a = 0.01, and D = 1.43, the meshing forces at different relative speeds are obtained. Then, the mean values of the meshing forces are shown in Table 1. Taking the parameters as ω = 0.5-4, µ m = 0.02-0.1, the tooth surface flash temperature is calculated, which is shown in Figure 16.  From Figure 16, it can be observed that the abscissa is the relative spe system, and different lines represent the tooth surface flash temperature friction coefficients. Then, through increasing the relative speed and the cient, the tooth surface flash temperature increases gradually. Since the does not increase monotonically, there are some fluctuations in the chang flash temperature. Then, the influence of friction coefficient μm on the system response studied. Taking the parameters as μm = 0.02, 0.06, 0.1, they are substituted i and the Poincare section diagrams are obtained, which are shown in Figure  It can be observed that when μm = 0.02, the system response is shown a regiment points, and the system is in a quasi-periodic motion state at this ti points in regions 1 and 2 is relatively close, and the system is relatively stab 0.06, a small number of points in regions 1 and 2 begin to diffuse outward. not obvious, that also indicates that the system transforms from a stable state. When μm = 0.1, many points in regions 1 and 2 break away from the s outward, and then, the system tends to be unstable and change to a chaotic It can be concluded that with an increase of the friction coefficient, the syste unstable. However, the trend is not very obvious. In order to conduct a mor obvious analysis, the root mean square diagrams of system displacement friction coefficients are calculated, as shown in Figure 18.
In the three cases of different initial wear amount (D = 1.43, 1.5, 1.6), rameters as μm = 0.02-0.12, D = 1.43, 1.5, and 1.6, the root mean square va displacement are calculated. It can be concluded that the tooth surface flas From Figure 16, it can be observed that the abscissa is the relative speed of the gear system, and different lines represent the tooth surface flash temperature with different friction coefficients. Then, through increasing the relative speed and the friction coefficient, the tooth surface flash temperature increases gradually. Since the meshing force does not increase monotonically, there are some fluctuations in the changing process of flash temperature.
Then, the influence of friction coefficient µ m on the system response (System III) is studied. Taking the parameters as µ m = 0.02, 0.06, 0.1, they are substituted into System III, and the Poincare section diagrams are obtained, which are shown in Figure 17. rises with the increase of friction coefficient. Then, the root mean square value of the vibration displacement of the system amplifies, which indicates that the system tends to be unstable.   It can be observed that when µ m = 0.02, the system response is shown as a set of four regiment points, and the system is in a quasi-periodic motion state at this time. The set of points in regions 1 and 2 is relatively close, and the system is relatively stable. When µ m = 0.06, a small number of points in regions 1 and 2 begin to diffuse outward. Although it is not obvious, that also indicates that the system transforms from a stable to an unstable state. When µ m = 0.1, many points in regions 1 and 2 break away from the set and diffuse outward, and then, the system tends to be unstable and change to a chaotic motion state. It can be concluded that with an increase of the friction coefficient, the system tends to be unstable. However, the trend is not very obvious. In order to conduct a more detailed and obvious analysis, the root mean square diagrams of system displacement with different friction coefficients are calculated, as shown in Figure 18.

Influence of Damping Ratio on the System Response
The influence of important factors in the dynamic wear model on the sys sponse (System Ⅲ) is studied in the above sections. Then, in order to analyze the in detail, the influence of the damping ratio on the system response (System Ⅲ) lyzed. Taking the parameters as μm = 0.06, μ = 0.3, ki = 1, ζh = 0.01, ζ1 = ζ2 = ζp = ζg = = eg = 0.01, e = 0.1, Fm = 0.105, Fa = 0.01, and D = 1.43, they are substituted into Sys Then, the phase diagram and the Poincare section diagram are obtained, as shown ures 19-21. Figure 19 shows the Poincare diagram and the phase diagram with damping = 0.01. Combined with Figure 19a,b, it can be observed that the system has expe chaotic, quasi-double-periodic, quasi-three-periodic, quasi-periodic, quasi-four-p and 1T-periodic motion states; when ζ = 0.05 (Figure 20), the phase diagram cha the limit cycle, and the set of points in the Poincare section diagram begins to c Then, combined with Figure 20a,b, it can be observed that the system has experien periodic, 3T-periodic, 1T-periodic, and quasi-double-periodic and motion state pared to the motion states when ζ = 0.01, the system changes from unstable quasi-p and chaotic motion to stable periodic motion; when ζ = 0.1 (Figure 21), the system p stable periodic motion states. Therefore, under the influence of dynamic wear, the i of damping ratio can also accelerate the energy loss of the system and make the tend to be stable. In the three cases of different initial wear amount (D = 1.43, 1.5, 1.6), taking the parameters as µ m = 0.02-0.12, D = 1.43, 1.5, and 1.6, the root mean square values of system displacement are calculated. It can be concluded that the tooth surface flash temperature rises with the increase of friction coefficient. Then, the root mean square value of the vibration displacement of the system amplifies, which indicates that the system tends to be unstable.

Influence of Damping Ratio on the System Response
The influence of important factors in the dynamic wear model on the system response (System III) is studied in the above sections. Then, in order to analyze the system in detail, the influence of the damping ratio on the system response (System III) is analyzed. Taking the parameters as µ m = 0.06, µ = 0.3, k i = 1, ζ h = 0.01, ζ 1 = ζ 2 = ζ p = ζ g = 0.01, e p = e g = 0.01, e = 0.1, F m = 0.105, F a = 0.01, and D = 1.43, they are substituted into System III. Then, the phase diagram and the Poincare section diagram are obtained, as shown in Figures 19-21. Due to different coordinates of Figures 19-21, the transition of motion states is not obvious. The system with relative speed ω = 2.5 is taken for analysis, and the Poincare section diagrams with ζ = 0.01, 0.05 and 0.1 are obtained, as shown in Figure 22. Then, it can be observed that the system gradually changes from quasi four periodic to single periodic, and the system tends to be stable. -50 The relative velocity of meshing point -50 The relative velocity of meshing point  Figure 19 shows the Poincare diagram and the phase diagram with damping ratio ζ = 0.01. Combined with Figure 19a,b, it can be observed that the system has experienced chaotic, quasi-double-periodic, quasi-three-periodic, quasi-periodic, quasi-four-periodic and 1T-periodic motion states; when ζ = 0.05 (Figure 20), the phase diagram changes to the limit cycle, and the set of points in the Poincare section diagram begins to contract. Then, combined with Figure 20a,b, it can be observed that the system has experienced 2T-periodic, 3T-periodic, 1T-periodic, and quasi-double-periodic and motion states. Compared to the motion states when ζ = 0.01, the system changes from unstable quasi-periodic and chaotic motion to stable periodic motion; when ζ = 0.1 (Figure 21), the system presents stable periodic motion states. Therefore, under the influence of dynamic wear, the increase of damping ratio can also accelerate the energy loss of the system and make the system tend to be stable. -50 The relative velocity of meshing point -50 The relative velocity of meshing point The relative velocity of meshing point -20 The relative velocity of meshing point -5 -20 The relative velocity of meshing point The relative velocity of meshing point  -50 The relative velocity of meshing point -50 The relative velocity of meshing point The relative velocity of meshing point -20 The relative velocity of meshing point -5 -20 The relative velocity of meshing point The relative velocity of meshing point  Due to different coordinates of Figures 19-21, the transition of motion states is not obvious. The system with relative speed ω = 2.5 is taken for analysis, and the Poincare section diagrams with ζ = 0.01, 0.05 and 0.1 are obtained, as shown in Figure 22. Then, it can be observed that the system gradually changes from quasi four periodic to single periodic, and the system tends to be stable.

Conclusions
(1) With the increase of the fractal dimension D, the uncertainty and the amplitude of backlash decrease. Then, the meshing force is reduced, and the wear is small when the

Conclusions
(1) With the increase of the fractal dimension D, the uncertainty and the amplitude of backlash decrease. Then, the meshing force is reduced, and the wear is small when the gear runs for the same time. The influence on the dynamic wear model is reduced, and the stability of the system considering dynamic wear is enhanced. Therefore, the stable backlash can slow down the tooth surface wear.
(2) The friction coefficient and rotational speed of the tooth surface have an important influence on the temperature of the tooth surface. With the increase of the friction coefficient, the flash temperature and root mean square value of displacement increase, and the system tends to be unstable.
(3) With the increase of the damping ratio, the system changes from unstable quasiperiodic and chaotic motion to stable periodic motion. The increase of damping accelerates the energy loss of the system, makes the system tend to be stable, and has the effect of reducing wear. Institutional Review Board Statement: Not applicable.

Informed Consent Statement: Not applicable.
Data Availability Statement: Some data, models, or code generated or used during the study are available from the corresponding author by request.

Conflicts of Interest:
The authors declare no conflict of interest.