{'fraction_non_alphanumeric': 0.024012886791039186, 'fraction_numerical': 0.00880347643665243, 'mean_word_length': 4.702841273232215, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Solution of Spur Gear Meshing Stiffness and Analysis of Degradation Characteristics

Excitation of time-varying meshing stiffness of gears is one of the main sources of vibration of gear system. The study of time-varying meshing stiffness has always been an important branch of gear dynamics. Distribution of loads on gear teeth when the gears are engaged is affected by many factors: method of manufacture, precision of manufacture, error of installation, modification of tooth profile, heat treatment of tooth surface. So it is difficult to calculate the time-varying meshing stiffness accurately. At present, the calculation of meshing stiffness is mainly divided into two methods: method of theoretical numerical calculation method and FEM. Many scholars have studied the solution of meshing stiffness of gears. Bu [1] used an improved method based on the linear programming method to calculate the meshing stiffness and distribution of loads, and calculated the meshing stiffness of a pair of inner and outer meshed helical gears at different helix angles. It is feasible and accurate to use the linear programming method to solve the meshing stiffness. But a slight disadvantage is very complicated. Tang [2] used the FEM to calculate the meshing stiffness of surface gears, and calculated stiffness of the spur gear and compared it with the existing literature. He verified the correctness of this method, but did not further prove it with a test-bed. Chang [3] proposed a method to calculate the gear meshing stiffness by combining the FEM with the elastic contact theory. This method used the FEM to predict the overall deformation of the object. At the same time, the elastic contact theory can be used to calculate the local contact deformation accurately. The insufficiency is not related to experimental verification. Based on FEM, Zhang [4] built a contact model of spur gear to calculate contact stress by simulation, which results were verified by means of calculation of Hertz equation of contact theory. Transmission of spur gear system was analysed and time-varing mesh stiffness using Ishikawa equation was calculated by Matlab software. Its disadvantage is that no mechanical test-bed is used to further prove the relevant conclusions. The meshing stiffness synthesis method of multi-point meshing condition was established by Zhang Jianyu [5], based on the inner kinetic relationship of the gearbox and phase relation between the planets. For the case of sun gear and ring with crack of different depths, the meshing stiffness distribution of cracked gear meshing with planetary gear was discussed, and finally the stiffness synthesis of the internal and external meshing of the same planet gear was completed. His results provide some insights for the study of crack stiffness. Lewicki [6] established the guidelines for spur gear design based on linear elastic fracture mechanics using finite element method for variety of gear tooth and rim configurations, to predict the crack propagation paths. In another study Lewicki and Ballarini [7] studied, the effect of gear rim thickness on the crack propagation direction, numerically and experimentally. Both the above studies use two dimensional finite element code FRANC for the simulation of the crack propagation, numerically and suggested that the crack propagation paths are continuous, smooth and straight with slight curvature. Various finite element models and analytical methods have been developed for calculation of gear mesh stiffness [8]. Zouari [9] used 3-D FEM to see the effects of the depth, the orientation of propagation and of the position on blending radius of the crack on the instantaneous mesh stiffness, consecutively, a 2-D FEM is used to analyze the variation of stress intensity factor on the tooth base according to crack depth, crack propagation angle, and the crack position. In this study, the finite element analysis software is used to solve the non-linear contact problem of gears. The hexahedron element can be used to realize mapping grid partition. Combining with the refined modeling method for local cells in contact bands, the contradiction between the finite element contact analysis accuracy and the calculation efficiency of gear can be solved effectively. The meshing stiffness of the gear is solved by the FEM, and the theoretical calculation is compared with the finite element calculation result. Finally, it is proved that the FEM is reliable for the meshing stiffness. The FEM can be applied to the dynamic analysis of the gear system. At the same time, the dynamic model of the transmission system of gear is established and the response of vibration of the dynamic simulation is obtained. The experimental results are used to verify the reliability of the simulation and the accuracy of calculation of the finite element meshing stiffness is illustrated.


Introduction
Excitation of time-varying meshing stiffness of gears is one of the main sources of vibration of gear system. The study of time-varying meshing stiffness has always been an important branch of gear dynamics. Distribution of loads on gear teeth when the gears are engaged is affected by many factors: method of manufacture, precision of manufacture, error of installation, modification of tooth profile, heat treatment of tooth surface. So it is difficult to calculate the time-varying meshing stiffness accurately. At present, the calculation of meshing stiffness is mainly divided into two methods: method of theoretical numerical calculation method and FEM. Many scholars have studied the solution of meshing stiffness of gears. Bu [1] used an improved method based on the linear programming method to calculate the meshing stiffness and distribution of loads, and calculated the meshing stiffness of a pair of inner and outer meshed helical gears at different helix angles. It is feasible and accurate to use the linear programming method to solve the meshing stiffness. But a slight disadvantage is very complicated. Tang [2] used the FEM to calculate the meshing stiffness of surface gears, and calculated stiffness of the spur gear and compared it with the existing literature. He verified the correctness of this method, but did not further prove it with a test-bed. Chang [3] proposed a method to calculate the gear meshing stiffness by combining the FEM with the elastic contact theory. This method used the FEM to predict the overall deformation of the object. At the same time, the elastic contact theory can be used to calculate the local contact deformation accurately. The insufficiency is not related to experimental verification. Based on FEM, Zhang [4] built a contact model of spur gear to calculate contact stress by simulation, which results were verified by means of calculation of Hertz equation of contact theory. Transmission of spur gear system was analysed and time-varing mesh stiffness using Ishikawa equation was calculated by Matlab software. Its disadvantage is that no mechanical test-bed is used to further prove the relevant conclusions. The meshing stiffness synthesis method of multi-point meshing condition was established by Zhang Jianyu [5], based on the inner kinetic relationship of the gearbox and phase relation between the planets. For the case of sun gear and ring with crack of different depths, the meshing stiffness distribution of cracked gear meshing with planetary gear was discussed, and finally the stiffness synthesis of the internal and external meshing of the same planet gear was completed. His results provide some insights for the study of crack stiffness.
Lewicki [6] established the guidelines for spur gear design based on linear elastic fracture mechanics using finite element method for variety of gear tooth and rim configurations, to predict the crack propagation paths. In another study Lewicki and Ballarini [7] studied, the effect of gear rim thickness on the crack propagation direction, numerically and experimentally. Both the above studies use two dimensional finite element code FRANC for the simulation of the crack propagation, numerically and suggested that the crack propagation paths are continuous, smooth and straight with slight curvature. Various finite element models and analytical methods have been developed for calculation of gear mesh stiffness [8]. Zouari [9] used 3-D FEM to see the effects of the depth, the orientation of propagation and of the position on blending radius of the crack on the instantaneous mesh stiffness, consecutively, a 2-D FEM is used to analyze the variation of stress intensity factor on the tooth base according to crack depth, crack propagation angle, and the crack position.
In this study, the finite element analysis software is used to solve the non-linear contact problem of gears. The hexahedron element can be used to realize mapping grid partition. Combining with the refined modeling method for local cells in contact bands, the contradiction between the finite element contact analysis accuracy and the calculation efficiency of gear can be solved effectively. The meshing stiffness of the gear is solved by the FEM, and the theoretical calculation is compared with the finite element calculation result. Finally, it is proved that the FEM is reliable for the meshing stiffness. The FEM can be applied to the dynamic analysis of the gear system. At the same time, the dynamic model of the transmission system of gear is established and the response of vibration of the dynamic simulation is obtained. The experimental results are used to verify the reliability of the simulation and the accuracy of calculation of the finite element meshing stiffness is illustrated.

Calculation of meshing stiffness of variable cross section based on cantilever beam
The gear teeth are equivalent to cantilever beams, as shown in Fig. 1. The gear tooth is a cantilever beam with one side fixed completely, and the intersection of the force F and the center line of gear is A. The middle of AO can be seen as many tiny cross-section rectangles as Si in Fig. 1. The expression of the gear of the meshing stiffness is as Eq.
where: n F is normal meshing force of tooth surface, N; n  is comprehensive deformation of gear teeth.
The comprehensive deformation of gear teeth includes the elastic deformation of contact caused by the local Hertz contact, the deformation of the gear teeth caused by bending of the gear teeth, the deformation of the wheel, the deformation of the bearing, the deformation of the shaft and the supporting structure, et al [10][11]. In order to simplify the calculation, the bending deformation, the deformation of the wheel and the elastic deformation of Hertz contact in the process of meshing are considered without considering the shaft and the bearing [12][13].
The bending deformation of the gear can be obtained by Eq. (2), and then the Eq. (1) can be used to solve the bending stiffness of the gear teeth.  v is Poisson's ratio; Sh is shear coefficient, mm. When the cross section is rectangular, its value is 1.2; Ai is area of rectangular section, mm 2 ; Ii is moment of inertia of rectangular section, mm; The gears will also be deformed in the process of meshing. The size of deformation can be obtained by using Eq. (3). Then, the gear stiffness is obtained according to the equation of the gear meshing stiffness: where: W is tooth width of gear, mm; Sf is root tooth thickness, mm; uf is distance between the line of force and the center line of the gear to the root circle, mm; L, M, P, Q are the values of polynomials of hf and θf can be obtained by Eq. (4) and Table 1.
where: hfi is ratio of the root circle radius to the hub circle radius, mm; θf is half of the central angle of the root thickness, mm; A, B, C, D, E and F are calculated constants. The values are shown in Table 1. The gear will be elastically deformed due to contact during meshing, and the stiffness of hertz contact of the gear formed by the elastic deformation can be obtained from Eq. (5).
Considering the three kinds of deformed stiffness of the main and driven wheels, the composite algorithm of gear meshing stiffness can be obtained, as shown in Eq. (6).
where: kb1, kb2 are bending stiffness of main and driven wheels, N/m; kf1, kf2 are deformed stiffness of main and driven wheels, N/m; kh is stiffness of Hertz contact, mm. The number of teeth participating in the meshing of the gears changes periodically with time. This causes the gear meshing stiffness to also exhibit periodic changes. The gear time-varying meshing stiffness is expressed as a Fourier series form with respect to the driving gear rotational angular displacement θ1 (°): where: k0 is average value of meshing stiffness; kj is coefficient of Fourier series; φj is phase angle; T=360/Z1, Z1 is the number of active gear teeth. The Fourier series of time-varying meshing stiffness with respect to time obtained by transformation is: where: ωe=2πZ1n1/60; n1 is number of driving gears, r·min -1 .

Calculation of Hertz theory of gear contact and analysis of contact stress based on FEM
The three-dimensional model of the gear is established by SolidWorks, and the gear parameters are shown in Table 2. The hexahedral element can be used to divide the gear mapping grid. In order to ensure the accurate contact and solution of tooth surface, the mesh of contact area of tooth surface is refined, and setting the mesh to 1.2 mm. The finite element model is obtained as shown in Fig. 2. The material properties of the finite element model are given under the analysis of the transient dynamics module by workbench, as shown in Table 3. Adding a pair of degrees of freedom to the main and driven wheels. The rotational speed of the main wheel is 110 rad/s, and the same rotational direction of 1050 Nm torque is added to the driven wheel. Considering the nonlinear contact behavior of tooth surface meshing, the type of tooth surface meshing is selected as frictional contact. The contact algorithm is chosen as augmented Lagrange method.
The overall equivalent stress results of the gear obtained are shown in Fig. 3.   Fig. 4, we conclude that the results of the two algorithms are very close, and there is a certain error between the simulation value and the theoretical calculation value. The maximum error is 4.9%, and within the general error range. From the analysis of metal materials, the material will be plastically deformed due to the applied load exceeding the elastic limit of the material, and then generate residual stress. Due to the existence of residual stress, the distribution of pressure of the contact surface will change, and then produce errors. The use of the finite element method to solve the meshing stiffness has high reliability and accurately reflects the actual meshing condition.

Gear dynamics modeling
The transmission system of spur gear generally includes substructures: gears, transmission shafts, bearings, and box. These substructures constitute the basic frame of gear transmission system [14]. It determines the kinematic and statics characteristics of gear transmission, and is also an important part of determining the dynamic characteristics of the system [15]. The meshing and coupling dynamic model of the gear pair is shown in Fig. 5.
According to the characteristics of both energy storage and energy dissipation of the gear meshing joint, it is equivalent to a spring-damper unit [16][17]. The stiffness and damping of the spring-damper unit represent the gear mesh stiffness and damping, respectively [18][19]. The model is an analysis model of lumped parameter based on discrete mechanical model, from which ordinary differential equations are derived. Regardless of the vibration of the transmission shaft, bearing and box, stiffness of the box, bearing and transmission shaft is expressed by equivalent value kpy, kqy and damping is expressed by equivalent value cpy, cqy.
, q q q q q p q q

Dynamics simulation of gear
The gear parameters of dynamic model are selected as shown in Table, and the active and driven wheels take the same parameters. The rotational speed of the driving wheel is 1140 r/min, and the meshing frequency is 361 Hz, that is to say, the corresponding meshing period is 0.0028 s. The time-varying meshing stiffness of the obtained gear is brought into the differential Eq. (10-13), and the Runge-Kutta method is used to solve the gear meshing dynamic response, as shown in Figs. 6 and 7.
During meshing of gear teeth, periodic response occurs in the process of transmission due to the periodic variation of meshing stiffness. The period is 0.0028 s and the frequency is 361 Hz, corresponding to the vibration condition of the primary engagement.

Dynamics experiment of gear
The gear dynamics experiment was carried out on a CL-100 closed power flow gear test stand. As shown in Fig. 8. The vibration signals in the gear meshing process were collected and dynamically analyzed by sensors and data acquisition systems mounted on the main test chamber. The experiment uses three phase variable frequency motor as the driving device of the system, and the rated speed is 1410 r/min under the voltage of 380 V. In order to achieve a relatively low output speed, a gear unit with a reduction ratio of 15 is provided at the motor output. At the same time, the motor can be adjusted by frequency converter to achieve adjustment within a certain speed range.
This study chooses magnetic powder brake to load and connect with the gearbox output shaft through the coupling. The magnetic powder brake uses the electromagnetic principle and the magnetic powder to transmit torque. In the working range, the loaded torque has a linear relationship with the exciting current, and the torque is kept constant at different rotational speeds. Therefore, the control of the torque can be realized by controlling the magnitude of the exciting current. The experiment uses a program-controlled steady flow source to control the size of the system loading. Fig. 9 is a schematic diagram of the experimental principle. Data Fig. 9 The schematic diagram of the experimental principle The data acquisition system includes data acquisition device and data acquisition program. In this experiment, the NI PXI series data acquisition device and Labview data acquisition program were used to build the experimental platform data acquisition system. Fig. 10 is a block diagram of the overall solution of the data acquisition system of experimental station. PXIe-4330 is connected with TB-4330 junction box to realize 1/4 bridge connection and complete the acquisition of strain signal. The PXIe-6361 collects analog voltage and counter signals, which correspond to torque sensor signals and encoders, respectively. The capture card is installed in the PXIe-1073 chassis and connected to the computer by data cable. Driver software installed on the computer NI Max and Labview data acquisition programs enable control of the data acquisition system. The mechanical structure and data acquisition equipment of the test bench constitute a complete test system.
In the process of the experiment, gear parameters and specific conditions are consistent with simulations. The sampling frequency of the test system is 34.6 kHz, the sampling interval is 10 minutes and the single sampling time is 60 s. Time domain of test signal is shown in Fig. 11, and partially magnified time domain is shown in Fig. 12  It can be seen from Fig. 11 that the period is 0.1 s, which corresponds to a periodic vibration shock with a frequency of 10 Hz, which corresponds to the frequency of the shaft where the drive wheel is located. It can be seen in Fig.  12 that due to the transmission error and the time-varying of the meshing stiffness during the gear transmission, a periodic vibration with a frequency of 600 Hz is caused, which corresponds to a meshing of the gears.

Comparison of experimental data and MATLAB simulation data
The signal measured by the experiment is compared with the signal derived from the simulation of Matlab software. The comparison chart is as Fig. 13-14.  In Fig. 13, it shows the time domain curve of the experiment and simulation. It can be seen from the peak of the two graphs that the experimental and simulation peak are ±350 mm/s 2 . The simulations and experiments coincide well with each other. In Fig. 14, the partial enlargement of the comparison chart is performed. The meshing period from the simulation and the test is 0.0028 s. In summary, there is a good agreement between the simulation results and the experimental results.

Gear crack model of different scales
Gear crack is one of the most important forms of gear failure. Using the stress intensity factor as the basis of crack propagation, the crack propagation angle is determined according to the maximum tangential stress criterion, and a three-dimensional finite element model is established and imported into ABAQUS. The initial crack was placed into the finite element model as shown in Fig. 15. In order to ensure the mesh quality of the tooth surface, the Hex Type Structured and the element type C3D8R were used to refine the mesh around the crack. In the process of setting the gear interaction, the contact pair is defined at the gear meshing point and the contact type is surface-surface contact, as shown in Fig. 16. At the same time, the axial parameters of the crack tip are averaged to determine the direction of propagation of the entire crack tip. The simulation to obtain the entire tooth fracture process is shown in Fig. 17. The thickness of the tooth root is 9mm. I choose 2, 4, 6, and 8mm and other crack lengths and tooth breakage status for stiffness analysis.  It can be seen from Fig. 18 that the occurrence of cracks causes the gear stiffness to decrease, and as the crack length increases, the gear stiffness decreases more. In the single-tooth meshing area, since the faulty tooth acts as a load-bearing tooth, the effect of the crack on the stiffness is most pronounced. In the double-tooth engaging-in and engaging-out area, the effect of the crack on the stiffness is different due to the different bearing conditions of the fault tooth.
In order to visually characterize the trend of the stiffness with crack length in different meshing area. For each crack length, the stiffness of the double-tooth engaging-in area, single-tooth meshing area, and double-toothed engaging-out area is averaged. Then compare with the average stiffness of the corresponding section of the complete gear and define the rate of change of stiffness α. The calculation formula is: where: n=0,1,2 respectively represent the double-tooth engaging-in area, single-tooth meshing area and double tooth engaging-out area.
' op K is the average stiffness of the complete gear for each interval. Taking the crack length as the abscissa and α as the ordinate, the results are shown in Fig. 19.
According to Fig. 19, in the single-tooth meshing area, as the crack length increases, the value of α increases, indicating that the stiffness variation of the single-tooth meshing area is more severe as the crack length increases. At the same time, through the curve of the single-tooth meshing area, it can be seen that the slope of the curve increases with the length of the crack, which means that with the increase of the crack length, the stiffness degradation rate of the single-tooth meshing area is accelerated. In the double-tooth meshing area,  is approximately linear with the crack length, and as the crack length increases, the meshing stiffness decreases, which is consistent with the actual situation. However, it can also be seen from Fig. 19 that the slopes of the curves in the two double-tooth meshing areas are different. The slope of the double-tooth meshing area increases first and then decreases with the increase of the crack length, which indicates that the speed of gear stiffness degradation in the double-tooth engaging-in area increases first and then decreases with the increase of the crack length. The change is most pronounced between crack lengths of 4-8 mm. In the double-tooth engaging-out area, the change of gear meshing stiffness degradation speed is a relatively stable process. Fig. 19 Effect of crack length on stiffness of different meshing area

Conclusions
1. Cantilever beam theory of material mechanics is used to calculate the time-varying meshing stiffness of gear, and the stiffness value calculated by finite element method is compared with the stiffness calculated by theory. The maximum error is 4.9%, and the cause of the error is analyzed. The reason is that the material will be plastically deformed due to the applied load exceeding the elastic limit of the material, and then generate residual stress. Due to the existence of residual stress, the distribution of pressure of the contact surface changes, and then produce errors. The error is within the general error range, which verifies the correctness of the FEM used in this paper.
2. It can be seen from comparison of time domain curve of gear meshing simulation and test that the meshing period is 0.0028 s and the peak value is ±350 mm/s 2 , which indicates that the experiment and simulation coincide well. In the meshing process of gear teeth, the periodic impact response with the meshing frequency is generated due to the time-varying meshing stiffness. It is indicated that the abrupt change of meshing stiffness is one of the main reasons for the vibration shock of gear, and then the abrupt change degree of meshing stiffness can be reduced as an important goal to improve the vibration of gear, so as to achieve the purpose of optimizing the vibration characteristics of gear.
3. The effects of different gear crack lengths on the meshing stiffness are compared. The results show that as the crack length increases, the stiffness decreases. In the singletooth meshing region, the effect of cracks on gear stiffness is most pronounced. In this interval, as the crack length increases, the gear stiffness decreases, and the deterioration rate of the meshing rigidity increases.
4. In the two meshing intervals of the front and rear of the tooth meshing area, the influence of the crack on the gear stiffness is approximately linear. As the crack length increases, the gear meshing stiffness decreases. However, the degradation rate of meshing stiffness is different in different meshing sections. The rate of stiffness degradation of the double-tooth engaging-in area increases first and then decreases, and is most pronounced between crack lengths of 4-8 mm. The stiffness degradation rate of the double tooth engaging-out area is relatively stable.