Numerical Calculation Method of Meshing Stiffness for the Beveloid Gear considering the Effect of Surface Topography

The tooth surfaces of beveloid gears have diﬀerent topography features due to machining methods, manufacturing accuracies, and surface wear, which will aﬀect the contact state of the tooth surface, thereby aﬀecting time-varying mesh stiﬀness between mating gear pairs. Therefore, a slice grouping method was proposed in this paper on the basis of potential energy to calculate the total meshing stiﬀness of beveloid gears with the surface topography. The method in this paper was veriﬁed by ﬁnite element method (FEM). Compared with the calculation results of this paper, the relative error is 5.9%, which demonstrated the feasibility and accuracy of the method in this paper. Then, the inﬂuence of parameters such as pressure angle, helix angle, pitch angle, tooth width, fractal dimension, and fractal roughness on meshing stiﬀness was investigated, of which results show that pressure angle, pitch angle, tooth width, and fractal dimension have an incremental impact on the mean value of mesh stiﬀness. However, the ﬂuctuating value of mesh stiﬀness has also increased as the pressure angle, tooth width, and pitch cone angle increase. Both the helix angle and the fractal roughness have a depressive impact on the total stiﬀness. But the diﬀerence is that, with the increase of the helix angle, the ﬂuctuation of meshing stiﬀness has been decreased. Conversely, with the increase of the fractal roughness, the ﬂuctuation of meshing stiﬀness has been increased.


Introduction
Involute beveloid gears, which were first proposed by Mettitt in 1954, have the variable profile shift modification coefficient along axes direction. Beveloid gears can be used to realize the transmission form of parallel shafts, intersected shafts, and crossed shafts, which can be used in four-wheeldrive transfer cases and marine gearboxes with a down angle. However, due to beveloid gears having different addendum circles, dedendum circles, and pressure angles on different cross-sections, it is difficult to calculate the total mesh stiffness of beveloid gears. Furthermore, due to machining methods, manufacturing accuracies, and surface wear, the tooth surface of beveloid gears is always not smooth, and the surface topography features will affect the contact state of beveloid gears, thereby affecting time-varying mesh stiffness between mating gear pairs and changing the dynamic characteristics of transmission systems. erefore, research on the time-varying meshing parameters of the beveloid gear is beneficial to improve the meshing characteristics of the beveloid gear and enhance the transmission stability of beveloid gears.
To accurately and efficiently calculate the time-varying meshing stiffness of gears, domestic and foreign scholars have done a lot of research. e current methods for calculating the time-varying meshing stiffness of gears generally include the potential energy method and the finite element method (FEM) from the literature in recent years. For the FEM, Tang proposed the gear meshing stiffness calculation method based on the finite element numerical calculation and gave the relationship between the modification parameters and the mesh stiffness of the modification gear [1]. Fatih Karpat developed an improved numerical method for the mesh stiffness calculation of spur gears with asymmetric teeth based on the theory of finite element analysis [2]. Hao applied the loaded tooth contact analysis methods to obtain the time-varying meshing stiffness of planets gear and analyze the influence of time-varying meshing stiffness on the load-sharing characteristics of the aeronautical two-stage five-branching planets gear train [3]. Xue et al. used a two-dimensional finite element model to calculate the torsional stiffness of the external gear pairs and the internal gear pairs of the planetary gear set [4]. For the analytical method, a so-called potential energy method which was first presented by Yang and Lin in 1987 can calculate the mesh stiffness conveniently and effectively, which divided the potential energy between the meshing teeth into Hertz potential energy, bending potential energy, and axial compression potential energy [5]. Based on Yang's theory, the shear potential energy of teeth was considered and this method was applied to the investigation of gear dynamics [6]. Wang et al. proposed an improved timevarying mesh stiffness model of a helical gear pair, which considers the influence of the axial mesh force component [7]. Chen and Shao [8] and Chaari et al. [9] considered the energy effects of other gear teeth and calculated the timevarying meshing stiffness of normal gear teeth and cracked gear teeth. Based on the unloaded tooth contact analysis (TCA) and the nonlinear Hertzian contact theory, Li et al. [10] presented a new method to calculate the mesh stiffness considering the effects of tooth profile modification and eccentricity error. Mengjiao et al. [11] proposed an improved method for solving the time-varying meshing stiffness of helical gear pairs based on the slicing method. In addition, Han et al. also investigated the influence of helix angle, tooth width, and friction coefficient on the time-varying meshing stiffness [12]. Song et al. proposed the potential energy-based slice grouping method to calculate the mesh stiffness for straight beveloid gears with parallel axes [13]. In recent years, some scholars have carried out research on the contact stiffness of rough surfaces based on the fractal model proposed by Weierstras and Mandelbrot. Li and Zhu proposed a new method for predicting the tooth surface wearing of involute gears based on an actual tooth surface model and an improved fractal method [14]. Mo and Gong conducted the precise model of complex tooth surface microtopography and analyzed the influence of tooth surface friction and meshing frequency on the dynamic characteristics of the system [15]. Liu and Zhang combined the W-M fractal model with the FEM to calculate the local contact stiffness of involute spur gears [16]. e motivation of this paper is to provide a numerical calculation method of meshing stiffness for the beveloid gear with the surface topography. e current researches on the meshing stiffness of beveloid gears were based on the assumption that the tooth surface is absolutely smooth. But in fact, due to machining methods, manufacturing accuracies, and surface wear, the tooth surface of beveloid gears is always not smooth. In view of this, the slicing idea was used in this paper to discretize the thickness of the beveloid gear teeth and comprehensively analyzed the influence of the tooth surface friction and the tooth surface topography on the meshing stiffness throughout the meshing operation. e slice grouping method was proposed on the basis of potential energy to calculate the total meshing stiffness of beveloid gears with the rough surface topography. e sections of this paper are arranged as follows. Section 1 is a brief introduction. In Section 2, the calculation model of mesh stiffness for beveloid gears with the surface topography was established based on fractal theory. In Section 3, a mathematical model was derived to describe local topography features of actual tooth surfaces of beveloid gears, and the proposed method in this paper was verified by FEM. In Section 4, the influence of the geometry parameters and surface topography features of beveloid gears on the meshing stiffness was investigated. Section 5 shows some crucial conclusions.

Calculation Model of Mesh Stiffness for Beveloid Gears Based on Fractal Theory
In the traditional potential energy method, the total potential energy stored in the mesh gear system was assumed to include five components: Hertzian energy U h , bending energy U b , shear energy U s , axial compressive energy U a , and base body potential energy U f . ey can be used to calculate Hertzian mesh stiffness kh, bending mesh stiffness k b , shear mesh stiffness k s , axial compressive stiffness k a , and filletfoundation stiffness k f , respectively. rough the knowledge of elastic mechanics, the five parts are separately calculated for the stiffness components. And then, the five stiffness components are combined to obtain the total mesh stiffness. However, because of the special tapered tooth shape of beveloid gears, the traditional potential energy method cannot calculate the mesh stiffness directly. Besides, the conventional potential energy method does not consider the influence of tooth surface friction and tooth surface topography. In summary, an improved potential energy method based on fractal theory was proposed to calculate the total meshing stiffness for beveloid gears of a rough surface in this paper. Because the axis-direction modification coefficient of beveloid gears changes linearly, the beveloid gear can be regarded as a continuous superposition of a certain number of slices with the same thickness and different modification coefficients. As shown in Figure 1, each piece of beveloid gears is regarded as a modified cylindrical spur gear. e mesh stiffness for each slice can be calculated using the potential energy method. Finally, the mesh stiffness for all the slices along the axial direction is combined to obtain the total mesh stiffness of the beveloid gear pairs.
To accurately calculate each slice's meshing stiffness, the teeth of beveloid gears are simplified into a variable crosssectional cantilever beam with variable linear along the axial direction. In Figure 2, the micrometeor area is A i , and the cross-sectional modulus is I i , the distance between any meshing point j. e gear origin O is x b . r b is the base circle radius. r a and r f are the addendum circle radius and the dedendum circle radius. y b is the half tooth thickness at any meshing point j. F n is the normal meshing force on the tooth surface. F x and F y are the partial force of the normal meshing force F n along the directions x and y, respectively.

Mathematical Problems in Engineering
Based on the potential energy method and the slice grouping method, the five components of the potential energy of each slice under the action of normal meshing force F i n can be expressed as where the contact stiffness is k i n , the bending stiffness is k i b , the shear stiffness is k i s , the axial compressive stiffness is k i a , and the fillet-foundation stiffness is k i f , respectively.
According to the cantilever beam model and the knowledge of elastic mechanics, the stored bending potential energy, shear potential energy, and compressive potential energy of each slice are expressed as where the component force in the x-axis direction can be expressed as e component force in the y-axis direction can be expressed as Equivalent bending moment can be expressed as x b x c

Mathematical Problems in Engineering
where μ(t) is the friction coefficient, α k is the working pressure angle at any meshing point k, and α 0 is the pressure angle at the pitch point.
Due to the relative sliding speed, the rolling ratio and other parameters are constantly changing. e tooth surface friction coefficient has changed with time, and the direction of tooth surface friction has changed at the pitch point. Xu and Kahraman [17] used the EHL model to conduct a large number of simulation tests on the friction coefficient between the gear teeth, and the tooth surface friction coefficient μ(t) can be expressed as where v 0 is the absolute viscosity and V ek is the entrainment velocity. S is the tooth surface roughness. b i (i � 1∼9) is the regression coefficient. SR k is the rolling ratio. P hk is the tooth surface contact pressure, and R k is the combined radius of curvature at any meshing point k. e calculation parameters are detailed in the literature [17]. Moreover, due to the fact that the effective contact part of the beveloid gear teeth can be regarded as the trapezoidal section in the axial direction, the calculation formula of I i , A i , and G can be expressed as where S ib and S is are the tooth thickness of the large section and the small section of beveloid gears. ] is Poisson's ratio. E e is the equivalent modulus of elasticity. B is the tooth width. e bending stiffness k i b of each slice can be expressed as erefore, the bending stiffness of beveloid gears can be expressed as Similarly, the shear stiffness k s and compressive stiffness k a of beveloid gears can be expressed as e fillet-foundation stiffness k f of beveloid gears can be calculated by the following equation: where μ f and S f are shown in Figure 2, L * , M * , P * , and Q * are constants related to R a , R f , and θ f . e calculation method is detailed in the literature [11].
Generally, Hertz's theory holds that the tooth surfaces in contact with each other are frictionless, and the contact deformation of the meshing area has a linear relationship with the normal meshing force. Due to machine body error, installation error, deformation error, and other factors, the tooth surface is generally uneven and rugged, and the machined surface can be characterized by continuity, nondifferentiability, and statistically self-affinity. As shown in Figure 3, the tooth surface comprises a series of randomly distributed asperities, and the actual contact area is only a tiny part of the theoretical contact area. And the asperities deform elastically, elastic-plastically, or plastically when two rough tooth surfaces contact each other.
Miao and Huang [18] proposed the expression equation of the total contact stiffness based on the fractal model and suggested critical indices to distinguish the deformation forms at each length scale. e elastic critical index can be defined as where int[] is the integer part of the value in the parentheses, K is the hardness coefficient related to the Poisson ratio of the material, H is the hardness of the material, and E is the elastic modulus. D is the fractal dimension, G is the fractal roughness, c is the spatial frequency of the rough surface, and c is generally taken as 1.5. Similarly, the plastic critical index can be defined as e asperities are in elastic deformation when n min < n < n ec , and the contact stiffness K nec in the elastic regime can be expressed by integrating into the whole contact surface as K nec � n�n ec n�n min a n a n+1 ′ k n n a r da r e asperities are in elastic-plastic deformation when n ec < n < n pc , and the contact stiffness K npc in the elasticplastic regimes can be expressed by integrating into the whole contact surface as K npc � n�n pc n�n ec a n a n+1 ′ k n n a r da r e asperities are in plastic deformation when n pc < n < n max , and the contact stiffness in the plastic regimes should take 0.
It follows that the total contact stiffness based on the fractal model can be given by erefore, the total mesh stiffness of the beveloid gear can be given by combining the five stiffness components as

e Mathematical Model of Beveloid Gears with the Surface Topography.
According to the space meshing principles and the processing method of beveloid gears, the tooth face equations of working tooth surface and the tooth root surface of beveloid gears are derived utilizing the virtual rack cutter's tooth face equations. As shown in Figure 4, the normal section of the virtual rack cutter consists mainly of two parts, including the fillet curves BC and the straight edges CD that generate the working tooth surface and the tooth root surface of the beveloid gear respectively.
(1) In the coordinate system S n , the tooth face equation of the straight edges CD of the virtual rack cutter can be expressed as Mathematical Problems in Engineering (2) In the coordinate system S n , the tooth face equation of the fillet curves BC of the virtual rack cutter can be expressed as follows: (22) Figure 5 shows a schematic of the coordinate relationship between the normal section of the virtual rack cutter and the generated beveloid gear. Based on the space meshing principles, the rack cutter's tooth face equations can be expressed by using a series of transformation matrices given by S n ⟶ S p ⟶ S c . en, the virtual rack cutter surfaces can be expressed in coordinate system S c as follows: e tooth face equations of the rack cutter in the coordinate system S 1 can be derived from the tooth surface equation of the rack cutter in coordinate system S c by using a series of transformation matrices given by S c ⟶ S 0 ⟶ S 1 . e tooth face equations can be expressed in coordinate system S 1 as follows: Mathematical Problems in Engineering rough the relative movement of the rack and pinion, the tooth surface equation of the beveloid gear can be derived from the tooth surface equation of the rack cutter in the coordinate system S 1 . And the beveloid gears working tooth surfaces can be expressed as follows: where x c , y c , and z c are the position vector of the working surface of the rack in the coordinate system S c . n → cx and n → cy are the unit normal vector of the working surface of the rack in the coordinate system S c . e rough tooth surface topography can be considered a series of randomly distributed asperities superimposed on the theoretical tooth surface's normal direction. erefore, the mathematical model of beveloid gears considering the surface topography characteristics can be expressed as e height of the asperities on the rough tooth surface can be given as where D x and D y are the fractal dimension of the rough surface in the directions x and y, G x and G y are the fractal roughness of the rough surface in the directions x and y, respectively, and n 1 and n 2 are the number of sampling points within a finite length of the rough surface in the directions x and y, respectively.
As shown in Figures 6(a)-6(e), the rough surface topography characteristics under different fractal parameters based on fractal theory were simulated in this paper. It can be seen obviously that as the value of the fractal dimension D increases, the simulated rough surface topography gradually becomes more complicated, and the fractal dimension is D positively correlated with the complexity of the rough surface. e fractal roughness G does not affect the complexity of the rough surface. Still, as the value of the fractal roughness G increases, the amplitude of the simulated rough surface topography gradually increases, and the fractal roughness G is positively correlated with the flatness of the rough surface.

Rough Surface Topography Measurement and Fractal Parameter Calculation
Material and Specimens. Geometric parameters of beveloid gear specimens are listed in Table 1. As shown in Figure 7(b), tooth surfaces of beveloid gear specimens present prominent topography characteristics after the long-term meshing operation. Tooth surface measurements of the specimens were performed by a coordinate measuring machine to accurately describe the tooth surfaces' topography characteristics, as shown in Figure 7(a). Measurement results of tooth surfaces of the specimens are shown in Figure 7(c). Extract the roughness value of the tooth width direction and the tooth profile direction from the measurement results, respectively, and the results are shown in Figure 8. It can be seen that the surfaces of beveloid gear specimens present prominent topography characteristics, which are a series of ravines and ridges in Figure 7(c). However, it needs to be emphasized that the fluctuation amplitude of the interface contour curve in the tooth profile direction is greater than the tooth width direction, which means that the interface contours in the tooth width direction are more smooth than the tooth profile direction.
According to Qi et al. [19], the fractal parameters could be obtained based on the structure-function method (SFM), the SFM considers the interface contour curve as a time series (x), and then the time series with fractal characteristic satisfies where t denotes the interval between data points, S(t) is a function of t, x is the abscissa on the contour curve, Z(x) is the contour height corresponding to the coordinate x, and [Z(x + t) − Z(t)] 2 represents the arithmetic mean of the difference. e double logarithmic plot of log (S) and log (t) is shown in Figure 9. en the fractal dimension D and fractal roughness G can be calculated by equations (34) and (35), and the calculation result is shown in Table 2.
where C s is scale coefficient (0 < C s < 1). α is the slope of the fitting line on the double logarithmic plot.

e Solid Model of Beveloid Gears with the Surface
Topography. e tooth surfaces of beveloid gears have different topography characteristics due to different machining methods, manufacturing accuracies, and surface wear, and they are always not smooth. And the rough tooth surface topography characteristics can be considered a series of randomly distributed asperities superimposed on the theoretical tooth surface's normal direction. Based on the fractal parameters in Table 3, the asperities' height on tooth surface with the rough topography can be calculated through equation (33). en the coordinates of the points on the rough tooth surface of beveloid gear could be obtained by superimposing the height of the asperities to the theoretical tooth surface's normal direction based on equation (31), shown in Figure 10(c). As shown in Figures 10(a) and 10(b), the solid models of the beveloid gear with the anisotropic surface topography were conducted by Solidworks software after surface reshaping of the spatial point sets.

Comparison and Verification of Meshing
Stiffness. e FEM was used to verify the mesh stiffness calculation model proposed in this paper. e main geometry parameters of the beveloid gear pair are shown in Table 3. e above 3D solid models were imported into the finite element analysis software ANSYS to conduct the mesh model for the beveloid   gear. To ensure the accuracy of the calculation, we refined the mesh of the tooth surface with topography characteristics, as shown in Figure 11. e material is defined with Young's modulus 209 GPa and Poisson's ratio 0.26. e rotating speed is 500 Rpm, and the loaded torque is 600 Nm. e type of element in the finite element analysis is Tetrahedrons, and the number of nodes is 1253538. e calculation flowchart diagram for the total mesh stiffness of the beveloid gear is shown in Figure 12. Figure 13 shows the distribution of the contact stress on the rough tooth surface, and it is not difficult to find that the rough surface contact area can be seen as a set of discrete contact areas. en, the mesh stiffness can be solved by [20] K � F n λe l , where F is the total contact force between the gear teeth, λ is the equivalent radius, λ � M/F, M is the total torque in the direction of axis, and e l is the transmission error. Figure 14 shows the time-varying mesh stiffness of the beveloid gear pairs with the surface topography based on the potential energy method and FEM method. e maximum value of single tooth meshing stiffness is 14.63 N/mm. e maximum value of single tooth meshing stiffness is 14.63 N/ mm, and the minimum value is 13.92 N/mm. e relative error is 5.1%. e maximum value of the total mesh stiffness is 26.89 N/mm, and the minimum value is 25.38 N/mm. e relative error is 5.9%.

The Influencing Parametric Analysis on the Time-Varying Meshing Stiffness
Time-varying mesh stiffness is one of the main reasons leading to the vibration in a gear transmission system. Analyzing the influencing factors of beveloid gears mesh stiffness has great significance for improving its meshing characteristics and transmission stability. Compared with traditional involute cylindrical gears, beveloid gears have different addendum circle, dedendum circle, and pressure    angle on different cross-sections, and the meshing characteristics are more complicated. e geometry parameters of beveloid gears mainly include pressure angle, helix angle, pitch angle, and tooth width. We define fluctuating value as the max-min of time-varying mesh stiffness divided by the mean value to evaluate the fluctuation degree. In this paper, the influence of the geometry parameters and surface topography of beveloid gears on the meshing stiffness was investigated by changing the geometry parameters shown in Table 3.    Figure 15(a) shows the impact of normal pressure angle on the total mesh stiffness of beveloid gear. It can be found that the total meshing stiffness has increased with the adding of the pressure angle, and the total mesh stiffness improved by 6% when the pressure angle increased from 17°to 21°. e reasonable explanation of this is that the tooth thickness at the root of the tooth increases.

Pressure Angle.
e radius of curvature of the tooth surface increases as the pressure angle increases. erefore, the total mesh stiffness of the beveloid gear has increased. However, as shown in Figure 15(b), the fluctuating value has increased when the mean value of mesh stiffness increases with the pressure angle increase. It may affect the smoothness of gear transmission and may cause vibration and noise problems.

Helix Angle.
e impacts of different helix angles on the total mesh stiffness are shown in Figure 16 Define dl and slice the tooth (i = 1... n) Total mesh stiffness:  angle, and the mean value of mesh stiffness decreased by 7% when the helix angle varies from 0°to 12°. It can be observed in Figure 16(b) that the single tooth contact zone increases and the double teeth contact zone decrease, and the fluctuating value of mesh stiffness has a noticeable reduction with the increasing helix angle. erefore, it can be concluded that a larger helix angle can effectively improve the stability of the beveloid gear transmission under the condition of meeting actual working torque.

Pitch Angle.
e effects of different pitch angles on the time-varying mesh stiffness of beveloid gears are shown in Figure 17(a). It is not difficult to find that the pitch angle has little effect on the mesh stiffness of beveloid gears, the pitch angle increases by 1°each, and the mean value increases by approximately 1.6%. However, as shown in Figure 17(b), the fluctuating value has increased with the adding of the pitch angle, which indicates that a larger pitch angle may cause the vibration and noise problems of the beveloid gear transmission.

Tooth Width.
e effects of different tooth widths on the time-varying mesh stiffness of beveloid gears are shown in Figure 18(a). It can be found that the total mesh stiffness increased in step with the adding of the tooth width, and the total mesh stiffness of beveloid gears improved by 6% when the tooth width increased from 22 mm to 30 mm. e reasonable explanation is that the length for the contact line is increased, which means that the adding of the tooth width will reduce the unit load under constant torque, so it contributes to enhancing the capability of tooth resisting deform. However, it can be seen from Figure 18(b) that the fluctuating value of mesh stiffness has an obvious increase with the adding of mean value, which means that the adding of tooth width may increase the vibration of beveloid gear transmission.

Fractal Dimension.
According to Mandelbrot's fractal theory, the value of fractal dimension D is between 1 and 2, and the number 2 means that the machined surface is absolutely smooth. Figure 19(a) shows the total mesh stiffness of beveloid gears when the fractal dimension D varies from 1.4 to 1.8. And it is not difficult to find from Figure 19(b) that the total mesh stiffness has increased with the adding of the fractal dimension D, and the mean value of mesh stiffness increases by about 1.4% with each increment 0.1 for the fractal dimension D. It needs to be emphasized that the  fluctuating value of mesh stiffness has an obvious decrease with the adding of mean value, which indicates that more smooth surface may decrease the vibration of gear transmission. Figure 20(a) shows the total meshing stiffness of beveloid gears with different fractal roughness G. Converse to the fractal dimension, the total mesh stiffness of beveloid gears has decreased with the adding of the fractal roughness G, and the mean value of mesh stiffness decreased by 8.6% when the fractal roughness varies from 2.0E−6 to 1.0E−5. And it can be seen from Figure 20(b) that the fluctuating value of mesh stiffness has an obvious increase with the adding of fractal roughness, which indicates that the roughness of the tooth surface is one of the causes of vibration and noise.

Conclusion
is paper concerns the numerical calculation method of meshing stiffness for the beveloid gear with the surface topography, and the influence of the geometry parameters and surface topography parameters of beveloid gears on the meshing stiffness was investigated. e following conclusions are drawn: (1) In this study, an improved potential energy method based on the slice grouping method was presented to calculate the mesh stiffness for beveloid gears with the tooth surface topography. e proposed method in this paper was verified by FEM. Compared with the calculation results of this paper, the relative error is 5.9%, which demonstrated the feasibility and accuracy of the method in this paper. (2) Based on fractal theory and the space meshing principles, a mathematical model was derived to describe local topography features of actual tooth surfaces of beveloid gears, and the fractal parameters were calculated by the SFM. e analytical results show that fractal dimension D has an obvious positive impact. e larger the fractal dimension, the smoother the tooth surface. Conversely, the larger the fractal roughness, the rougher the tooth surface.
(3) e analytical results show that pressure angle, pitch angle, tooth width, and fractal dimension have an incremental impact on the total stiffness of beveloid gears. However, the fluctuating value of mesh stiffness has also increased as the pressure angle, tooth width, and pitch cone angle increase. It is worth emphasizing that, with the increase of the fractal dimension, the mean value of meshing stiffness has been improved, but the fluctuation of meshing stiffness has been decreased. (4) Both the helix angle and the fractal roughness have a depressive impact on the total stiffness of beveloid gears. But the difference is that, with the increase of the helix angle, the fluctuation of meshing stiffness has been decreased, which means that the increase of the helix angle can improve the gear transmission's stability. Conversely, with the increase of the fractal roughness, the fluctuation of meshing stiffness has been increased, which indicates that the roughness of the tooth surface is one of the causes of vibration and noise.
Data Availability e raw/processed data required to reproduce these findings cannot be shared at this time as the data also form part of an ongoing study.