Study on multi-clearance nonlinear dynamic characteristics of herringbone gear transmission system under optimal 3d modification

The tooth surface 3d modification technology is the most effective method to realize the reduction of vibration and noise of gear system (GS) by comprehensively considering tooth profile modification and axial modification. Moreover, the system nonlinear factors such as backlash and bearing clearance will make the motion of GS have complex variability. Therefore, we will combine 3d modification technology with contact performance and dynamic characteristics to optimize the meshing characteristics of herringbone gear pair through tooth contact analysis technology, loaded tooth contact analysis technology and antlion optimizer with error minimum amplitude as optimization objective, determine three nonlinear dynamic model of herringbone gear for numerical solution. Finally, this paper studies the influence of different system parameters on the system nonlinear dynamic characteristics under 3d modification and non-modification from both local vibration characteristics on the diagrams of time history, phase, spectrum and Poincare map and global vibration characteristics on the diagrams of system bifurcation and maximum Lyapunov exponent. The results show: After optimizing modification, the system vibration is obviously reduced. With the increase of input power, input speed, backlash and static transmission error (STE), the system has jump, while there is no jump with the change of bearing clearance and damping ratio. Compared with other parameters, the changes of input speed and STE make system have complex bifurcation characteristics. The 3d modification can eliminate the jump and make system motion more regular.

Abstract The tooth surface 3d modification technology is the most effective method to realize the reduction of vibration and noise of gear system (GS) by comprehensively considering tooth profile modification and axial modification. Moreover, the system nonlinear factors such as backlash and bearing clearance will make the motion of GS have complex variability. Therefore, we will combine 3d modification technology with contact performance and dynamic characteristics to optimize the meshing characteristics of herringbone gear pair through tooth contact analysis technology, loaded tooth contact analysis technology and antlion optimizer with error minimum amplitude as optimization objective, determine three nonlinear dynamic model of herringbone gear for numerical solution. Finally, this paper studies the influence of different system parameters on the system nonlinear dynamic characteristics under 3d modification and non-modification from both local vibration characteristics on the diagrams of time history, phase, spectrum and Poincare map and global vibration characteristics on the diagrams of system bifurcation and maximum Lyapunov exponent. The results show: After optimizing modification, the system vibration is obviously reduced. With the increase of input power, input speed, backlash and static transmission error (STE), the system has jump, while there is no jump with the change of bearing clearance and damping ratio. Compared with other parameters, the changes of input speed and STE make system have complex bifurcation characteristics. The 3d modification can eliminate the jump and make system motion more regular. The gear rotation angle for axial modification boundary line y 1 and y 3 Maximum modification amount at root and top of tooth y 2 and y 4 Modification length at root and top of tooth. y 5 and y 6 Maximum modification amount of the lower and upper of tooth width. z hf and z ha Modification lengths of the lower and upper of tooth width y 7 The length without modification in axial. d n , d m , d p and d q Modification amounts of the top, root, lower and upper of tooth a ms Modification coefficient of tooth profile u 1 Involute expansion angle of whole tooth profile u 0 , h 0 Starting point of different modified areas k 0 The half angle of base circular groove E 0 The center distance between grinding wheel and gear center E 0 0 Standard center distance a mk Axial modification coefficient h 1 Gear rotation angle S m The coordinate system fixed on gear mounting table S 1 The coordinate system fixed on gear blank S t The coordinate system fixed on grinding wheel ''±'' and '';' ' The upper symbol indicates right tooth surface of the groove, and the lower symbol indicates the left tooth surface of the groove. d The axis non-parallelism error DE Center distance error of driven gears M Coordinate system transformation matrix L The submatrix removing the last row and last column of M. w ik Initial clearance of ellipse center point i before contact w jk Tooth surface clearance at any point j before contact Base radius of passive gear a n Normal pressure angle b The helix angle DT e The normal loaded comprehensive deformation obtained by LTCA Dd The normal GTE obtained by TCA. v s Speed difference of mesh-in point d s Maximum impact deformation k s Meshing stiffness of mesh-in point J 1 and J 2 Moment of inertia e m The average value of errors of gear pair e r Error amplitude of gear pair f s1 (t) and f s2 (t) The meshing impact excitation of left and right helical gear pairs of herringbone gear. k 12L (t) and k 12R (t) The normal time-varying meshing stiffness of left and right helical gear pairs of herringbone gear c 12L and c 12R Normal meshing damping of left and right helical gear pairs of herringbone gear. e 12L (t) and e 12R (t) The STEs of left and right helical gear pairs of herringbone gear 2b 12L and 2b 12R The backlash between left and right helical gear pairs 2c 1L , 2c 1R ,2c 2L and 2c 2R Each radial bearing clearance. f The damping ratio k m The average meshing stiffness m 1 and m 2 The equivalent meshing mass of left and right helical gear pairs M D , C D , K D , q and F D Mass matrix, damping matrix, stiffness matrix, displacement vector and force vector. C D The dimensionless damping matrix K D The dimensionless stiffness matrix F D The dimensionless load array. I A unit matrix P The input power N The input speed

Introduction
Herringbone gear not only has the characteristics of helical gear transmission, but also can self-balance axial force generated in the transmission process. It has the advantages of large contact ratio, large loading capacity and steady movement [1], so is widely used in various fields. Due to machining and installation error, the long-term running wear, there will be inevitable clearances between teeth, and there are also other clearances in the gear system (GS). Nonlinear factors, such as backlash and bearing clearance, can cause the GS to show strong nonlinear characteristics, which lead to chaotic motion, high noise, jump and bifurcation. Therefore, it is of great practical significance to study the nonlinear dynamic characteristics of the GS [2]. The gear transmission system will cause vibration and noise due to the large load fluctuation in the tooth number alternation area. On the other hand, due to the bending and torsional deformation of tooth after being loaded and tooth machining error, the tooth contact is uneven along the tooth width direction, resulting in the partial load phenomenon. Tooth profile modification can compensate base pitch deviation in the actual meshing process of gears, thus slowing down fluctuation of meshing stiffness and reducing meshing impact [3]. The axial modification can effectively improve uneven load along contact line and avoid edge contact, thus increasing bearing capacity of meshing tooth [4]. In contrast, the 3d modification is a comprehensive modification method combining tooth profile modification and axial modification, which can not only eliminate impact, but also avoid edge contact and improve tooth meshing performance [5].
Tooth contact analysis (TCA) is an advanced numerical simulation technology of gear meshing process under light load, which can comprehensively consider the factors, such as modification, installation error and manufacturing error, etc., and obtain tooth surface imprint and transmission error to measure the transmission quality of gear pair [6]. However, when gear pair is under heavy load, tooth surface imprint and transmission error are different from results of TCA, so it is necessary to apply load to simulate meshing process, that is, loaded tooth contact analysis (LTCA). It can comprehensively consider gear geometry error and installation error, etc., and obtain the loaded transmission error (LTE) and tooth surface load distribution coefficient, which is a bridge connecting gear geometry design and mechanical analysis [7].
Kahraman et al. [8] constructed a nonlinear model of three degree of freedom GS with bearing clearance and backlash, studied the nonlinear frequency response characteristics with the time-varying meshing stiffness, and found the transition frequency, the subharmonic response and the chaotic channel. Blankenship et al. [9] established multiple degrees of freedom gear model to analyze the influence of internal and external excitation on the dynamic characteristics of the GS, considering the influence of the backlash and other direction backlash on the model. Kim et al. [10] built the nonlinear model of three degrees of freedom GS and analyzed dynamic response of system with backlash. Zhou et al. [11] analyzed the influence of input speed, backlash and damping coefficient on the system dynamic characteristics based on torsional model of spur gear pair under the excitation of time-varying meshing stiffness. Farshidianfar et al. [12] built the nonlinear dynamic model of the GS, and then studied influence of backlash, time-varying stiffness and external excitation on bifurcation characteristics and chaos of system by solving the dynamic model analytically. Li et al. [13] proposed a more realistic dynamic meshing force model based on the coordinated deformation of gear teeth, which was suitable for the modified and nonmodified GS and verified by the finite element method. The research showed that vibration amplitude of GS in bifurcation region, quasi-periodic region and chaotic region was greater than that of the geometric transmission model. Shi et al. [14] put forward a dynamic model of hypoid gear, calculated the time-varying meshing stiffness based on dynamic meshing force, and studied the influence of meshing stiffness on gear system dynamic response. The result showed that the meshing stiffness had obvious influence on dynamic characteristics of hypoid gear. Considering the effects of nonlocal mesoscopic deformations nonlinear characteristic of confined granular systems, Gonzalez et al. [15] presented a nonlocal formulation of contact mechanics that accounted for the interplay of deformations due to multiple contact forces acting on a single particle. The calculation results of the formulation were the same as finite-element simulations and experimental observations. Wang et al. [16] adopted multiscale simulation of molecular dynamics to investigate the nonlinear forced vibration of PMMA/ CNT composite plate under coupled temperature and aerodynamic load conditions. The partial differential equation of composite plate was derived by nonlinear strain-displacement relation and Hamilton's principle, which was discretized into a set of coupled ordinary differential equations by applying Galerkin method and was solved by the fourth-order Runge-Kutta method. Yadav et al. [17] studied the nonlinear vibration characteristics of the spring-dashpot mechanisms. The system qualitative behavior with deteriorated conditions of low viscous damping, higher backlash and increased dry friction was established through numerical simulations. When existing dry friction and backlash at the same time, the comprehensive influence of system parameters on the transition to chaotic motion was analyzed. Liu et al. [18] studied the nonlinear forced vibrations of functionally graded material sandwich cylindrical shells with porosities on an elastic substrate. The partial differential equations of the system were obtained by an energy approach based on the Donnell's nonlinear shallow shell theory and Hamilton's principle. Then, the multi-degree of freedom nonlinear ordinary differential equations were built by Galerkin method. Finally, the pseudo-arc length continuation method was utilized to perform the bifurcation characteristics of system and analyze the effects of system parameters on its vibration characteristics. Further, Liu et al. [19] also proposed a new solution approach to solve the nonlinear forced vibrations of functionally graded piezoelectric shells in multi-physics fields. Based on Hamilton's principle and the Donnell nonlinear shallow shell theory, they obtained the motion equations. The multi-mode Galerkin method and the pseudo-arc length continuation method were used to solve the nonlinear vibration characteristics of the multi-degree of freedom system. The results showed that the external applied voltage, temperature change, external excitation and other parameters had important roles on nonlinear vibration response and bifurcation analysis of FG piezoelectric shells with micro-voids. In addition, Liu et al. [20,21] studied the nonlinear forced vibration of composite cylindrical shells in multiple physical fields and the nonlinear forced vibration of rotating shell under harmonic excitation in thermal environment. The research showed that the nonlinearity of shell structure was material nonlinearity, which was caused by the nonlinearity of the relationship between stress and strain of materials. Nonlinear analysis of shell structures need to require the calculation of kinetic energy, strain energy and work done by external forces, and then the partial differential dynamic equations of shell structures were obtained based on Hamilton principle. Then, the equations were discretized by Galerkin method to form the corresponding system coupling degrees of freedom nonlinear ordinary differential equations. Finally, the pseudo-arc length continuation algorithm were used for numerical simulation. Wang [22] established the relationship between tooth surface modification, tooth surface load distribution and system dynamic model, and proposed a multi-objective modification optimization method with noise reduction and load distribution, and verified it by an example. The results showed that the root mean square (RMS) of acceleration after being optimized was reduced by nearly 4 times. Based on TCA and LTCA technology, Wang et al. [23] put forward a tooth surface optimization modification method to improve comprehensive characteristics of gears with the minimum tooth surface flash temperature and minimum RMS of vibration acceleration as the objectives, and tooth surface contact was obviously improved after optimization. Yang et al. [24] proposed a dynamic excitation calculation method of 3d modified herringbone gear pair considering actual contact state, and studied the impact of modification parameters on timevarying meshing stiffness excitation and friction excitation. After modification, the vibration and noise of herringbone gear system were obviously improved. Nishino [25] proposed a 3d modification method to study the relationship between contact conditions and meshing stiffness of helical gear, and the meshing performances of helical gear were improved after modification. Jian et al. [26] studied vibration characteristics of the GS with various nonlinear factors, and obtained the system dynamic vibration forms such as the periodic response, subharmonic and chaotic response by adjusting gear system parameters. Li et al. [27] used Poincare mapping, Lyapunov exponent spectrum, spectrum analysis and other means to comprehensively judge the chaotic vibration form of gear transmission system. Lin et al. [28] analyzed the system global bifurcation behavior when the parameters such as meshing frequency, backlash, comprehensive transmission error and damping ratio changed in torsional-parallel gear transmission system, and investigated the evolution law of maximum Lyapunov exponent (MLE) with system parameters. Using diagrams of Poincare map, the phase and spectrum to analyze locally the influence of parameters on dynamic characteristics.
Based on the above analysis, we find that the study on the nonlinear dynamics of gear system only considers backlash, and the excitation only includes stiffness excitation or impact excitation. There is little research on the nonlinear dynamic characteristics of multi-clearance GS under the comprehensive consideration of three internal excitations of system, especially for herringbone gear. The research of gear modification technology mostly involves the meshing performance of gear pair, and rarely studies the influence of the optimal modification parameters on nonlinear dynamic characteristics of GS. In addition, the global influence of system parameters on the dynamic characteristics of GS is also rare, and the research on the influence of system parameters on chaotic motion state is also insufficient. There is no literature study on whether optimal modification parameters meet the requirements of nonlinear vibration characteristics of the GS.
The 3d modification of herringbone gear tooth surface can eliminate the edge contact state of tooth surface under actual working conditions, and different modification ways can change meshing marks and geometric transmission errors, thus affecting the timevarying meshing stiffness and transmission errors of gear system, and further affecting the system internal excitation, resulting in complex and varied vibration characteristics of gear system. Actually, there are many kinds of clearances in GS, but the existence of these clearances is reflected in gear dynamic equation as strong nonlinear terms. Therefore, it is necessary to study the system nonlinear vibration characteristics with multiple clearances under optimal modification. Based on the above shortcomings, the forming grinding principle proposed in this paper can realize the 3d modification of herringbone gear mainly through the two directions of the grinding wheel's spiral rise along the tooth direction of the teeth and its movement along the radial parabola of the teeth. In addition, different grinding wheel profiles and motion trajectories can achieve different gear modification methods. The terminal profile equations of tooth profile sectional modification and the tooth profile of grinding wheel are determined. The radial motion equation of grinding wheel along the workpiece during axial sectional modification is derived. Based on these, the 3d modified tooth surface equation is determined. Combining TCA and LTCA technology, the optimal modification parameters are obtained by the Antlion optimizer (ALO) with the minimum amplitude of the loaded transmission error (LTE) as the optimization objective. On this basis, the internal system excitation is determined, and the bendingtorsion-axial coupled multi-clearance nonlinear dynamic model with time-varying meshing stiffness excitation, meshing impact excitation and error excitation is established. Finally, from the global perspective of system, based on the diagrams of system bifurcation and the maximum Lyapunov exponent (MLE), and from the local perspective of system, based on diagrams of time history, phase, spectrum and Poincare map, the influences of system parameters on nonlinear dynamic characteristics of the GS under the unmodified and modified conditions are studied, respectively, which provides a basis for the selection and optimization of system parameters.
The 3d modification of herringbone gear is realized through the grinding forming principle, which improves the theoretical level of gear modification and the application value of modification technology. On this basis, combined with TCA and LTCA technology, three kinds of internal excitation of GS are determined, and the influence of system parameters on the system nonlinear vibration characteristics under multi-clearance is studied, which can improve the dynamic and matching characteristics of transmission system and box body, so as to control the noise of GS, realize the purpose of vibration and noise reduction of transmission system, prolong the service life of gear and reduce the influence of system error, and improve the engineering application value of herringbone gear transmission system. On the basis of gear modification, the system nonlinearity is studied from two aspects of local vibration characteristics and global vibration characteristics of GS, which means combining gear modification technology with system nonlinearity, thus enriching the nonlinear theory of gear system.

Tooth surface modification method
Form grinding is an efficient grinding method for tooth profile modification and axial modification of gear. Based on the requirements of actual working condition, bearing capacity and dynamic behavior, the tooth surface of herringbone gear can be modified by changing the modification parameters, and the 3d modification can be realized from a micro perspective [29]. It includes tooth profile modification and axial modification, as shown in Fig. 1. In Fig. 1, P 1 represents 3d modified tooth surface, and P 2 represents theoretical involute tooth surface. u o , u p , u q and u a are the involute expansion angles at tooth profile modification boundary line. h o , h n , h m , and h max are gear rotation angle for axial modification boundary line, respectively. y 1 and y 3 are maximum modification amount at root and top of tooth, y 2 and y 4 are modification length at root and top of tooth. y 5 and y 6 are maximum modification amount of the lower and upper of tooth width. z hf and z ha are modification lengths of the lower and upper of tooth width. y 7 is the length without modification in axial. d n , d m , d p and d q are modification amounts of the top, root, lower and upper of tooth, respectively. In this paper, by controlling involute expansion angle and gear rotation angle corresponding to boundary line of tooth profile modification and axial modification, the division areas of 3d modified tooth surface can be realized, that is, combination of the two parabola and a straight line modification design or parabola modification design of whole tooth surface can be completed.

Equation of 3d modified tooth surface
Because herringbone gear can be regarded as ''double helical gear'', the tooth profile modification forming combining parabola line at both tooth top and tooth root and straight line can be realized by superimposing the modification amount DH on generating line of involute, as shown in Fig. 2. s'f' is standard tooth profile, and sf is modified tooth profile. cg is standard involute, u c and u g are expansion angles at c and g, respectively. The relationship between modification amount and expansion angle is showed as: where, a ms is modification coefficient of tooth profile, and r b is base circle radius. u 1 is involute expansion angle of whole tooth profile, which takes different values in different modified areas. u 0 is starting point of different modified areas, u 0 = u c at tooth root and where, k 0 is half angle of base circular groove, '' ± '' and '' ; '', upper and lower symbol represent right and left involute.
Changing motion trajectory of grinding wheel along radial direction of workpiece to realize sectional axial modification, as shown in Fig. 3. When gear turns h 1 , grinding wheel moves a distance p 1 h 1 along axial z 0 , here, p 1 is gear pitch, and its feed rate is l x along radial x 0 . The trajectory form of grinding wheel along the axial is the parabola, the straight and the parabola, among which on and mi are parabola line, and nm is straight line. h o , h n , h m and h max are gear rotation angle at point o, n, m and i along axial, respectively. h is the axial distance between grinding wheel axis and starting position of modification, that is, h = p 1 h 1 . E 0 is center distance between grinding wheel and gear center, that is, Then radial feed rate l x of grinding wheel along workpiece can be expressed as: where, a mk is axial modification coefficient, h 1 is gear rotation angle, which takes different values in different modified areas. h 0 is starting angle of different modified areas, h 0 = h n at tooth lower and h 0 = h m at tooth top. In addition, modification amount of grinding wheel along the radial needs to be transformed into the normal amount of tooth surface according to projection transformation relationship between the two directions of gear tooth surface.
In this paper, the active gear is modified while the passive gear is not modified (that is, the active is used to trim the root instead of the passive to trim edge). When tooth surface is 3d modified, the grinding wheel moves in two directions: spiraling up along the axial of tooth and making parabolic motion along the radial of tooth. Therefore, according to the terminal section profile of grinding wheel, the tooth surface equation of 3d modification can be obtained by the coordinate transformation. The coordinate relation of single helical tooth of herringbone gear is shown in Fig. 4. The coordinate system S m is fixed on gear mounting table, S 1 is fixed on gear blank, z 1 passes through the axis of mounting table and blank, and S t is fixed on grinding wheel. The grinding wheel makes a spiral motion around z m and a radial parabola motion along x m to form a 3d modified tooth surface. When coordinate system S 1 rotates h 1 around z m , the origin of S t moves in two directions: moving distance h = p 1 h 1 along z m and negative feed distance l x along x m .
As shown in Fig. 4, 3d modification vector of herringbone gear with single helical tooth is showed as: Equation (4) is expanded to obtain: After simplifying Eq. (6), we can get: Here, The unit normal vector of 3d modified tooth surface is: It is worth noting that there exists c m between coordinate system S t and S m , that is, the installation angle of grinding wheel, which can be obtained according to gear helix angle b. Therefore, 3d modified tooth surface equation of herringbone gear is obtained by changing the helix angle b direction and adjusting center distance E 0 in the same way as above.
3 Contact analysis of 3d modified tooth surface

Tooth contact analysis(TCA) of herringbone gear
In order to numerically simulate two intermeshing herringbone gears, the herringbone gear pair is regarded as two pairs (I and II) of single helical gears meshing. According to the solution method of helical tooth contact analysis [30], tooth contact analysis of tooth pairs I and II is carried out by transforming meshing coordinate system to reference coordinate system, as shown in Fig. 5. As shown in Fig. 5, S f is the fixed reference coordinate system located at the midpoint of tooth groove, S 0 i and S i (i = 1,2,3,4) are fixed coordinate system and rotation coordinate system located at the width midpoint of helical gear of tooth pairs I and II, respectively, and S m is the auxiliary coordinate system. When the herringbone gear is p and center distance error DE of driven gears are inevitable, which are expressed by auxiliary coordinate systems S 00 3 and S 00 4 , respectively. In coordinate system S f , the position vector and normal vector of tooth surface of herringbone gear pair are showed as: where, M is coordinate system transformation matrix, and L is the submatrix removing the last row and last column of M.
During meshing process of herringbone gear pair, the two tooth surfaces are in continuous tangent contact, and the contact process is expressed as: where, ''1,2'' means active gear, ''3,4'' means passive gear, and ''f'' means that all vectors are in fixed coordinate system S f . When herringbone gear enters meshing and exits meshing, the edge contact occurs on both sides of the tooth surface and the tooth top, and the contact process is expressed as: where, r fi represents position vector of edge contact point between tooth side of active gear and tooth top of active gear,or fi =ou ci and or fi =oh ci represent the tangent directions at the edge of tooth side of active gear and tooth top of active gear. The geometric transmission error (GTE) of herringbone gear pair is the variation between the actual rotation angle and theoretical rotation angle of the passive gear with rotation angle change of the active gear [31], that is expressed as: where, Z 1 and Z 2 are tooth number of active and passive gear, u j (0) and u i (0) are initial positions of active and passive gear.

Loaded tooth contact analysis (LTCA) of herringbone gear
The contact process of herringbone gear pair is simplified into the geometric model of cross-section contact of left and right single helical gear pairs along long axis of contact ellipse, as shown in Fig. 6. i and j represent instantaneous contact ellipse center and any discrete point. w ik represents initial clearance of ellipse center point i before contact, w jk represents tooth surface clearance at any point j before contact. p ik and p jk represent unit normal load of contact point corresponding to the long axis of ellipse, d jk represents tooth surface clearance after deformation at any contact point j under load, k = L, R.
The mathematical model of loaded tooth contact problem composed of the displacement coordination equation, the force balance equation and non-embedded condition of herringbone gear is showed as follows: where, [w] k is initial clearance of tooth surface before being loaded, which is composed of inter tooth clearance and normal clearance.
[d] k is tooth surface clearance after being loaded.
[Z] is tooth rigid body displacement.
[L] k [P] k represents the elastic deformation obtained by multiplying normal flexibility matrix and normal load.
The loaded transmission error (LTE) in a meshing period is the angular displacement that converts [Z] to the meshing line, which is shown as: where, r b2 and b are base circle radius and helix angle of the passive gear.

Optimization model of modification parameters
The optimization of the optimal modification parameters of tooth surface is an iterative improvement process of nonlinear contact problems, which is a process of changing contact state of tooth surface through micro-geometric iteration of different modification parameters until a specific goal is met [32]. The loaded transmission error reflects deviation between actual meshing and ideal meshing of gears.
The larger its value, the smaller stiffness excitation, which is direct excitation of vibration. Thus, this paper takes the minimum amplitude of loaded transmission error as the optimization objective. The modified way adopts two parabolas and one straight line in tooth profile and axial. Because the axial modification adopts symmetrical modification method, six modification parameters are taken as optimization variables.
Variables: x ¼ ½y 1 ; y 2 ; y 3 ; y 4 ; y 5 ; y 7 Objective : maxðDT e Þ À minðDT e Þ j j ð 15 À 1Þ Constraint condition : s:t where, x is 3d modification parameter. F(x) is the loaded transmission error amplitude calculated by Eq. (14). C min and C max are minimum and maximum modification amount of tooth profile at tooth root and tooth top. D min and D max are minimum and maximum modification length of tooth profile at tooth root and tooth top. S min and S max are minimum and maximum modification amount at the upper end and lower end of axial direction. G min and G max are the minimum and maximum values of non-modified length of axial direction (axial modification is symmetrical, so only the non-modified length is taken). In view of poor global search ability of traditional algorithms in solving multi-peak practical problems, seyedali [33] first proposed a novel nature-inspired algorithm called Antlion Optimizer (ALO) in 2015. It has two types of individuals: Ant and Antlion, and realizes five steps: the random walk of ants, building traps, entrapment of ants in traps, catching preys, and re-building traps. The proposed algorithm is Tooth groove width (mm) 50 benchmarked in three phases: 19 mathematical functions, three classical engineering problems and the shapes of two ship. The research results showed that the ALO can provide very competitive results in terms of the improved exploration, local optima avoidance, and convergence. Moreover, it has merits in solving constrained problem with diverse search spaces. The Antlion optimizer optimizes the problem by using the hunting relationship between antlion and ants, and completes the search of the search space with the help of the ants walking randomly around the antlion. The Antlion catches ants to realize population evolution. In the process of evolution, the antlion learns from elite antlion to ensure the diversity of population. The antlion represents the solution of function to be optimized, and the ants with high hunting fitness update and save the approximate optimal solution [34].
Taking basic parameters of herringbone gear pair in Table 1 as the research object, based on the obtained 3d modified tooth surface, combined with TCA and LTCA technology, this paper adopts the ALO to optimize the modification parameters with taking the minimum amplitude of the LTE as the optimization objective. The rated load is 500Nm, input speed is 4000 rpm, the maximum number of iterations is 100, the number of variables is 6, and the objective function is fobj = @F(x). The iterative optimization process of modification parameters is shown in Fig. 7. The results of modification parameters after optimization are shown in Table 2. The tooth surface of active gear is evenly divided into nodes along tooth height and axial direction, and tooth profile modification curve, the axial modification curve and the fitted 3d modified tooth surface after optimization are shown in Fig. 8.
Taking basic parameters of herringbone gear pair in Table 1 as research object, the tooth surface contact simulation of the standard, non-parallelism error d = 0.5 00 and center distance error DE = 0.1 mm, and 3d modification is calculated to get meshing impression and GTE, as shown in Fig. 9, 10 and 11.
As shown in Fig. 9, the standard tooth surface is mainly in contact with tooth surface and the impression is symmetrical equal, and the GTE is relatively  There is obvious tooth edge contact, and the difference of the GTE between the left and right tooth is obvious. The contact occurs partial load and will cause vibration. As shown in Fig. 11, the 3d modified tooth contact with d = 0.5 00 and DE = 0.1 mm does not have tooth side and tooth top edge contact. The unmodified contact path is distributed along pitch circle, and the GTE of left and right tooth surfaces is symmetrical equal and parabolic. In contrast, the 3d modification can not only improve contact performance of tooth surface, but also eliminate the influence of installation error.
As shown in Fig. 12, when d = 0.5 00 and DE = 0.1 mm, the variation curves of comprehensive LTE on the meshing tooth surface of the active gear before and after 3d modification optimization. The amplitude of the LTE of the standard tooth surface is 4.194 00 . After modification, the amplitude of the LTE is 1.157 00 , and its fluctuation amplitude is reduced by 72.41% and becomes relatively gentle.

Multi-clearance dynamic model of herringbone gear pair
Analyzing vibration characteristics of GS under internal and external disturbances can not only provide reliable parameter selection for practical engineering, but also provide a basis for attenuating vibration and noise. In gear dynamics system, the loaded teeth will produce the elastic deformation, and different elastic deformation may cause system nonlinear characteristics. However, when calculating the time-varying meshing stiffness of gear pairs, it can be obtained by loaded tooth contact analysis (LTCA) technology, and the tooth elasticity is considered in the time-varying meshing stiffness. Considering that the plastic deformation of the tooth will cause interference and fracture between teeth, this paper ignores the nonlinearity it brings. In addition, the friction effect between teeth can be ignored under the high-speed operation of gear. Therefore, based on TCA and LTCA technology of the modified tooth surface, the bending-torsion-axial   The change curves of LTE before and after 3d modification coupled multi-clearance nonlinear dynamic model with stiffness excitation, impact excitation and error excitation is established in this paper.

Stiffness excitation
The resistance to the loaded tooth deformation is called the stiffness. The section calculates timevarying meshing stiffness of gear pair based on the calculation results of modified tooth surface TCA and LTCA technology. Considering the size effect caused by modification, it can be considered in the timevarying meshing stiffness. When calculating the LTE, the load shared by multiple pairs of meshing teeth and the initial clearance of tooth surface should be considered at the same time. Therefore, the LTE includes the GTE and comprehensive elastic deformation of the tooth surface [35]. According to the definition of stiffness, it is expressed as: where, T is output torque, r b2 is base radius of passive gear, a n is normal pressure angle, b is helix angle, DT e is the normal loaded comprehensive deformation obtained by LTCA. Dd is the normal GTE obtained by TCA.

Impact excitation
Because modification destroys conjugate characteristics of the standard involute, the mesh-in point of meshing tooth pair will deviate from the theoretical meshing line. The specific position of off-line mesh-in impact, whether on the active gear tooth surface or passive gear tooth surface, cannot determine its position through the geometric relationship of gear pair. Therefore, based on the tooth contact analysis, this paper uses contact simulation model [36] of meshing tooth with base pitch error caused by being loaded to calculate accurate position and meshing stiffness at mesh-in point based on LTCA technology. At the same time, calculates the relative speed of mesh-in point with tooth surface modification and base pitch error, and finally calculates the meshing impact force.
According to the theory of the impact mechanics, the impact kinetic energy is equal to the elastic potential energy, and its expression [37] is showed as: where, v s is speed difference of mesh-in point, n is obtained by data fitting, d s is maximum impact deformation, k s is meshing stiffness of mesh-in point, J 1 and J 2 are moment of inertia, r b1 and r b2 are radius of base circle of the active and passive gear, respectively. Finally, according to the formula of maximum impact force: F s = k s d s , the meshing impact force is obtained as:

Error excitation
Due to the lack of manufacturing level, tooth profile deviation and tooth pitch deviation will occur, which cause gear tooth to deviate from theoretical motion trajectory. The design errors with the manufacturing transmission error and the randomness are called static transmission error (STE) in this paper. In order to facilitate calculation without losing accuracy, the simple harmonic function is used to represent the STE [38], as follows: eðtÞ ¼ e m þe r sinðw n t þ uÞ ð 19Þ where, e m and e r are the average value and the amplitude of comprehensive errors of the gear pair.
Generally, e m = 0,e r ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffi , f a and f b are base pitch error and tooth profile error, which are determined by gear accuracy grade. w n is the mesh frequency. u is initial phase angle of error, usually, u = 0.

Establishing the bending-torsion-axial
coupled multi-clearance nonlinear dynamic model of herringbone gear pair The nonlinear factors, such as backlashes, bearing clearances, time-varying meshing stiffnesses and other nonlinear factors, can cause the GS to show strong nonlinear characteristics, resulting in chaotic vibration, large noise, jump, accuracy reduction and other phenomena. Therefore, comprehensively considering the effects of stiffness excitation, impact excitation, error excitation, backlashes and bearing clearances, this paper adopts the lumped parameter method to establish the twelve degree of freedom bending-torsion-axial coupled nonlinear vibration model of herringbone gear pair, as shown in Fig. 13. As shown in Fig. 13 are the mass and moment of inertia of the left and right helical gears (each including half of relief groove) of the active and passive herringbone gears. T 1L and T 1R are input torque, T 2L and T 2R are load torque, which are regarded as left and right equality. f s1 (t) and f s2 (t) are the meshing impact excitation of left and right helical gear pairs of herringbone gear. k 12L (t) and k 12R (t) are the normal time-varying meshing stiffness of the left and right helical gear pairs of herringbone gear. c 12L and c 12R are the normal meshing damping of the left and right helical gear pairs of herringbone gear. e 12L (t) and e 12R (t) are the STEs of left and right helical gear pairs of herringbone gear. 2b 12L and 2b 12R are the backlashes between the left and right helical gear pairs, respectively. 2c 1L , 2c 1R , 2c 2L and 2c 2R are each radial bearing clearance. k 1LRz , k 2LRz , c 1LRz and c 2LRz are the axial equivalent stiffness and damping of the middle relief groove of the left and right helical gears, respectively, which are equivalent to the solid short shaft to calculate. k iy , k iz , c iy and c iz (i = 1L, 1R, 2L, 2R) are the radial and axial equivalent bearing stiffness and damping, respectively. b L and b R are helix angles of left and right helical gears. r b1 and r b2 are base radius of active and passive gear.
The generalized displacement array q of gear system is expressed as: Here, y i , z i and h i (i = 1L,1R,2L,2R) are the translational vibration displacement and angular displacement of the concentrated mass points O 1L , O 1R , O 2L and O 2R of the active and passive herringbone gears in the y and z directions, respectively.
Combined with Fig. 13, according to Newton's laws of mechanics, the system dynamics equations are showed as: The tangential and axial dynamic meshing forces on the left and right sides are expressed as: The relative vibration displacements in the meshing line direction of the left and right end faces are written as: The normal meshing damping c 12L and c 12R of left and right helical gear pairs of herringbone gear are expressed by [39] as followed: Here, f is the damping ratio, k m is the average meshing stiffness.
The clearance functions on the meshing lines of left and right helical gear pairs are showed as followed [40]: where, k i are the relative displacement and 2b i are the backlash.
The clearance functions f i (y i , c i ) of the radial bearing clearances are written as followed [41]: where y i are radial displacement, 2c i are the bearing clearance.
Since the GS has rotational degrees of freedom, the dynamics equations have the rigid body displacement in the torsion direction, so the stiffness matrix of the equations is a positive semidefinite matrix. If the equations are solved numerically directly, the iteration may not converge and the correct response cannot be obtained. In order to eliminate rigid body displacement of system on rotation angle, the relative displacement method is introduced to reduce the rigid body displacement of dynamic equations.
Taking q 1 ¼ r b1 h 1L À r b2 h 2L ,q 2 ¼ r b1 h 1R À r b2 h 2R as the relative displacement, the torsional rigid body vibration displacement in Eq. (21) can be transformed into: Where m 1 ¼ I1LI2L are equivalent meshing mass of left and right helical gear pairs [42].
Then, combining with translational vibration displacement, the nonlinear dynamic equations of system can be written as: From the derivation process of the equation, the relative displacement method is equivalent to a coordinate change of the original dynamic equation. The transformation process does not change inherent characteristics of the GS itself. Therefore, the results obtained by this method are completely loyal to the original system and truly reflect the dynamic characteristics of original dynamic model. After simplification, the dynamic equations of herringbone gear system with 10 degrees of freedom can finally be obtained as followed: where M D , C D , K D , q and F D are mass matrix, damping matrix, stiffness matrix, displacement vector and force vector.
Then Eq. (29) is subjected to dimensionless treatment. Define dimensionless time s = tx n and displacement nominal scale b c . Here, x n can take the first natural frequency of single helical gear pair, and b c can take the unit displacement of the same order of magnitude as the static transmission error. The dimensionless displacement q¼b c d, dimensionless velocity _ q¼b c x n _ d and dimensionless acceleration € q¼b c x 2 n € d can be obtained. Then they are brought into the dynamic differential equation group Eq. (29). After simplification, the dimensionless bending-torsion-axial coupled nonlinear dynamic equation of herringbone gear pair with 10 degrees of freedom can be obtained: where D F D is the dimensionless load array.
Assuming that I is a unit matrix and uðsÞ ¼ ½dðsÞ d ! ðsÞ T is a new state vector, Eq. (30) can be rewritten as: According to the above derivation, Eq. (31) can be directly solved by using the fourth-fifth-order adaptive Runge-Kutta iterative solver ode45 provided by MATLAB programming, and the steady-state response of nonlinear dynamics of system can be obtained by ignoring its initial response. The numerical iteration method is mainly the direct numerical integration method in the time domain, which can accurately solve the steady-state and transient responses of system equations. In addition, this method has high adaptability to dynamics considering various system ð28 À 2Þ factors, and can accurately solve the steady-state and transient responses of equations in gear system.

Influence of different parameters on gear nonlinear dynamic characteristics
Because of various nonlinear factors, the vibration state of herringbone gear system is rather complicated. The system dynamic characteristics will change as the different system parameters change [43]. Therefore, when a single system parameter changes and other parameters are set unchanged, by combining the global analysis method and the local analysis method, the nonlinear dynamic characteristics of the GS before and after 3d modification are studied, and the corresponding diagrams of bifurcation, maximum Lyapunov exponent (MLE), time history, phase, spectrum and Poincare map are obtained to reflect the system dynamic characteristics.

Influence of input power on system dynamic characteristics
It is of great guiding significance to study the influence of external excitation on the dynamic behavior of the GS. By controlling the external motor current or motor voltage and input torque, the system can maintain a stable operation state and suppress system vibration.
The input power is provided by the motor, and the external excitation is different when the power value is different. Considering the actual engineering needs and ignoring the output torque fluctuation, when the power value is different and the input speed is constant, the input torque value is different according to T = 9550ÁP/n. Therefore, taking basic parameters in Table. 1 as research object, d = 0.5'',DE = 0.1 mm (the values are same when analyzing the following parameters), the input power P varies from 10 to 2000 Kw, the others are shown in Table 3. The change curves of the root mean square (RMS) values of the relative comprehensive vibration accelerations of left and right meshing lines of herringbone gear pair before and after 3d modification with the change of P are obtained, as shown in Fig. 14. According to Fig. 14, with the increase of P, the RMS values of acceleration before and after 3d modification change from big to small and then tend to be stable. The jump occurs at about 550 Kw before modification and about 420 Kw after modification, and the RMS values of vibration acceleration after modification are greatly reduced. The diagrams of system bifurcation and the MLE of dimensionless comprehensive relative vibration displacements of left and right meshing lines of herringbone gear pair with the change of P before and after 3d modification are obtained, as shown in Fig. 15. As can be seen from Fig. 15(a), with the increase of P from 10 Kw, the system shows chaotic motion, and the MLE is   When P reaches about 640 Kw, the system enters a stable single periodic motion, the displacement increases continuously, and the MLE is less than 0. As can be seen from Fig. 15(b), after modification, with the increase of P from 10 kW, the system shows chaotic motion, and the MLE is greater than 0. When P reaches about 360 Kw, the system bifurcates into twice periodic motion, the MLE changes from positive to negative, and there is jump at about 420 Kw. When P increases from 420 to 920 Kw, the system enters chaotic motion, and there is a periodic window. After about 920 Kw, the system has been in single periodic motion, and the MLE is less than 0. In contrast, after modification, with the increase of P, the system motion state undergoes the change of ''chaos-double period-chaos-single period'', showing rich nonlinear characteristics. In addition, the modification can eliminate the jump, but increase chaotic interval. Therefore, when P must be more than 920 Kw, chaotic motion can be effectively avoided.
In order to further describe system vibration characteristics shown in Fig. 15, the local analysis method is used to analyze corresponding response under the specific input power P after modification through the diagrams of time history, phase, spectrum and Poincare map. As shown in Fig. 16, when P = 50 Kw, the time history diagram has no periodicity, phase trajectory keeps winding and fills a closed area, spectrum diagram is a continuous spectrum band, and the Poincare map diagram is composed of dense points, so the system is chaotic motion. As shown in Fig. 17, when P = 1500 Kw, the time history diagram has obvious periodicity, phase diagram is a closed curve, spectrum diagram is obviously discrete, and the Poincare map diagram has a attractor, so the system is a single periodic motion.

Influence of input speed on system dynamic characteristics
Considering that the actual input speed fluctuation will change the meshing frequency, meshing error and centrifugal force of gear, which may make the GS have very complex dynamic characteristics, the influence of input speed on system dynamic characteristics should be studied. Studying the influence of input speed on the GS under the complex excitation can provide the cognition for system dynamic characteristics in engineering practice and verify the necessity of gear modification. Therefore, when writing the dynamic program, the input torque and output torque are controlled to remain invariant at the rated speed, so that only the input speed can take different values to analyze its influence on system dynamic characteristics. Take basic parameters in Table 1 as research object, the input speed n varies from 100 to 8000 rpm, the others are shown in Table 4.
The change curves of the RMS values of the relative comprehensive vibration accelerations of left and right meshing lines of herringbone gear pair before and after 3d modification with the change of n are obtained, as shown in Fig. 18. As can be seen from Fig. 18, with the increase of n, the system appears multiple resonance peaks before and after 3d modification, which corresponding to system resonance frequency and 1/2 resonance base frequency respectively, and appears several jumps at low speed. The modification can reduce the low-speed vibration of system, and increase the high-speed vibration of system, which has a great impact on the high-speed vibration.
The diagrams of global and local bifurcation and the MLE of dimensionless comprehensive relative vibration displacements of left and right meshing lines of herringbone gear pair with the change of n before and after 3d modification are obtained, as shown in Fig. 19-20. As can be seen from Fig. 19, 20, with the increase of n, the system shows very complex and changeable bifurcation characteristics, as follows: In Fig. 19(a), when n increases from 100 to 1000 rpm, the system shows chaotic motion and the MLE is greater than 0. The MLE is less than 0 after a short multiple cycle at about 1000-1440 rpm. Then the system enters a short chaotic motion until 1750 rpm. From 1750 to 1940 rpm, it enters a short period motion with a peak, and then enters a chaotic motion again. Until about 2200 rpm, the MLE is greater than 0. Then  the system diverges into a short twice periodic motion, then enters chaotic motion, then enters single periodic motion until 2800 rpm, and then enters chaotic motion until 3015 rpm. The system is in single periodic motion since about 4000 rpm. At about 4650 rpm, it appears large resonance peak, and peak value is about 96.43 m/s 2 . In Fig. 19(b), after modification, when n increases from 100 to 1000 rpm, the system shows chaotic motion and the MLE is greater than 0. However, the MLE is less than 0 when it undergoes multiple cycles from 1000 to 1478 rpm. Then the system enters chaos until 1750 rpm. From 1750 to 1955 rpm, it enters a short period with a peak, then enters a chaotic movement again until about 2200 rpm, and the MLE is greater than 0. Then the system diverges into a short twice periodic motion, then enters chaotic motion. The system enters single periodic motion until 2789 rpm, then enters periodic motion until 3042 rpm. It has been in single periodic motion since about 4000 rpm. At about 5060 rpm, the system appears a large resonance peak, and the peak value is about 108.3 m/s 2 . In contrast, the periodic motion exists in chaotic motion before modification, and the system motion is very complex. After modification, the chaotic region of system will increase, but periodic motion in chaotic motion will be eliminated to make the system motion regular. In addition, over 4000 rpm, the modification can eliminate system bifurcation characteristics.
To sum up, the change of input speed has a great impact on the system, and the system nonlinear characteristics are very sensitive for input speed, which shows complex and changeable nonlinear vibration characteristics. During the change of n from 100 to 8000 rpm, the system undergoes the chaotic motion, complex periodic motion, and appears many jumps and resonance peaks. The channel into chaotic motion can be through the jump or multiple period phenomenon, etc.
Similarly, for Fig. 19, 20, the local analysis method is used to further analyze system motion characteristics under specific parameter after modification. As shown in Fig. 21, when n = 1500 rpm, time history diagram has no periodicity, the phase diagram shows that the several phase lines are intertwined, the spectrum diagram is a continuous spectrum band, and Poincare map diagram shows dense points in pieces, at which time the system moves in chaos. As shown in Fig. 22, when n = 2200 rpm, the time history diagram has no obvious periodicity, the phase diagram randomly fills a closed area, the spectrum diagram is a continuous spectrum band, and the Poincare map diagram is composed of patches of dense points, at which time the system moves in chaos. As shown in Fig. 23, when n = 7000 rpm, the time history diagram is periodic, phase diagram is a closed curve, spectrum diagram has the obvious discreteness, and Poincare map diagram has an attractor, so the system moves in a stable single period.
6.3 Influence of clearances on system dynamic characteristics.

Backlash
In the actual production, the backlash is inevitable. It can not only facilitate lubrication, but also make up for the adverse effects of manufacturing and installation errors on the system transmission. With the change of the difference between the relative meshing vibration displacement along the meshing line and the backlash, that is, the actual meshing displacement along the meshing line is different, the gear pair will have different meshing states. Therefore, taking basic parameters in Table 1 as research object, the backlash varies 2b from 0 to 100 lm, and other parameters of system remain unchanged, as shown in Table 5, to study the influence of the backlash on the system dynamic characteristics. The change curves of the RMS values of the relative comprehensive vibration accelerations of left and right meshing lines of herringbone gear pair before and after 3d modification with the change of 2b are obtained, as shown in Fig. 24. From Fig. 24, with the increase of 2b, the RMS values first increase, then decrease, and then keep increasing, with a peak and a jump at about 17 lm and 23 lm, respectively. In  contrast, the vibration acceleration is decreased after modification. The diagrams of global and local bifurcation and the MLE of dimensionless comprehensive relative vibration displacements of left and right meshing lines of herringbone gear pair with the change of 2b before and after 3d modification are obtained, as shown in Fig. 25-26. As can be seen from Fig. 25-26, with the increase of 2b, the system shows complex bifurcation characteristics when the backlash is small, and finally it keeps increasing in a chaotic state. As follows: In Fig. 25(a), when 2b increases from 0 lm to 3 lm, the system shows multiple periodic motion, and the MLE is less than 0. when 2b increases from 3 to 7 lm, the MLE is greater than 0 and the system motion is chaotic. When at 7 lm, the system bifurcates into the twice periodic motion and the MLE is less than 0. Until 10 lm, it enters chaotic motion and keeps the state. When the MLE is equal to 0, the system shows quasi-periodic motion. In Fig. 25(b), after modification, when 2b increases from 0 to about 4 lm, the system is chaotic motion including a fivefold periodic window and the MLE is greater than 0. At about 4-5 lm, the system bifurcates into multiple periodic motion and the MLE is less than 0. The system is in chaos again at about 5 lm-9 lm, and the MLE is greater than 0. Then the system bifurcates into two times periodic motion at about 9 lm until it enters chaotic motion at 12 lm and remains the state. To sum up, with the increase of backlash 2b, the modification affects system motion state, increases the system chaotic region, and makes the system in a variety of periodic motion states.
Similarly, the local analysis method is used to further analyze the system vibration characteristics shown in Fig. 25 and 26 under the specific parameter after modification. As shown in Fig. 27, when 2b is about 2 lm, the time history is periodic motion, the phase diagram is multiple closed curves, the spectrum diagram is discrete, and there are five attractors in Poincare map diagram, so the system moves in five times period motion. As shown in Fig. 28, when 2b is about 20 lm, the time history diagram has no periodicity, phase diagram fills a closed region irregularly, spectrum diagram has a continuous spectrum band, and Poincare map diagram shows dense points in pieces, so the system motion is chaotic.

Bearing clearance
Considering the actual gear transmission system, due to the different strength and stiffness of gearbox and shaft, the installation error and loaded deformation of box-bearing and bearing-shaft will make the GS have the bearing clearances. The difference between bearing clearance and transverse vibration displacement in the y direction of the established model changes, that is, the actual vibration displacement along the y direction is different, which will also cause the meshing vibration displacement along the meshing line to be different, thus affecting the meshing state of gear pair. Therefore, taking basic parameters in Table 1 as research object, the bearing clearance 2c varies from 0 lm to 100 lm, and the others remain unchanged as shown in Table 6, to study the influence of bearing clearance on the system dynamic characteristics.
The change curves of the RMS values of the relative comprehensive vibration accelerations of the left and right meshing lines of herringbone gear pair before and after 3d modification with the change of 2c are obtained, as shown in Fig. 29. From Fig. 29, with the increase of 2c, the RMS values of acceleration increase continuously, and there are no peak and jump before and after modification. In contrast, the system motion state becomes complicated as 2c increases and the RMS values decrease after modification.
The diagrams of the global and local bifurcation and the MLE of the dimensionless comprehensive relative vibration displacements of left and right meshing lines of herringbone gear pair with the change of 2c before and after 3d modification are obtained, as shown in Fig. 30-31. As can be seen from Fig. 30 and 31, with the increase of 2c, the system shows very complex bifurcation characteristics when 2c is small, finally keeps increasing and is in a chaotic state before and after modification. As follows: In Fig. 30(a), when 2c is about 0 lm-4 lm, the system moves from a single period to a chaotic motion and the MLE is greater than 0. From 4 to 6 lm, the system bifurcates into 2 times periodic motion and the MLE is less than 0. From 6 to 9 lm, it enters chaotic motion again and the MLE is greater than 0. The system is a short 2 times periodic motion from 9 to 10 lm, and from 10 to 50 lm, the system is always in chaotic motion during existing periodic motion and the MLE   is greater than 0. In Fig. 30(b), after modification, the system motion state is relatively regular with the increase of 2c. From 0 to 1.40 lm, the system moves in multiple periods and the MLE is less than 0. From 1.40 to 1.90 lm, the system is the transient chaotic motion. From 1.90 to 3 lm, the system bifurcates into twice periodic motion and then enters chaos. Until about 5 lm, it becomes twice periodic motion, then enters chaos, and then bifurcates into twice periodic motion. Until about 11 lm, the system re-enters chaos and the MLE is greater than 0. From 11 to 50 lm, the system moves from chaotic motion to chaotic motion through twice period and keeps chaotic motion since then.
To sum up, the system motion state is superimposed by periodic motion and chaotic motion before modification. After modification, the system motion state changes alternately in multiple periodic motion and chaotic motion and changes are relatively regular.
Similarly, for Figs. 30 and 31, the local analysis method is used to further analyze system motion characteristics under specific parameter after modification. As shown in Fig. 32, when 2c = 5 lm, the time history diagram is obviously periodic motion, phase diagram is a closed curve, the spectrum diagram is discrete, and Poincare map diagram has two attractors, at which time the system moves twice periodically. As shown in Fig. 33, when 2c = 40 lm, the time history diagram is not periodic, phase diagram is irregular and intertwined, the spectrum diagram has continuous spectral band, and there are dense points in Poincare map diagram, at which time the system motion is chaotic.

Influence of STE on system dynamic characteristics
Considering the actual machining accuracy on gear, the meshing errors of gear pair are generally decomposed into transmission errors consisting of tooth profile deviation and tooth pitch deviation. They mainly refer to the meshing error of gear teeth, which are called static transmission error (STE). The differences between them and relative vibration displacement on the meshing line of gear pair represent the actual vibration displacement, which can be introduced into the dynamic equation. Studying the influence of static transmission error by changing error amplitude on the system dynamic characteristics can find the influence of different errors on the vibration characteristics when the gear teeth are machined in practice, which will provide a theoretical basis for improving the machining accuracy of the gear. Therefore, taking basic parameters in Table 1 as research object, the error amplitude e r varies from 0 to 20 lm, and the others remain unchanged, as shown in Table 7.
The change curves of the RMS values of relative comprehensive vibration accelerations of left and right meshing lines of herringbone gear pair before and after 3d modification with the change of e r are obtained, as shown in Fig. 34. In Fig. 34, with the increase of e r , the RMS values increase continuously and have no peak before and after modification. There has an upward jump at about e r = 2.77 lm, while after modification, there is a jump at about e r = 3.40 lm. In contrast, the RMS values greatly decrease after modification, which indicates that the modification can reduce the impact of the STE on the dynamic characteristics of the GS.
The diagrams of global and local bifurcation and the MLE of dimensionless comprehensive relative vibration displacements of left and right meshing lines of herringbone gear pair with the change of e r before and after 3d modification are obtained, as shown in Figs. 35-36. From Fig. 35, 36, the system relative vibration displacements on the meshing line keep increasing, which shows the complex and changeable bifurcation characteristics, and finally keep increasing in a single cycle motion. As follows: In Fig. 35(a), when e r is at about 0 lm-1.33 lm, the MLE is less than 0 and the system is in multiple periodic motion. At about 1.33 lm-4.5 lm, the MLE is greater than 0, and the system is in chaotic motion. At about 4.5 lm, the system bifurcates into a twice periodic motion, finally enters a chaotic motion, which contains a periodic window, and the MLE is greater than 0. Until at about 9 lm, the system bifurcates into a short twice periodic motion, and then enters a chaotic motion at about 10.95 lm. Finally, it bifurcates into a single cycle motion at about 13 lm and remains unchanged, and the MLE is less than 0. In Fig. 35(b), after modification, when e r is at about 0 lm-4 lm, the MLE is less than 0, and the system moves with multiple cycles. At about 4 lm-5.5 lm, the system occurs chaotic motion, which contains multiple periodic windows. At about 5.5 lm-6 lm, the system bifurcates into briefly periodic motion, and then enters chaos again. Until at 9 lm, the system is periodic motion and the MLE is less than 0. At about 10 lm, the system enters chaotic motion again, and becomes   64 lm, the system enters chaos again, then bifurcates into a single period motion at about 17 lm and keeps unchanged, and the MLE is less than 0. To sum up, the modification can increase the system period window, delay system to enter chaotic motion. The probability of system entering chaotic motion with twice periods is relatively large. Similarly, for Fig. 35 and 36, the local analysis method is used to further analyze the system motion characteristics under specific parameters after modification. As shown in Fig. 37, when e r = 7 lm, the time history diagram is not periodicity, phase diagram has a lot of closed curves, which are irregularly intertwined, spectrum diagram has an obvious continuous spectrum band, and there are patches of dense points in Poincare map diagram, at which time the system is chaotic motion. As shown in Fig. 38, when e r = 19 lm, the time history diagram has periodicity, phase diagram is a closed curve, the spectrum diagram is a discrete spectrum line, and there is an attractor in Poincare map diagram, at which time the system moves in a single period.
6.5 Influence of damping ratio on system dynamic characteristics There are many damping factors such as gear meshing damping, shaft damping and bearing damping in gear transmission system, which are the unique energy consumption parameters of the system. In the study of gear dynamics, the meshing damping between gear pairs is mainly considered, and different meshing damping is usually characterized by different damping ratios f. Therefore, taking basic parameters in Table. 1 as the research object, this paper analyzes the influence of damping ratio f on system dynamic characteristics under various excitations before and after 3d modification, to obtain the importance of damping in the GS and the motion state of the system under modification when the damping changes. Thus, the damping ratio f varies from 0 to 0.2, and the others keep unchanged, as shown in Table 8. The change curves of the RMS values of the relative comprehensive vibration accelerations of left and right meshing lines of herringbone gear pair before and after 3d modification with the change of f are obtained, as shown in Fig. 39. From Fig. 39, with the increase of f, the RMS values of system decrease  continuously, and have no peak and obvious jumping. In contrast, with the increase of f, the system vibration acceleration decreases after modification. The global bifurcation diagrams of the dimensionless relative comprehensive vibration displacements of the left and right meshing lines of herringbone gear pair with the change of damping ratio f before and after modification are obtained, as shown in Fig. 40. As shown in Fig. 40, with the increase of f, the relative vibration displacements of the meshing line of gear pair decrease continuously. In addition, when f is at about 0-0.05, the system is chaotic motion. When f is at about 0.05-0.084, the system bifurcates into twice periodic motion, and then keeps in single periodic motion. After modification, the system is chaotic motion when f is at about 0-0.027. When f is at about 0.027-0.06, the system bifurcates into twice periodic motion, and then keeps single periodic motion. In contrast, with the increase of f, the chaotic region of system becomes larger, but the chaotic time decreases and the periodic motion time increases after modification.
Similarly, for Fig. 40, the local analysis method is used to further analyze system motion characteristics under specific parameters after modification. As shown in Fig. 41, when f is 0.05, the time history diagram has periodicity, phase diagram has two closed rings, spectrum diagram has discreteness, and Poincare map diagram has two attractors, so the system moves twice as periodically. As shown in Fig. 42, when f is 0.1, time history diagram has obvious periodicity, phase diagram is a closed curve, spectrum diagram has a single discrete spectral line, and   Poincare map diagram has an attractor, at which time the system moves in a single period.

Conclusions
The 3d modification method of herringbone gear by form grinding proposed in this paper only removes the tooth surface materials from the microscopic perspective, and can be applied to the power transmission system in the fields of ship, aviation, etc. The improvement of meshing gear teeth from the vibration source can not only improve the meshing characteristics, but also eliminate the system transmission error, improve the bearing capacity and reduce vibration. In addition, compared with rack cutters and other machining gears, this method can effectively realize 3d modification of tooth surface, with high machining precision and efficiency, improving gear wear resistance, having good smoothness and cooling, and easy for practical operation. Because the tooth surface friction is not taken into account in the machining error, it is necessary to consider the relationship between friction and modification to study the effects of different system parameters on tooth surface micromorphology and friction, so as to further understand the role of friction effects in gear nonlinear dynamics under complex excitation with multi-clearances, which is a certain challenge.
In this paper, the nonlinear dynamic characteristics of herringbone gear transmission system under multiclearance and complex excitations before and after 3d modification are studied. The combination of the global vibration characteristics in the parameter domain on the diagrams of bifurcation and the maximum Lyapunov exponent, and the local vibration characteristics under the specific parameters on the diagrams of time history, phase, spectrum and Poincare map can fully understand the dynamic characteristics of herringbone gear transmission system before and after modification under complex excitations, and provide a basis for the selection and optimization of gear system parameters. The research results are as follows: (1) The standard tooth surface contact has tooth side edge contact, the axis non-parallelism error causes the contact path to deviate from pitch circle distribution, and the center distance error affects the length of contact line. The 3d modification can eliminate the tooth edge contact, and reduce the influence of installation error on the contact performance of tooth surface. The LTE amplitude of standard tooth surface is 4.194 00 , while the amplitude is 1.157 00 after optimizing modification. The error fluctuation amplitude is reduced by 72.41% and becomes relatively gentle.
(2) With the increase of input power, input speed, backlash and static transmission error, the system exists the jump, while with the increase of bearing clearance and damping ratio, there is no jump. The 3d modification can achieve the purpose of reducing vibration and noise of system.
(3) Compared with other parameters, the change of input speed and static transmission error makes system have complex and changeable bifurcation characteristics. The 3d modification can reduce the low-speed system vibration, but increase the highspeed system vibration, which has a great influence on the high-speed vibration. With the increase of static transmission error, the modification can increase periodic window of system and delay the system into chaotic motion. The probability of system entering chaotic motion with twice period is relatively large. (4) With the increase of input power, the modification can eliminate the jump, but increase the chaotic interval. With the increase of backlash, the modification increases the chaotic region of system and makes system in a variety of periodic motion states. Before and after modification, the system shows complex bifurcation characteristics when backlash is very small, and finally increases to chaotic motion. With the increase of bearing clearance, the modification makes the system motion alternate between multiple periodic motion and chaotic motion, and makes the system motion more regular. With the increase of damping ratio, the chaotic region of system becomes larger after modification, but chaotic time decreases and the periodic motion time increases.