Fractal‐based dynamic response of a pair of spur gears considering microscopic surface morphology

The meshing surfaces of a gear pair are rough from a microscopic perspective and the surface topography will affect the dynamic response. To study the influence of real surface topography on the gear system dynamic performance, this paper establishes a 3‐degree of freedom transverse‐torsional dynamic model with regard to the morphology of the interface. By fractal theory, the expression of backlash between gears is modified based on the height of asperities. The time‐varying stiffness is calculated according to the fractal method rather than assuming a constant, which is more realistic. The dimensionless dynamic differential equations are established and solved with surface topography affected backlash function and time‐varying stiffness. The dynamic response of the gear system with respect to fractal dimension and fractal roughness is analyzed.


| INTRODUCTION
As a kind of widely used transmission mode, gears have been studied extensively. The dynamics of gear pair is the center and key content of gear research. Among the generally applied gear dynamic models, the gear pair is always simplified as two elastomeric contact bodies.
In the process of meshing, the backlash and meshing stiffness are the main nonlinear dynamic parameters. Accordingly, their accurate modeling is the key to solve the dynamic response of gears.
Traditionally, the backlash has been treated mathematically as a constant or a steady random number conforming to the normal distribution. 1 proposed a nonlinear backlash model composed of a third-order polynomial function, which was applied to a single degree of freedom (DOF) dynamic model. Considering the radial backlash as a control parameter, Sheng et al. 8 studied its effects on unstable gears and chaos of the system. Also, the wear of the surface may lead to the change of tooth thickness at the meshing region and finally the variation in backlash. In this respect, Onishchenko 9 took instantaneous contact temperature into consideration and proposed a tooth wear model. Feng et al. 10 established an updatable wear model in which system vibration was calculated considering the initial wear, which subsequently affects the current wear. The model formed a closed-loop update of wear.
However, only setting an appropriate constant for the backlash from a mathematical point of view has high randomness, which cannot reflect the influences of surface topography in reality.
In terms of meshing stiffness, it has become a consensus to adopt time-varying meshing stiffness, which has different calculation methods.
Kaharman and Singhr 11 had earlier introduced time-varying stiffness and studied its influence on gear system dynamics. Wang et al. 12 Figure 1, the roughness of tooth surfaces exists certainly. The height of microscopic asperities will influence the backlash of meshing gears and the surface features will affect the normal contact stiffness of the interface, which ultimately influences the mechanics and dynamic behavior of gears. Therefore, the study of gear dynamics needs to consider the influence of meshing surface morphology to obtain more accurate results.
In this study, considering the microscopic topography of meshing teeth surfaces, the nonlinear dynamic response of spur gear based on a transverse-torsional dynamic model is studied. Both the backlash function and meshing stiffness of the gear pair are derived with regard to the influence of surface topography by the fractal method.
The system dynamic performance under various fractal parameters that represent the surface topography is analyzed and discussed. The results show that morphology has a significant effect on the dynamic response.
where ȳ p is the displacement of the pinion center and ȳ g is the displacement of the gear center.
The torsional DOF is equivalent to the relative displacement of gears, that is where e t ( ) denotes the time-varying comprehensive error that originates from manufacturing, 3  Considering Equations (1) and (2), the differential equations governing the coupled transverse-torsional motion are expressed as where m e is the equivalent mass of the gear pair decided by is the rotational inertia of pinion and gear respectively, c m is the damping coefficient of meshing gears, that is where b is the backlash of the gear pair.

| Classical expressions of nonlinear backlash and meshing stiffness
In Equation ( The meshing stiffness of the gear pair is essentially composed of five parts: bending stiffness K b , shear stiffness K s , axial compression stiffness K a , the fillet-foundation deformation stiffness K f , and contact stiffness K h . One of the commonly used methods to calculate stiffness is the potential energy method, which treats gear tooth as a cantilever beam. Based on the principle of virtual work and considering the tooth profile, the work borne by each meshing point under the meshing force was integrated to obtain the virtual work possessed by the whole tooth. Also, the potential energy of the beam consists of five components corresponding to the five parts of stiffness. Subsequently, stiffness is calculated by corresponding potential energy. To simplify the calculation, a common method is to fit the meshing stiffness by Fourier series 20,21 where K am is the average meshing stiffness, which is always set as a suitable constant. ω t ϕ cos( + ) m m is the first-order harmonic component of the stiffness and stands for the variant with time, ϕ m is the initial phase.
Besides, the backlash in Equation (5) is always supposed as a fixed value or a steady random number that meets the normal distribution.
But in fact, the backlash is related to the machining error, and the machining process also determines the surface topography, hence the gear backlash must also be affected by the surface morphology.
Thus, the above methods of calculating backlash and stiffness cannot reflect the influence of microscopic features between meshing interfaces. Since the contact stiffness of two interfaces is under the influence of microscopic features, the method to assume stiffness as a constant will ignore the effects, resulting in reduced accuracy and limited application scope. In this study, we will propose a modified expression of meshing stiffness and backlash to reflect the influence of surface topography.

| Effects of surface morphology on the backlash
The fractal theory is a way to describe a rough surface with fractal characteristics, which obeys the cross-scale self-affinity and selfsimilarity of the surface in reality. Different from the statistics method, the fractal method has the advantage of being independent of instrument resolution. Meanwhile, fractal dimension D and fractal roughness G are the main parameters to determine , n denotes the frequency index that determines the length scale of asperities. Equation (7) cannot reflect the effects of fractal roughness G. To take G into consideration, Equation (7) can be modified as where L is the sampling length of the fractal surface. For meshing gears, L means actual contact length and can be derived as being the pressure angle.
Based on Equation (8), the backlash of gear pair considering surface topography is derived as where b 0 is the initial backlash as an appropriate constant, D G , Substituting Equation (10) into Equation (5), the modified backlash function is derived as Equation (11) indicates that the values of backlash are affected by surface topography directly. The details are illustrated in Figure 4.
In Figure 4A, the backlash of a single tooth, i.e., the height of asperities, decreases with fractal dimension D. Since a higher fractal dimension means a smoother surface, the height of asperities will be less under a big value of D. The result in Figure 4A is consistent with reality. Nevertheless, in Figure 4B, as the backlash is the value of

| Effects of surface morphology on meshing stiffness
Based on fractal theory, the stiffness of contact interfaces is the sum of all asperities at a microscopic scale. The stiffness of a single asperity in elastic deformation mode is where E is the equivalent elastic modulus, a is the real contact area of a single asperity.
For a single asperity, the elastic deformation mode occurs at initial contact, and then the transition into plastic deformation can take place. The contact area and critical deformation to divide the two stages are given as 24 respectively, where q is the coefficient of Poisson's ratio, σ y is the yield strength. The height of a single asperity is expressed as When the critical deformation exceeds the asperity height, the contact will be in pure elastic deformation mode. Combining Equations (13) and (14), the critical length scale l c can be derived as The distribution function of contact area between a pair of engaged gear surfaces is 25 where a L denotes the maximum contact area of a single asperity, a denotes the single asperity contact area, λ m is the coefficient of the curved interface, that is where D is equivalent fractal dimension of meshing surfaces decided by the fractal dimension of gear pair as G is the equivalent fractal roughness of meshing interfaces by a L denotes the largest contact area of a single asperity, that is 27 Substituting Equation (21) into Equation (18), the contact stiffness is finally expressed as It can be seen that in Equation (22), contact stiffness is decided by given surface topography. As the time variation of stiffness is caused by bending stiffness K b , shear stiffness K s , axial compression stiffness K a , and the fillet-foundation deformation stiffness K f , the contact stiffness can be applied as the average stiffness in Equation (6). And the modified meshing stiffness is Equation (23) shows that the time-varying stiffness is under the influence of surface topography. For various fractal dimensions and fractal roughness, i.e., the changeable morphology, the stiffness has obvious changes. The detailed change is displayed in Figure 5. As illustrated in Figure 5A

| SIMULATION AND DISCUSSION
To study the dynamic response under different surface topography, fractal parameters D and G are used to determine the microscopic morphology. The basic parameters of the gear pair for simulation are provided in Table 1. Meanwhile, the gear system is assumed to be subject to a light load to avoid the possible topography changes caused by overloading. The main dynamic parameters are set as follows, , as shown in Figure 6, the meshing interfaces are very rough, and the displacement varies from about 0.8 to 1.5. In the phase diagram, there is a close but wide band. The vibration frequency has only one wave crest in the FFT spectrum which is around 1.5. The points distribution appears to be relatively less concentrated but still not divergent. From above, the gear system is in quasiperiodic stage for rough surfaces at D = 1.1.
From Figure 7, as D p and D g grow to 1.3, i.e., the surface gets smoother, the gear system translates into an evident chaotic region. The wave shape has a high degree of similarity as G changes. (3) Another fractal factor G has more gentle effects on the dynamic response. The main waveform of the vibration remains the same, but the amplitude of the lower order oscillation decreases gradually when G increases. The range of velocity decreases as G grows. But the chaos characteristic has little change.
In this paper, the coupling of meshing stiffness and backlash under microscopic roughness is not considered, which can be researched in depth. Besides, the influences of morphology on other contact behavior like lubrication, friction, and stress are still worth to be studied. On the other hand, it will be a good research direction to apply the three-dimensional fractal method to gears.

CONFLICT OF INTEREST
The authors declare that there are no conflict of interest.

DATA AVAILABILITY STATEMENT
The datasets generated during the current study are available from the corresponding author on reasonable request.