Nonlinear response of gear system based on a novel backlash model with coupling dynamic change and surface topography

Backlash is one of main nonlinear internal excitation factors in gear transmission system and therefore has been widely concerned. Most existing models of backlash are based on a random constant, which ignore the dynamic characteristics of backlash itself and the effects of surface topography. To model the backlash precisely, in this paper, the constant part of backlash is revised through average height of all asperities in contact region related to surface roughness by fractal method. Simultaneously, the dynamic part is modeled considering the displacement of gear center motion that comes from shaft deformation in coupling dynamic meshing. A complete backlash model consisting of the two parts is established subsequently and a corresponding close-loop algorithm is proposed to solve system dynamics by coupling mesh stiffness and time varying pressure angle. Through time history charts, phase portraits and Poincare mapping as well as frequency spectrograms, calculation results clearly demonstrate the comprehensive effects of dynamic backlash on the nonlinear dynamics involving vibration amplitude, frequency and chaotic characteristics of a spur gear pair. The effects of surface topography on backlash and system nonlinear response including vibration amplitude and chaotic features are also analyzed, therefore dynamic backlash and surface topography are important factors that cannot be ignored in gear issue. The comparison with experimental data as well as other previous models is conducted to verify the superiority of proposed model.

to display the effects of surface morphology on backlash and compared it to normal constant backlash. Huang et, al. [15] analyzed the chaotic characteristics of gear system based on fractal backlash. Although the value methods of gear backlash through the above models are different, the models of backlash are all in a constant form.
Meanwhile, as the dynamic characteristics change in meshing process, the constant backlash is not accurate enough and dynamic model of backlash has been proposed. Chen et, al. [16] proposed a dynamic backlash model due to gear center motion with an invariant pressure angle. Yang et, al. [17] considered the backlash of planetary gears as an alternating meshing mechanism and represented backlash as dead zone of input torque. Li et, al. [18] established a dynamic backlash model with and without the effects of teeth wear respectively and analyzed the difference of system response. Liu et, al. [19] introduced a dynamic backlash model related to geometric eccentricity into lateral-torsional-rocking coupling gear model. Wang and Zhu [20] further applied the dynamic backlash model to planetary gear train system and helical angle was taken into account. Wang et, al. [21] proposed a simplified dynamic backlash model composed of two gears' eccentricity, which has a form of angular displacement. Yi et, al. [22] researched the dynamics of gear involving dynamic backlash under time-varying pressure angle. Shi et, al. [23] calculated comprehensive backlash of involute gear pair considering various kinds of deformation and contact ratio.
Although model of backlash has been investigated by different methods, a simultaneous consideration to the physical meaning of both invariant and variable parts of the backlash model is really rare. The absence of either morphology or dynamic backlash will lead to incomplete backlash modeling and limited nonlinear response accuracy. In this paper, the motivation is to establish a complete dynamic backlash model with regard to effects of gear center motion and microscopic morphology at tooth surface. The constant part of proposed model employs fractal method to obtain effective mean value of surface roughness and subsequently backlash is decided. This method can avoid the result that average height of all microscopic asperities always be zero no matter how rough the surface is due to the equivalent plane of fractal theory and truly reflect the influences of all asperities entering contact. The time varying part of backlash originates from gear center motion, which has inevitable influences on backlash. The relationship between gear center displacement and variation of backlash is derived geometrically. Corresponding to the characteristics of dynamic backlash, a close-loop algorithm is established to update the backlash value. The nonlinear response of system with proposed dynamic backlash model is solved and validation by experiment is made.

Modeling of dynamic gear backlash
In essence, as Fig. 1 shows, since there is always a difference between the actual and ideal tooth thickness along the pitch circle for various reasons, the backlash has inevitable existence. In a pair of meshing gears, the existence of the initial backlash will affect the operation of the system and in turn cause the backlash itself to change.
actual tooth thickness ideal tooth thickness As a result, backlash between meshing gears consists of constant and dynamic components. The constant part originates from initial invariant parameters such as manufacture or installation errors and is affected by surface topography. While the time varying part roots in the change of related parameters in actual dynamic meshing, among which the most important factor is the movement of gear center [16] [22]. The complete backlash can be expressed as Herein b s denotes the constant part, b dy (t) is the time-variant part.

Revised invariant part of backlash
Generally, the constant part of backlash is always set as a fixed value that confirms to standard normal distribution. The selection of the median value in a normal distribution is totally a random process, which leads to no physical basis for this behavior. In fact, from a microscopic point of view, the roughness of the gear surface is the factor that affects the backlash of the gear by changing the distance between the corresponding contact areas of the two teeth. As Fig. 2 illustrates, ideal backlash is b 0 between two perfectly flat surfaces. But in the actual meshing process, the meshing surfaces are rough and the profiles are decided by asperities of microscopic level, as Z p and Z g demonstrates. Hence from a microscopic point of view, it is the inevitable asperities on the surface of gears that are subjects in contact. This means that higher the asperities are, shorter the distance need the two surfaces cross to complete contact, thus reducing the backlash. Therefore, the actual value of backlash between rough surfaces will be decided by the average height of all asperity peaks of each surface i.e. the surface topography. Fractal theory, which has self-affinity and self-similarity across scales, is a widely used method contact region l Zp Zg b0 Fig. 2 Backlash under effects of surface roughness to describe surface morphology of machined interface. Based on fractal theory, the surface profile can be described by modified W-M function as [24] where l is the length scale of a single asperity, G denotes the fractal roughness, λ is the factor of surface frequency density, which is 1.5 for most surfaces. k is the level of asperity and k min is decided by λ k min = 1/l a , l a is the maximum base diameter of asperity. It can be seen that, the heights of asperities are decided by fractal dimension D and fractal roughness G. In practice, D stands for the amount of space taken up by profile waviness, G stands for height scale parameter from microscopic level, thus a smoother surface with fractal characteristics corresponds to a large D value and a small G value. Accordingly, to demonstrate the effects of surface morphology, Chen et, al. [14] has employed fractal theory to describe asperity heights as where t is sampling time [14], k is the level of single asperity. Although the value of backlash is influenced by surface topography in Eq. (3), the above fractal backlash model neglects the effects of fractal roughness G. At the same time, because only one asperity is calculated in the sampling time, the influence of the overall morphology of the contact area could not be adequately reflected. To improve the weakness, here we establish a new model of fractal backlash. As mentioned above, the backlash value depends on the average of all asperity heights. However, in fractal theory, the description of profile is based on an equivalent plane, which means simply taking the mean values of all asperities heights will always lead to a result of zero even the surface roughness is different completely, as displayed in Fig. 3. The red, green and blue lines represent rougher, medium and smoother surface topography respectively. In the meshing process, asperities of rougher surface have higher peaks and will enter contact early i.e. a small backlash. Similarly, smooth surfaces require more distance to touch at microscopic level i.e. a large backlash. Thus the value of backlash ought to be the difference between the initial clearance and the average height. But as the base line shows, all surfaces with various roughness in fact have a same base, which means take simply the arithmetic mean value as average asperity height will always be zero. Hence the effects of surface topography are hidden. To resolve this problem, a method that truly reflects the average height of all rough peaks in the contact area is needed to determine the backlash.
In this paper, we propose a method to calculate the arithmetic mean of the absolute values of all asperities and get effective mean value, which can be expressed as where L denotes the sampling length in fractal theory, while for meshing gears, L can be regarded as the length of contact region. By Eq. (4), b i represents the effective mean height of all asperities at meshing surface. Clearly, its value will change with fractal dimension D and fractal roughness G, which presents the influence of surface topography. Meanwhile, by assuming the initial backlash b 0 , the constant part of backlash is the difference of initial backlash b 0 and fractal backlash b i of each surface by Eq. (5) indicates that the constant part of backlash is determined by surface topography through fractal factors D and G. For a rougher surface, the asperity heights get higher but the backlash of gears decrease reversely. Besides, as a rougher surface means higher asperities peaks and therefore a higher average height, b i ought to have negative correlation with D.
By simulation, the actual situation of constant backlash values under various fractal factors is presented in Fig. 4. The effective average heights of asperities is represented by b i , while the value of constant backlash is defined by b s . It can be seen that, D has very strong effects. For rough surface with D = 1.2, because of high average asperity height, constant backlash presents a small value that is smaller than 20µm whatever G is. For very smooth meshing surface with D = 1.8, the backlash almost equals initial backlash and average asperity height gets very small. Effects of fractal roughness G has same trend as D that smoother surface with lower G leads to a bigger constant backlash. But the influences of fractal roughness G are relatively gentle compared to effects of D.

Time-variant part of backlash
In the meshing of gears, microscopic contact behaviors such as tooth deformation, thermal expansion or geometric eccentricity [23] cause the backlash to change i.e. the occurrence of dynamic backlash. Among the time-varying parameters that affect backlash, the most important one is the movement of gear center. When the shaft and bearing of gears are not considered as totally rigid bodies, the tiny deformation allows the gear to move in a small area i.e. has additional translational DOF (Degree of Freedom). Once the motion of gear center occurs, the meshing point will change and lead to a dynamic backlash. Fig. 5 demonstrates this phenomenon. The relationship between the change of backlash ∆ b and the variation of center distance ∆ a can be expressed as where ∆ a is the difference of distance between two centers due to motion, ∆ p is the distance in original tooth profile between ideal mesh point and the point corresponding to actual mesh point. The proof is as follow. It can be seen from Fig. 5 that, b 0 is the initial backlash and b is actual backlash under center motion. Since the movement of the two centers is relative, here we view one as fixed and the other moving for simplicity. Tooth profile drawn in dot line denotes the initial Along the actual backside-LoA, timely backlash b can be divided into two line segments, AB and BC. In reality, the variation of tooth thickness along the tooth profile is very small i.e. AE will be very fraction and AB will be close to DE in length. Since AB can be regarded as equal to b 0 , which means length of BC is the main reason of ∆ b and thus solving for ∆ b = b − b 0 means the same thing as determining the length of BC. For involute spur gear, the backside LoA is always orthogonal to the actual tooth profile. Also because actual tooth profile is parallel to initial tooth profile, arc BD is orthogonal to AC and right triangle BCD can be found. In triangle BCD, the length of BD (∆ p) can be approximated as arc length Herein R Pb denotes the radius of base circle of pinion, β is the rolling angle which equals the sum of evolving angle and pressure angle that decided by β = tan α. The subscript 0, c here is corresponding to initial point and actual point.
Subsequently, ∆ p i.e. length of side BD can be derived as ∆ p ≈ ⌢ BD. Back to right triangles BCD, where motion between points C and D mainly comes from overall movement caused by the shift of the center. Therefore, the length of the hypotenuse CD equals the motional distance of gear center i.e. ∆ a. Afterwards, the value can be obtained by the Pythagorean theorem and Eq. (6) can be proved.
The difference of gear center distance ∆ a is still an unknown parameter, thus the displacement of gear center need to be solved. To analysis the difference of gear center, as shown in Fig. 6, a 0 is original or ideal distance between two gears that equals R p + R g , where R p and R g is reference circle radius of pinion and gear respectively. a is the actual center distance after motion. At initial situation, as shown in Fig. 6, the gear centers are assumed to be on the x axis. After the translation motion occurs, the displacements of pinion and gear are x p , x g and y p , y g along x axis and y axis respectively. The real time distance a can be conveniently derived as Subsequently, the difference on distance of gear centers ∆ a is Substituting Eqs. (9) and (7) into Eq. (6) and finally the difference of backlash ∆ b is derived as Form Eq. (10), the displacement x i (i = p, g) and y i (i = p, g) is apparently time varying, which means that the center motion is time varying, and thus the difference of backlash ∆ b is also time varying. As mentioned earlier, the change in backlash is caused by the center motion of the gear, hence the dynamic part of the backlash b dy (t) corresponds to the time varying change ∆ b by b dy (t) = ∆ b. At this point, both parts of Eq. (1) are modeled, and a complete dynamic backlash model is obtained. The final expression is Eq. (11) illustrates that the proposed backlash model is a complete model contains both dynamic and constant part. The information of surface topography, gear center displacement and geometric profile of tooth are all embodied in the proposed model, which fits the physical reality better compared with the traditional backlash model. As displayed in Fig. 7, a dynamic model of gear system coupling DOF of gear centers translation due to flexibility of shafts is established involving dynamic backlash, time varying pressure angle and meshing stiffness.
Based on basic mechanical principles, the dynamic equations of gear pair can be expressed as where I i (i = p, g) is rotational inertia, T i (i = p, g) denotes external torque, m i (i = p, g) is the mass, c i j (i = p, g; j = x, y) and k i j (i = p, g; j = x, y) denote the damping and stiffness of pinion as well as gear along axis respcetively, α is the timely pressure angle, F m is the meshing force between gears, that is F m = c m δ + k m f (δ ) (13) where k m is the meshing stiffness, c m is the meshing damping factor determined by c m = 2ξ c (k m /m e ) 1/2 , ξ c is damping ratio and obeys 0.03 < ξ c < 0.17, m e is the equivalent mass, f (δ ) is nonlinear function of backlash in form of piecewise as [22] f = where B(t) is the complete dynamic backlash given in Eq. (11).
In the coupling dynamic model of gear pair, the dynamic transmission error (DTE) δ of system is expressed as [16][22] where e(t) is the static transmission error, which is always expressed as [9][25] e(t) = e r cos(ω h t + φ 0 ), where e r is amplitude, ω h is meshing frequency and φ 0 denotes initial phase.
Note: Here pressure angle α is not ideal meshing angle α 0 = 20. That is because of gear center motion, pressure angle also becomes dynamic, which is contradictory to constant pressure angle in most previous literatures. The actual pressure angle is the one between the actual meshing line and velocity of pitch point. With the effects of gear center motion, dynamic pressure angle can be expressed as [19] where a is the dynamic real-time gear center distance, R i (i = p, g) denotes radius of reference circle.
As for the meshing stiffness, its nonlinearity is usually demonstrated by fitting trigonometric functions [25], that is where k am denotes average meshing stiffness, ε is ratio of fluctuation to mean value, ω h is meshing frequency of gears and φ 0 is the initial phase. The value of k am can be equal to contact stiffness, which is always calculated based on Hertz theory that ignores the influences of surface topography. In our previous study [26], the fractal theory is employed to obtain the contact stiffness between curved surfaces with roughness, which can be applied here to calculate the meshing stiffness. The expression is where φ is material property by H/E, H is hardness, E is elasticity modulus, q is coefficient of Passion's ratio as q = 0.454 + 0.41v, λ c is the coefficient of curved surfaces expressed by (3FR/4E ) R and R is the equivalent curvature radius of gear pair by R = R p R g /(R p + R g ), a L is the largest contact area of a single asperity related to force applied on surface [26], l a is the largest asperity base diameter and l s is the smallest asperity diameter. Based on Eq. (18), taking k c as k am , the meshing stiffness can be obtained by Eq. (17).
After solving the pressure angle and meshing stiffness, the dynamic equations can be nondimensionalized and subsequently solved. The equivalent mass of meshing pairs is m e = (R 2 p /I p + R 2 g /I g ) −1 , and the nominal frequency is ω n = k m /m e , nominal dimension is b n . Accordingly, the dimensionless form of all the parameters can be derived as δ * = δ /b n , x * i = x i /b n (i = p, g), y * i = y i /b n (i = p, g), ζ = c m ω n m e , κ px = k px /m p ω 2 n , κ pm = k m /m p ω 2 n , ζ px = c px /ω n m p , ζ pm = c m /ω n m p , κ py = k py /m p ω 2 n , ζ py = c py /ω n m p , κ gx = k gx /m g ω 2 n , κ gm = k m /m g ω 2 n , ζ gx = c gx /ω n m g , ζ gm = c m /ω n m g , κ gy = k gy /m g ω 2 n , ζ gy = c gy /ω n m g , T * = T i /k m b n (i = p, g), T ah = e/b n . Then the nondimensionless dynamic equations are derived as 4 Nonlinear response of gear pair

Closed loop algorithm of solution
Different from traditional models with backlash of constant value or sinusoidal function, calculation process of proposed model needs to update the dynamic backlash values in each step of calculation responding to the obtained center displacement by past step. Both the time varying pressure angle and actual center displacement are updated and afterwards the actual backlash is calculated accordingly and substituted reversely into dynamic equation. Subsequently, next step of solving system dynamic equation can be continued. The whole algorithm forms a closed loop in which the center displacement of gear is determined by the initial backlash and the actual backlash in turn determines the center displacement. Also the surface topography factors will be taken into algorithm. The flowchart is demonstrated in Fig. 8.
By the above algorithm, the dynamic backlash values of system can be obtained. In Fig. 9, the time-domain curve of dynamic backlash changed with fractal dimension D under G = 1 × 10 −11 m is displayed. Obviously, the value of the dynamic backlash deviates significantly from the constant backlash and is always greater than the constant backlash. The reason is that the influence of gear center motion is considered, which makes the dynamic part larger. In detailed, when D = 1.2, average value of dynamic backlash locates at about 5.5 and constant backlash is only 4.5. With D = 1.4, the values of dynamic backlash are fluctuating around 7.2 while the constant backlash is only 4.9. While as D = 1.6, the average dynamic backlash is about 9.2 and the constant backlash is almost 5. Besides, the dynamic backlash has a quasi-periodic form due to periodicity of gear center motion whatever D is, but waveforms are various under different morphologies. The value of dynamic backlash itself also grows with fractal dimension D. In addition, these curves have relative various forms under different morphologies, indicating that constant part and dynamic part has coupling effects to each other and hence both are indispensable. By the comparison, the significant distinction of dynamic backlash and constant backlash proves the importance of dynamic characteristics in clearance modeling and the superiority of proposed model.  In Fig. 10, to illustrate the gear system response with dynamic backlash, the time history of DTE is presented under different surface topography and dimensionless frequency ω * = 0.2. In this paper, we assume that meshing surfaces of gear and pinion have same fractal factors, thus the fractal dimension and fractal roughness used to describe roughness hereinafter all mean D = D p = D g and G = G p = G g .
One of the most striking things is that, system DTE under dynamic backlash is always larger than DTE under constant backlash for D = 1.2, 1.4 and 1.6 and when D = 1.8, except for the wave trough, the value DTE with dynamic backlash is also larger. As Fig. 10(a) shows, with a rough surface that D = 1.2, the system response under dynamic and constant backlash has really small difference. The amplitude under dynamic backlash is slightly larger than the one under constant backlash with the same vibration form. The vibration both fluctuate around 70 and have close amplitude scale, meanwhile system under goes periodic motion. With the surface gets smoother, the situation has changed at D = 1.4, as displayed in Fig. 10(b). The most significant change is that system vibration decrease drastically to a mean value of about 22. Besides, the disparity becomes relative bigger and achieves about 3 at wave peak. But the form of wave still keeps a same form, which is very similar to situation at D = 1.2. When D = 1.6, as demonstrated in Fig. 10(c), great distinction appears between dynamic and constant backlash. The vibration form is totally different. The average value of dynamic backlash is about 12 while the one corresponding to constant backlash is only about 7. Besides, the curve of dynamic backlash is a four-periodic motion while constant model is quasi-periodic motion. Fig. 10(d) shows the time history of system with D = 1.8 i.e. very smooth. The system has different form under two models and even the chaotic characteristics change. The dynamic model becomes chaotic while the constant model still maintains periodic. Both models have same mean value at about 0 but the dynamic model has a larger scale of amplitude. Through above, the system response amplitude becomes different because of the dynamic backlash. The influences contain the form and the mean value of vibration as well as chaos. Meanwhile, an interesting phenomenon can be found that for rough surface with less initial backlash has conversely large amplitude. The amplitude of vibration gets about 115 for D = 1.2, while decrease to 30 at D = 1.4. The tendency shows monotonous and the max amplitude is only about 18 and 16 with D = 1.6 and D = 1.8 respectively. The reason maybe the compensation effects of dynamic clearance, which changes the amplitude. Besides, the meshing stiffness and damping are also affected by surface topography. For smooth surface with high stiffness, the contact property is enhanced and leads to small amplitude of DTE. From a practical point of view, the more precise gear with smoother surface should have smaller error, which is exactly consistent with our conclusion, indicating the correctness of proposed dynamic backlash model under the influence of topography.
Dynamic backlash also has effects on the vibration frequency of gear system. Corresponding to time history charts, the spectrum diagram derived from FFT (Fast Fourier Transformation) is shown in Fig. 11. At relatively rough situation D = 1.2 and D = 1.4, there is no clear change at the component of main frequency between dynamic and constant backlash, the only change is a small increase in the amplitude corresponding to the main frequency at ω * = 0.2 and ω * = 0.4. While for D = 1.6, as shown in Fig. 11(c), vibration frequency of gear system with dynamic backlash has more frequency components at about ω * = 0.59 and ω * = 0.75. Besides, the amplitude also grows by a greater amount. For D = 1.8, more and more complex frequency components appear under dynamic backlash compared to constant one. The amplitude has also become larger for different frequencies. The above phenomenon is consistent with time history that vibration amplitude grows under dynamic backlash and has more chaotic form. The variant may due to the introduction of dynamic backlash that forces the coupling between torsion frequency and gear center translation motion frequency. Meanwhile, all the situations of different roughness keep same that the value of main frequency has no change, which indicates that dynamic backlash has no effects on the frequency of main torsional vibration itself. The effects of surface topography on gear system response are demonstrated in Fig. 12 by the form of phase portraits and Poincare mapping. In Fig. 12(a), it can be seen that all the curves are thin closed lines and have completely different shapes, which indicates that gear system under all morphology undergoes periodic motion at given dimensionless frequency. But both the DTE and relative velocity shows decreasing trend with D. There is another phenomenon that curves of D = 1.2 and D = 1.4 have relative similar shape, but as D gets larger i.e. a smoother surface, the shape becomes completely different. This phenomenon declares that surface topography has greater influences at a smoother and more precise condition.
Comparing with fractal roughness G as a variable, the phase portraits is shown in Fig. 12(b). All the curves have same form even G has changed by an order of magnitude. But G still has effects, which are mainly reflected in the range of vibration amplitude. As G increases, the vibration amplitude also increases, in a minor degree. Given a larger G means higher roughness, this situation conforms to the reality.
From Poincare mapping, it can be seen that although system at different surface topography all presents periodic, but the points corresponding to D = 1.8 and G = 1 × 10 −12 m show better clustering. This also conforms to the rule that the smoother the tooth surface, the more precise the gear and more systematic the motion are.
From above, fractal dimension D has significant effects on gear dynamics while effects of G are more gentle. Therefore the consideration of morphology is essential for modeling of backlash.

Comparison and experimental validation
In this work, to verify proposed model and compare with available literatures, both the DF (Dynamic Factor) and RMS (Root Mean Square) of gear system are calculated and compared with experimental data.
DF means the ratio of maximum dynamic meshing force and statics root stresses [10] [25], that is where F m is meshing force from Eq. (13), F s is static force decided by average torque on the base radius [22].
The compared model is Chen's model based on fractal backlash from Ref. [25] and Kahraman's model that contains classic constant backlash from Ref. [10]. The experimental test rig and main parameters which can be found in Ref. [10] are presented in Tab. 1.
To make proposed model comparable, the effects of surface topography are eliminated to accord with the actual experimental parameters and the main distinction lays on the dynamic characteristics of backlash. The result is displayed in Fig. 13. All the models have same trend with experimental data. Proposed model has obvious advantage of approaching experimental data at relative low (0.2˜0.8) and high frequency (0.9˜1.2). Besides, the proposed model has the most approximation at the discontinue jump action at amplitude from both jump frequency and amplitude aspect. By contrast, predicted values of Chen' model and Kahraman's model both are much bigger than the experimental data. Since the difference comes from backlash model, the proposed dynamic backlash can be validated and its superiority is also demonstrated.
Herein A i denotes the amplitude of ith harmonic, n is the number of harmonics taken into account. The compared model is Cao's model that based on force-dependent meshing stiffness with constant backlash from Ref. [27] and Kahraman's experimental data from Ref. [28]. The test rig in Ref. [28] is a four square gear dynamic setup contains a precision gear pair. By accelerometers mounted tangentially to gear, the torsional vibration amplitude can be measured. Since Cao's model is based on stiffness and ignores the effects of backlash, the emphasis of comparison can focus more on the influence of the presence or absence of backlash. The experimental parameters are presented in Tab. 2 [28].
Overall, calculation values of proposed model have a better approach as illustrated in Fig. 14. Cao's model has an obvious tendency that is larger than the experimental data. While the proposed model shows a slight fluctuation around the experimental data. Meanwhile, the predicted values of proposed model are very close to the experiment data at low frequency (0.6˜0.8). By the results, the correctness and superiority of proposed model can be proved.

Conclusions
In this paper, a novel dynamic backlash model composed of constant part and time-varying part is established. The constant part is decided by surface topography through fractal method, which displays the effects of mean height of all asperities in contact. The time varying part comes from gear center motion at multi-DOF condition. The relation between dynamic backlash and gear center displacement is analyzed and the dynamic pressure angle is also considered. Incorporating meshing stiffness considering morphology, the system response is solved by a proposed close-loop algorithm. Subsequently, the effects of dynamic backlash and surface topography are analyzed. The values of RMS and DF calculated by proposed model and experimental data as well as other models are compared finally to verify proposed model. The main conclusions are as follow: (1)There shows great distinction between values of the dynamic form of backlash and constant one. The amplitude always gets larger and the quasi-periodic variation occurs. The dynamic backlash is also affected by surface topography through the values of constant backlash and the coupling effects of constant and dynamic backlash. By comparison, the proposed dynamic backlash has obvious difference with traditional constant backlash and its necessity is proved.
(2)Nonlinear response of gear transmission system is under comprehensive effects of proposed dynamic backlash. The influenced aspects include amplitude, frequency of DTE and relative velocity as well as chaotic characteristics. Dynamic backlash will enlarge the amplitude of vibration obviously. By dynamic backlash, the motion of gear center has a coupling effects on torsional vibration indicated through the change of frequency. For a smoother surface, the effects of dynamic backlash get relative stronger, in both frequency and amplitude.
(3)In terms of surface topography, the values of backlash in gear meshing are under direct influences of morphology. It also can be found that DTE of gear system has positive correlation with surface roughness. Therefore, improving surface smoothness will be helpful for realizing precise transmission of gears.
(4)From comparison with experimental data involving RMS and DF, the proposed model is verified. The calculation results of proposed model show good approximation to reality, proving the necessity of dynamic characteristics of backlash in modeling of gear meshing process.

Data availability statements
All data generated or analysed during this study are included in this published article.

Conflict of interest
The authors declare that they have no conflict of interest.

List of symbols
Alphabets a Actual gear center distance a 0 Initial gear center distance A i Amplitude of ith harmonic a L Largest contact area of single asperity B(t) The complete gear backlash b 0 Initial gear backlash b i Effective mean value of asperity height b n , x * , y * , T * , T ah Dimensionless dynamic parameters of gear system b s The constant part of gear backlash b dy (t) Average stiffness k c Contact stiffness k i j(i = p, g; j = x, y) Damping coefficient of bearing k m Meshing stiffness L Length of contact region l Length scale of asperity l a Largest base diameter of single asperity l s Smallest base diameter of single asperity m e Equivalent mass of gear pair m i (i = p, g) Mass m i (i = p, g) Stiffness coefficient m i j(i = p, g; j = x, y) Damping coefficient of bearing n The number of harmonics taken into account q Coefficient of Passion's ratio R i (i = p, g) Radius of reference circle R pb Radius of base circle of pinion T i (i = p, g) External torque v Passion's ratio z(l) Profile of surface morphology z(t) Fractal backlash Z p , Z g Profile of asperities Greek symbols α Actual pressure angle β Rolling angle ∆ a Change of gear center distance ∆ b Change of backlash Meshing frequency ω n Natural frequency of gear φ Material property ε Ratio of stiffness fluctuation ξ Damping ratio ξ , δ * , κ Dimensionless dynamic parameters of gear system Subscript 0 Initial c Actual g Gear p Pinion x x axis y y axis