Study on bifurcation characteristics of multi-clearance bending torsional coupling gear transmission based on EHL

In order to study the influence of tooth surface friction on the non-linear bifurcation characteristics of multi-clearance gear drive system, a 6 degree-of-freedom bending torsional coupled vibration model was established. The time-varying mesh stiffness, backlash, support clearance and damping were considered comprehensively in this non-linear vibration model. Loaded tooth contact analysis was used to calculate the time-varying mesh stiffness. Based on the elasto-hydrodynamic lubrication, the time-varying friction coefficient was calculated. Runge–Kutta numerical method was used to solve the dimensionless dynamic differential equation. Using phase diagram, Poincaré diagram, time history diagram, and spectrum diagram, the influence of tooth surface friction on bifurcation characteristics was studied. The results show that the system undergoes a change from 1-periodic motion, multi-periodic motion, to chaotic motion through bifurcation and catastrophe when the speed changes independently. When the friction coefficient of tooth surface changes from 0, 0.05 to 0.09, the chaotic motion of the system is suppressed. Similarly, with the increase in tooth friction, the chaotic motion characteristics are suppressed. Tooth surface friction is the main factor affecting chaotic motion. With the increase in friction coefficient of tooth surface, the chaos characteristic does not change obviously and the vibration amplitude decreases slightly.


Introduction
Tooth surface friction is one of the key factors that affect the vibration and noise of the gear transmission system. In order to clarify the non-linear dynamic behavior of the gear transmission system under the real working condition, many scholars have done a lot of work on the research of the non-linear dynamic model of the gear with friction. Azar and Crossley 1 considered the time-varying mesh stiffness, tooth surface friction, tooth separation, and gear modification and established the differential equation of motion of spur gear transmission system. Rao 2 explored the coupling characteristics of friction coefficient and backlash and found that the change in backlash will affect the friction coefficient. Yang and Lin 3 established the differential equations of motion of spur gear transmission system including involute profile, energy loss, time-varying mesh stiffness, material elastic modulus, backlash and tooth surface friction. Kahraman and Singh 4 established a 3-degree-of-freedom (DOF) gear system nonlinear model considering the non-linear factors such as support clearance, backlash, and time-varying mesh stiffness. Menon and Krishnamurthy 5 established a dynamic model including tooth surface friction and tooth side clearance and applied neural network compensation method to explore its dynamic impact on the transmission system. Velex and Cahouet 6 established the finite element model of spur gear and helical gear with friction, based on Coulomb friction theory and finite element theory. Howard et al. 7 established the dynamic model of spur gear transmission system with 16 DOFs including the factors affecting the friction of tooth surface. Vaishya and Singh 8-10 established a dynamic model considering friction and simplified the time-varying mesh stiffness to rectangular wave. Wang et al. 11 analyzed the chaos and bifurcation phenomenon of the gear transmission model with 2DOFs including friction. Yuan et al. 12 studied the application of fractal description in the characterization of engineering surface and wear particles. Lundvall et al. 13 used the method of multi-body dynamics to establish a 3-DOF motion differential equation including friction. He et al. 14 established a multi-DOF dynamic model including friction of tooth surface based on Coulomb friction theory. Liu and Parker 15 established the coupling vibration mechanical model of a single pair of gear system including friction, stiffness, and damping. Tang et al. 16 established a single pair of gear transmission system dynamics model considering friction, clearance, stiffness, and other factors. Zhang et al. 17 established a 6-DOF dynamic model including friction and solved the influence of tooth surface roughness on friction by numerical method. Lu et al. 18,19 use the normal distribution function to describe the random backlash and analyze the dynamic characteristics of single-pair gear transmission system. Wen et al. 20 studied the response of a single-DOF gear system with random backlash under white noise excitation. Wang 21 established the dynamic model of herringbone gear including tooth friction and analyzed the influence of tooth friction on chaos and bifurcation of the system. Ma et al. 22 have studied the research progress of the dynamics of the gear transmission system with cracks. Bao et al. 23 studied the dynamic characteristics of variable pitch external gear system with time-varying friction coefficient. Zhou et al. 24 studied the non-linear dynamic characteristics of the coupling gear rotor bearing system under the internal and external excitation. Lu et al. 25 studied the non-linear vibration characteristics of the gear rotor bearing system. Wang et al. 26 studied the establishment method of time-varying mesh stiffness model of helical gear pair considering axial meshing force. Yin et al. 27 studied the dynamic modeling and analysis of singlestage gear system with tooth impact and friction. Li et al. 28,29 studied the influence of the backlash with fractal characteristics on the dynamic characteristics of the gear bearing system. Chen et al. 30,31 introduced the fractal theory and studied the influence of the gear backlash with fractal characteristics on the dynamic characteristics of single-stage gear system. Shi and Zhao 32 studied the vibration characteristics of the second-stage helical cylindrical gear rotor system and believed that the high-frequency modal frequency of the gear rotor system changed greatly with the rotating speed due to the existence of gyroscopic effect. Zhang et al. 33 studied the stiffness calculation and influence factor analysis of spur gear pair considering timevarying friction.
From the above studies, some research has been done on the non-linear characteristics of gear transmission systems with considering the friction of the tooth surface. However, in most of these studies, the actual meshing changes in the tooth surface loaded contact process are less considered. Therefore, in order to study the influence of tooth surface friction on the bifurcation characteristics of the system, through considering the time-varying mesh stiffness, time-varying friction coefficient, fractal backlash, support clearance and damping, and a 6-DOF multi-clearance bending torsional coupling vibration model of the gear transmission system are established. The loaded tooth contact analysis (LTCA) is used to calculate the time-varying mesh stiffness. The time-varying friction coefficient is calculated based on the elasto-hydrodynamic lubrication (EHL). The fractal theory is introduced into the dynamic model, and the W-M function is used to describe the change trend of the backlash. The Runge-Kutta numerical method is used to solve the dimensionless dynamic model, and the phase diagram, Poincare´diagram, and bifurcation diagram of the system response are obtained. The influence of tooth friction on the bifurcation characteristics of the system is analyzed when the parameters change. The numerical method is used to solve the problem, and the bifurcation characteristics of the system are studied under the influence of tooth friction.
In Figure 1, the angular displacement of driving gear and driven gear is represented by u s and u p , respectively. The lateral displacement is represented by X s and X p , respectively. The longitudinal displacement is represented by Y s and Y p , respectively. The base circle radius is represented by r bs and r bp , respectively. The stiffness, half support clearance, and damping in the x and y directions of the driving gear are represented by K sx , K sy , b sx , b sy , C sx , and C sy , respectively. The stiffness, half support clearance, and damping in the x and y directions of the driven gear are represented by K px , K py , b px , b py , C px , and C py , respectively. Here, i is x and y. The comprehensive meshing error, time-varying mesh stiffness, half backlash, meshing damping, and friction coefficient of the gear pair are expressed by e(t), K sp (t), b h , C sp , and m, respectively. The friction arm of driving gear and driven gear is represented by L s and L p , respectively. The eccentricity errors of driving gear and driven gear are expressed by e s and e p , respectively. The initial phase angles of eccentricity errors are a s and a p , respectively. The phase angle between the driven gear and the driving gear is d. The linear displacements x s and x p produced by the angular displacements of the driving gear and the driven gear on the meshing line can be shown in equation (1) x s = u s r bs The comprehensive meshing error includes the static transmission error and eccentricity error of driving gear and driven gear. The harmonic period expansion of meshing transmission error by Fourier function can be expressed as follows In equation (2), E is the average value of static transmission error. a s and a p are the initial phase angles of the bth harmonic component, respectively. l is the number of harmonic terms expanded by Fourier. e s and e p are the eccentricity error. v s and v p are rotation frequency of driving gear and driven gear (rad/s). The non-linear functions of backlash and support backlash can be expressed in equation The dynamic meshing force of gear pair is defined as the sum of elastic restoring force and damping force, which can be expressed as follows In equation (4), X is the relative displacement of the gear pair's meshing point.
The damping coefficient is expressed as follows In equation (5), j is the relative damping ratio of the meshing pair, m s is the mass of the driving gear, m p is the mass of the driven gear, and the relative damping ratio of the gear transmission system is generally 0.03 to 0.17.

Establishment of dynamic backlash b h (t) based on fractal theory
The W-M model of gear pair meshing based on the fractal theory is shown in Figure 2. A is the highest point of surface profile of gear 1. B is the highest point of surface profile of gear 2. The research in Yin et al. 27 shows that the machined surface has fractal characteristics, and the tooth surface is also machined surface. The fractal formula of backlash can be defined by W-M function as follows b h (t) = X + ' In equation (6), b h (t) is the backlash varying with the time. l b is the characteristic scale coefficient. D is the fractal dimension, and only two-dimensional fractal is considered. k is the contact coefficient. t is time.
Using equation (6), the simulation of backlash under different fractal dimensions is shown in Figure 3. Here, l is 1.5 and D is 1.1-1.9, respectively. The larger the fractal dimension D, the more complex the micro morphology of W-M simulation. The curve changes more dramatically with the increase in fractal dimension D. The larger the fractal dimension D, the smoother the tooth surface and the smaller the backlash. The average backlash b hm with different fractal dimension D is shown in Table 1.

Calculation of time-varying tooth surface friction based on EHL
The schematic diagram of gear force analysis is shown in Figure 4.
In the process of gear transmission, the tooth surface will be affected by the dynamic meshing force F N along the direction of the meshing line and the friction force F f perpendicular to the meshing line. XY is the theoretical meshing line length. AD is the actual meshing line length. r as and r ap are the top circle radius of driving gear and driven gear. r s and r p are the pitch circle radius of driving gear and driven gear. v s and v p are the rotation frequency of driving gear and driven gear, rad/s. a n is the meshing angle.
From the geometric relationship of gear transmission, we can get the length XY of theoretical meshing line, the length AD of actual meshing line, the distance XA between the start end of actual and theoretical meshing line, and the distance YD between the end of actual and theoretical meshing line. The calculation expression is as follows XY = r s + r p À Á sin a n AD = ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi r 2 as À r 2 bs q + ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi r 2 ap À r 2 bp q À r s + r p À Á sin a n XA = r s + r p À Á sin a n À ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi r 2 ap À r 2 bp q YD = r s + r p À Á sin a n À ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi r 2 as À r 2 In equation (7), P b is the pitch of gear. Z s and Z p are the number of teeth of driving gear and driven gear, respectively. Therefore, the expression of the friction arm can be derived as follows In equation (8), L s and L p are the time-varying friction arms of the driving gear and driven gear, respectively, mod (Á) is the complementary function. The relative velocity u e of two meshing contact tooth surfaces at any position is shown in equation (9) Since the direction of friction force will change at the node, the direction coefficient l can be used to represent the direction change. The direction coefficient is defined as In equation (10), sgn(Á) is a symbolic function.
The calculation formula of Coulomb's law friction model is as follows In equation (11), m 0 is the amplitude of friction coefficient.
The time-varying sliding friction force of tooth surface according to Coulomb's law can be expressed as follows In equation (12), m is the friction coefficient.
In aviation standard, the empirical formula of friction coefficient at any meshing point K synthesized according to the actual test results is as follows In equation (13), n 0 is the initial viscosity of the lubricating oil (Pa/s). W t is the tangential load per unit tooth width (N/mm). u sp is the entrainment velocity (m/s). R is the comprehensive curvature radius of the gear pair (mm). R a is the mean value of tooth surface roughness (mm). Ignoring the effect of roughness, R a is equal to zero. The friction coefficient is related to the geometric shape of the teeth, the surface roughness of the teeth, the relative sliding velocity of the teeth, the contact stress between the teeth, and the lubrication of the system, so it is difficult to accurately predict the friction coefficient. When the lubrication effect of tooth surface is considered, the gear contact problem becomes the linear contact EHL problem. Therefore, the friction coefficient is calculated by EHL theory.
By solving Reynolds equation, viscous pressure equation, dense pressure equation, energy balance equation, and load equation, the oil film pressure distribution and film thickness distribution are obtained. Combined with Ree-Eyring fluid equation, the oil film shear stress distribution is calculated, and then, the friction force and friction coefficient between two contact surfaces are obtained. It is assumed that the tooth surface is completely smooth, that is, the gear is in the state of complete EHL. In the lubrication analysis, the influence of temperature effect is ignored, and it is considered that the oil film thickness has not changed.
The Reynolds equation is In equation (14), x is the tangent direction of tooth surface. r is the density of lubricating oil (kg/m 3 ). h is the oil film thickness (m). p is oil film pressure (Pa). n is the viscosity of Ree-Eyring fluid (PaÁs). u sp is the entrainment velocity (m/s), u s , u p are the tangential velocity of lubricating oil along gear axis (m/s). x in and x out are the entrance and exit positions of the calculation area in the x-direction. r* and (r/n) e are the equivalent parameters. The equivalent parameters are defined in equation (15) r n e = 12 n e r 0 e n 0 e À r 00 e r Ã = ½r 0 e n e (u s À u p ) + r e u s =u sp The boundary conditions for the differential equations are p(x in ) = 0, p(x out ) = 0, p(x in ) . 0, and (x in \ x \ x out ). n* is the equivalent viscosity of lubricating oil (PaÁs).
The film thickness equation is as follows In equation (16), h 0 is the central film thickness (m) of the two gears without meshing deformation, which can be determined by the load balance condition. R is the comprehensive curvature radius (m). d(x) is the elastic deformation displacement (m) due to the contact pressure of the tooth surface. E# is the comprehensive elastic modulus (Pa) of the two gears. Sr is the sliding roll ratio when the gear is engaged. C d is the viscosity pressure coefficient of Cameron's viscosity pressure formula (m 2 /N), C d = x/15, and x is the viscosity pressure coefficient of Barus' viscosity pressure formula (m 2 /N). The viscous pressure equation and the dense pressure equation can be expressed as follows n = n 0 expf( ln n 0 + 9:67)½À1 + (1 + 5:1 3 10 À9 p) z 0 g r = r 0 1 + 0:6 3 10 À9 p 1 + 1:7 3 10 À9 p ! z 0 = x=½5:1 3 10 À9 3 ( ln n 0 + 9:67) In equation (17), n 0 is the initial viscosity of the lubricating oil (PaÁs). r 0 is the initial density of the lubricant (kg/m 3 ).
The load balance equation can be expressed as follows x 0 = pm n cos a n (e À 1) In equation (18), W is the equivalent cylinder external load (N). c w (t) is the load factor of gear. Taking the meshing line of gear transmission as the coordinate axis and the node as the origin, x in and x out are the coordinates of entering and leaving the meshing point, respectively, m n is the modulus and e is the gear coincidence degree.
The constitutive equation of Ree-Eyring fluid can be expressed as follows In equation (19), t is the shear stress of oil film in x-direction (Pa). u is the velocity of lubricating oil in x-direction (m/s). t 0 is the shear stress of Ree-Eyring fluid (Pa). The formula of velocity and velocity gradient derived from the motion equation of fluid particle can be expressed in equation (20) The friction force on the contact tooth surface of the gear and the velocity distribution of the lubricant are shown in Figure 5.
Here, the sliding friction is mainly generated in the contact area (2b, b). F 1 and F 2 are the friction force on the surface of the two tooth surfaces, respectively. The sliding friction coefficient of the gear based on the EHL model can be expressed as follows In equation (21), Àb tj z = h dx. The above mathematical model is dimensionless. Then, the multi-grid method is used to solve the pressure, and the numerical method is used to solve the friction coefficient and friction.
In order to verify the correctness of the proposed sliding friction coefficient, this article gives the friction coefficient calculated by EHL empirical formula in Xu et al. 34 and makes a comparison. Calculation model expression of friction coefficient based on EHL empirical formula proposed by Xu et al. 34 is as follows In equation (22), P h is the maximum contact Hertz stress (GPa). R rms is the root mean square (RMS) value of roughness. f (Sr, P h , n m , R rms ) is calculated in equation (23) as follows f (Sr, P h , n m , R rms ) = b 1 + b 4 Sr j jP h log 10 (n m ) + b 5 e À sr j jP h log 10 (n m ) + b 9 e R rms ð23Þ The values of b 1 , b 2 , ..., b 9 are shown in Table 2.
Maximum Hertz contact stress can be expressed as In equation (24), E# is the comprehensive elasticity modulus, GPa, and calculated as follows In equation (25), m s and m p are Poisson's ratio of driving gear and driven gear, respectively. E s and E p are the elastic modulus of driving gear and driven gear, respectively.
The instantaneous velocity u s and u p at any meshing point is calculated as follows In equation (26), r s and r p are the curvature radius. s is the distance between the instantaneous meshing point and the actual meshing starting point, m.
The sliding roll ratio Sr, the tooth surface relative sliding velocity u e , the entrainment velocity u sp , and the rolling velocity u r are calculated in equation (27) Here, regardless of the influence of roughness, RMS value of roughness R rms = 0.0 mm. Mobil dte24 antiwear lubricating oil is selected, and its relevant parameters are shown in Table 3.
In order to compare the calculation results of different friction coefficient models, the power is 100 kW and the input speed n p is 5000 r/min. The change curve of tooth surface friction coefficient with slip roll ratio sr is shown in Figure 6. Figure 6(a) and (b) shows the curve of sliding friction coefficient with the entrainment velocity at 1.0, 2.0, and 3.0 m/s. EHL model comprehensively considers the load distribution, relative sliding velocity of tooth surface, gear surface morphology, and lubrication conditions. In this article, the friction coefficient calculation model based on EHL is used. From Figure 6, due to the calculation of EHL model without considering surface roughness, and assuming that the lubrication state is complete EHL, the value of sliding friction coefficient obtained by EHL model is smaller than that of Coulomb model. The coefficient of sliding friction calculated by EHL model increases first and then decreases with the decrease in sliding roll ratio Sr. The coefficient of sliding friction decreases with the increase in entrainment velocity u sp .

Calculation of time-varying mesh stiffness based on LTCA
LTCA method is used to calculate the time-varying mesh stiffness, and the results are compared with those of Ishikawa method. LTCA model of driving gear and driven gear is shown in Figure 7.
Under the load P, the driving gear goes through an approach Z. Due to the tooth deformation, contact load will become distributed. After contact deformation, the state can be described by the following equation for the tooth pair k Table 2. Parameters of EHL friction coefficient calculation model.   The elastic deformation of the gear teeth under normal load, including the geometric bending deformation d B of the gear teeth, the elastic deformation d M at the root of the gear teeth, and the contact deformation d C at the meshing point, is obtained. The relationship between the comprehensive elastic deformation Du sp [T sp (k)] and the transmission torque T sp (k) of the gear pair is as follows In equation (29), a, b, and c are the constant terms in the formula, which are determined by the geometric parameters of the gear and the coordinates of the meshing position, T sp (k) (k = 1,2, ...,5) is the torque of the kth meshing position in a meshing cycle. d B = a, d M = bT sp (k), and d C = cT sp (k) 2/3 . Through the finite element calculation and the gear geometric contact simulation analysis program, when the gear pair is in the meshing position k, solving equation (27), the load transmission error Dd T (0.1T sp (k)), Dd T (0.5T sp (k)), and Dd T (0.9T sp (k)) are obtained. The results are brought into equation (28), and the coefficients a, b, and c are obtained, and the function relation between the loaded deformation of the driving gear and driven gear and the nominal load is obtained. The timevarying mesh stiffness of driving gear and driven gear is obtained. The following equation (30) is satisfied with the integrated time-varying mesh stiffness K sp of driving gear and driven gear with the corresponding angular error Du sp (T sp (k)) The discrete value of the time-varying mesh stiffness calculated by LTCA is expanded into a harmonic periodic function with v as the fundamental frequency by  Fourier series. Its mathematical form can be described as follows In equation (31), K m is the average value of the fluctuation of the meshing stiffness. v is the meshing frequency of the gear pair, and the meshing period is T=2p/v. l is the order of the harmonic term expanded by the Fourier series. K a is the stiffness amplitude of the ath (a = 1, 2, 3, ..., l) harmonic component. c a is the initial phase angle of meshing of the ath (a = 1, 2, 3, ..., l) harmonic component. The calculation formula of Ishikawa method is that the contact gear is equivalent to the cantilever beam composed of rectangle and trapezoid to calculate the deformation, and the dangerous section of the gear is defined as the length of rectangle, the specific calculation will not be described here. Figure 8(a) shows the comparison curve between LTCA calculation results and Ishikawa calculation results. The time-varying mesh stiffness including the friction factor of tooth surface is calculated, as shown in Figure 8(b). The time-varying mesh stiffness of the gear considering the friction of the tooth surface is larger than that without considering the friction, and the variation range is about 8%. This is because when the friction factor of tooth surface is considered, the friction makes the meshing stiffness of tooth pair increase in the meshing stage and decrease in the meshing out stage. The larger the friction coefficient is, the greater the variation of stiffness will be.

Establishment of differential equation of bending torsional coupling vibration based on LPM
Based on the Lumped Parameter Theory (LPM), a multi-clearance coupled bending torsional vibration model is established, and the differential equations of the system are obtained as follows (equation (32)). Here, it is assumed that the transmission shaft has enough stiffness, ignoring the torsional stiffness and torsional damping of the shaft, and ignoring the elastic deformation of the shaft In equation (32), u = d À a n , u is the projection angle. I s is the moment inertia of driving gear. I p is the moment inertia of driven gear. T in is the input torque. T out is the output torque. In order to facilitate comparison and analysis, making the solution of the system in the same order of magnitude, dimension time t = v n t and displacement nominal scale b c are introduced.
Here, v n = ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi K m (I s r 2 p + I p r 2 s )=(I s I p ) q . The dimensional displacement, velocity, and acceleration are X = X b c , _ X = _ X v n b c , and € X = € X v 2 n b c , respectively. The dimensional backlash, dimensional rotational speed, dimensionless angular velocity can be expressed as The dimensional comprehensive meshing error of gear pair can be expressed as follows (equation (33)) The non-linear function of dimension clearance of gear pair can be expressed as follows (equation (34)) The dimensionless differential equation of coupled bending torsional vibration of the system can be expressed as follows (equation (35)) In equation (37), m es = I s =r 2 bs , m eq = I p =r 2 bp , F in = T in =r bs , F out = T out =r bp .

Study on the bifurcation characteristics of the system
The parameters of the gear are shown in Table 4. The four-to five-order adaptive variable step Runge-Kutta numerical method is used to solve the problem. The bifurcation diagrams under different bifurcation parameters are obtained and the bifurcation characteristics of the non-linear vibration of the system are analyzed.
Bifurcation characteristics of the system with the change in dimension speed In this system, the given dynamic backlash b h = 40 mm and the support clearance b sx = b sy = 10 mm and b px = b py = 10 mm. When the damping ratio is 0.05, the friction coefficients are 0.0, 0.05, and 0.09, respectively.
The bifurcation characteristics with the change in dimension speed O are shown in Figure 9.
From Figure 9(a), when O is in the [0.58,0.79] range, the system is one-periodic motion; when O is in the [0.79,0.88] range, the motion state of the system shows jumping phenomenon; when O increases to 1.63, the system bifurcates into two-periodic motion; when O increases to 1.72, the system bifurcates into threeperiodic motion and multiple periodic motion; when O is in the [1.81,2.84] range, the system enters the chaotic motion state, and A1 is used to represent the chaotic region; when O continues to increase, the chaotic state of the system gradually weakens and the system enters the periodic state again.
From Figure 9(b), when O is in the [0.58,0.80] range, the system is one-periodic motion; when O increases to 0.81, the motion state of the system shows jumping phenomenon; when O increases to 1.63, the system bifurcates into two-periodic motion; when O increases to 1.83, the system bifurcates into three-periodic motion; when O is in the [1.99,2.62] range, the system enters the Table 4. System parameters.

Parameters
Driving gear Driven gear Module m n (mm) 3.5 3.5 Pressure angle a n (°)  Figure 10. Continued chaotic motion state, and A2 is used to represent the chaotic region; when O continues to increase, the system moves from chaotic state to 2-periodic state again and then finally stabilizes to 1-periodic state. From Figure 9(c), when O is in [0.58,0.91], the system is one-periodic motion; when O increases to 0.92, the motion state of the system shows jumping phenomenon; when O increases to 1.68, the system bifurcates into two-periodic motion; when O is increased to 1.91, the two-periodic motion state of the system shows jumping phenomenon; when O increases to 2.12, the two-periodic motion of the system enters into threeperiodic motion; when O is in the [2.42,2.69] range, the system will enter the chaotic motion state; when O continues to increase, the system moves from chaotic state to two-periodic state again, and then, it stabilizes to one-periodic state when O is equal to 2.81.
From Figure 9(a)-(c), with the increase in friction coefficient, the chaotic region gradually decreases, and the system exits from chaotic motion and returns to single periodic motion when the speed O = 2.62. When the speed is low (O \ 1.8), the system has a stable motion state, and with the continuous increase in rotation speed, the chaotic characteristics of the system are obvious when there is no friction (m = 0.0). With the increase in friction coefficient of tooth surface, the system drives the bifurcation inside the chaotic region by the bifurcation outside the chaotic region, and there is a short period of multiple periodic motions before entering the chaotic region. The increase in friction coefficient will increase the tooth surface resistance of the gear pair, and the area of the chaotic region will be significantly reduced. Figure 10 shows the phase diagram, Poincare´diagram, time history diagram, and spectrum diagram corresponding to the bifurcation characteristics at different rotating speeds when the friction coefficient is 0.05. The whole process of the system from one-period motion to chaos motion and then to one-period motion is analyzed. Here, the red solid dot represents Poincared iagram. Figure 10(a) shows that the system is in stable one-periodic motion. Figure 10(b) shows that when dimension speed O = 1.79, the frequency of bending vibration includes not only meshing frequency but also other continuous frequency components and the system enters aperiodic motion state. Figure 10(c) shows that when O = 1.96, the amplitude of the low-frequency region increases, indicating that the system is in an aperiodic motion state and tends to enter chaotic motion. In Figure 10(d), the ordered structure inside the phase diagram shows that the system is in quasi- periodic motion. Figure 10(e) shows that the sideband component appears near the meshing frequency, and the system enters into the chaotic motion state.
From Figure 10(a) to (e), the system moves from one-periodic motion to chaotic motion after multiple periods, and the phase diagram changes from regular circle to irregular shape and the Poincare´diagram gradually increases from the single vibration period point to the discrete multi-vibration period change. In the low-frequency part of fast-Fourier transform (FFT) spectrum, the amplitude of the system increases continuously, and the dynamic characteristics of the system become worse. Figure 10(e) and (f) shows that the system returns from chaotic motion to stable one-periodic motion state, the low-frequency amplitude of the system in FFT spectrum diagram decreases, and in Poincare´diagram, the range of vibration period number decreases gradually, and the dynamic characteristics of system are improved. The change in dynamic characteristics in Figure 10 is consistent with the performance of dynamic bifurcation characteristics in Figure 9.
Bifurcation characteristics of the system with the change in dynamic backlash Take the dimension speed O equal to 2.4, the support clearance b sx = b sy = 10 mm, b px = b py = 10 mm, and the relative damping coefficient z = 0.05, the friction coefficient is 0.0, 0.05, and 0.09, respectively. The bifurcation diagram of the system dimension-relative displacement with the change in the backlash b h is shown in Figure 11.
From Figure 11(a), when b h is in [0,3.72] range, the system is one-periodic motion; when b h is equal to 3.73, the system bifurcates into two-periodic motion; when b h is in [12.8,13.3] range, the system has multiple periodic motion and two-periodic motion; when b h is increased to 13.4, the system enters into chaotic motion.
From Figure 11(b), when b h is in [0,3.61] range, the system is one-periodic; when b h is equal to 3.62, the system bifurcates into two-periodic motion; when b h continues to increase, the system experiences multiple periodic motion and then enters into chaotic motion. From Figure 11(c), when b h is in [0,3.60] range, the system is one-periodic motion; when b h is equal to 3.61, the system bifurcates into two-periodic motion; when b h continues to increase, the system enters multiple periodic motion.
From Figure 11(a) to (b), the system moves from single-periodic motion to double-periodic motion and then enters chaotic motion through rapid change. The critical backlash does not change with the increase in friction coefficient. When the backlash is less than 4 mm, the system moves in single-periodic motion. When the backlash increases to 11 mm, the system suddenly enters chaotic motion. From Figure 11(b) and (c), the vibration amplitude changes little with the increase in backlash. With the increase in backlash, the system will produce irregular vibration due to the excessive backlash. The vibration amplitude will continuously increase with no tooth surface friction. With the increase in tooth surface friction, the chaotic region is divided from one block to three blocks. In order to avoid the chaotic movement of the system, the value of backlash should be in the range of 0-11.6 mm.
When the backlash changes from 2 to 16 mm, Figure  12 shows that the system moves from one-periodic motion to two-periodic motion and moves from the irregular vibration to chaos. From Figure 12(a) to (c), it is shown that in the low-frequency part of the FFT spectrum, the amplitude increases continuously, and the dynamic characteristics of the gear pair become worse. When the backlash b h =16 mm, the FFT spectrum shows continuous other frequency components. Under the action of an unstable attractor, the system moves from period motion to chaos motion. It is also proved that the two-periodic motion is a way to enter chaos. When the dimension speed O = 2.4, the relative damping ratio is 0.05, 0.07, and 0.09, respectively. The bifurcation diagram of the system's dimension relative displacement changing with the friction coefficient m is shown in Figure 13.
From Figure 13(a), when the friction coefficient m is in [0,0.43] range, the system is chaotic motion; when m is in [0.43,0.69] range, the system includes two chaotic bands, which are accompanied by periodic motion and quasi-periodic motion; when the friction coefficient continues to increase, the system is stable in chaotic motion. From Figure 13(b), the main motion of the system is chaos motion. When the friction coefficient increases to 0.85, the system will be accompanied by periodic motion and quasi-periodic motion, and the state of chaos motion will be weakened. From Figure  13(c), when the friction coefficient is equal to 0.55, the system gradually changes from multi-period to threeperiod and finally stabilizes at two-period.
When the system dimension-rotational speed O = 2.4 and the relative damping ratio is 0.05, and the friction coefficient of the tooth surface changes from 0.46 to 0.69, the whole process of the system changing between chaotic motion and period doubling motion is shown in Figure 14. Figure 14(a) shows that the system is in chaotic motion state, and Figure 14(b) shows that under the action of stable attractor, the system changes from chaotic motion to periodic motion with the increase in friction coefficient of tooth surface and then   enters into chaotic motion in the form of cataclysm. As shown in Figure 14(c), as the friction coefficient of tooth surface continues to increase, the system changes from the same cataclysm form and finally returns to the chaotic motion state. In the low-frequency part of FFT spectrum, with the primary and secondary changes in the three states of hybrid motion, periodic motion, and chaotic motion, the amplitude first increases and then decreases, and the dynamic characteristics of gear pair are in good and bad states.

Conclusion
In this article, a non-linear dynamic model of 6-DOF bending torsion coupling gear transmission system is established, which considers the time-varying friction coefficient of tooth surface, the time-varying mesh stiffness, and the fractal backlash. The bifurcation characteristics of the system are analyzed and studied under the influence of the parameters such as dimension speed, backlash, and friction coefficient of the tooth surface. The conclusions are as follows: 1. The innovation of this article is to introduce EHL theory and LTCA technology into the establishment of bending-torsion coupling nonlinear vibration model of transmission system. The EHL model comprehensively consider the load distribution, the relative sliding velocity and rolling velocity of the tooth surface, the surface morphology of the gear, and the influence of lubrication conditions. EHL and LTCA methods are used to accurately obtain the friction coefficient and meshing stiffness of the time-varying tooth surface, which reflect the real meshing characteristics of the tooth surface from a micro-perspective and more accurately reflect the influence of various factors on the non-linear characteristics of the system. By using the fractal theory, dynamic backlash more accurately reflects the influence of the real tooth surface parameters on the non-linear bifurcation characteristics. 2. When the rotational speed changes independently, the chaos is obvious without the friction of the tooth surface. With the increase in tooth surface friction coefficient, the bifurcation behavior outside the chaotic region of the system will cause the bifurcation behavior inside the chaotic region, and the chaotic behavior of the system will be suppressed. 3. When the backlash changes independently, the vibration amplitude of the system increases continuously without tooth surface friction, and the non-linear characteristics of the system are obvious. With the increase in friction coefficient of tooth surface, the system gradually changes from chaotic behavior to quasi-periodic behavior and then to double-periodic behavior, which is suppressed. In order to avoid the chaotic movement of the system, it is very important to ensure the rationality of the backlash value. The value of backlash should be in the range of 0-11.6 mm. 4. When the tooth friction changes independently, with the increase in the tooth friction, the chaotic behavior of the system does not change significantly, and the vibration amplitude of the system decreases slightly. With the increase in damping coefficient, the tendency of regional partition of chaos behavior and the effect of chaos suppression are more significant.

Declaration of conflicting interests
The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.