A New Model for Studying Nonlinear Vibration Behaviors of Gear-Bearing System

: A new model for nonlinear vibration behaviors of gear-bearing system is proposed in this work. For presenting the nonlinear excitation from bearing compliance, the enhanced bearing force excitation model containing two kinds of bearing stiffnesses, which are mean stiffness for the load transfer capacity and alternating stiffness for the disturbance resisting ability, is developed. Considering other dynamic excitations including mesh stiffness, contact pressure angle, center distance, unbalance force caused by static and dynamic eccentricity, an advanced iterative numerical method is introduced, which can timely and accurately update the excitations caused by load-dependent and time-varying nonlinearities inside of the system. The constant bearing stiffness and time-varying bearing alternating stiffness models are introduced and compared with the enhanced bearing excitation force model. The parametric resonant regions and system nonlinear periodic motion states are studied and compared for different bearing supporting models. The effects from internal and external excitations on the system nonlinear vibration behaviors are investigated.


Introduction
Gear-bearing transmission system plays a vital role in the rotating machine applications.The investigations for the vibration behaviors of gear-bearing rotating system have been developed step by step and keep attracting attention for the past decades.From the very beginning, Özgüven H. and Houser D. [1] developed a single DOF nonlinear model for gear mesh pair which considers the time-varying mesh stiffness and damping, other forms of errors including pitch, profile and runout errors, backlash and tooth profile modifications.Kahraman A. and Singh R. [2] compared the results from their analytical model of spur gear pair with experimental data.While, the models they provided or used above only had the loaded static transmission error for presenting the nonlinearity of mesh stiffness as system excitation.Then, Kahraman A. and Singh R. [3] considered the bearing compliances and clearances by introducing the piecewise linear functions for the gear systems; the parametric studies including the effect from bearing stiffness to the mesh stiffness ratio on the system dynamic behaviors were conducted.Very soon after that they [4] used the digital simulation and multiple scales method to obtained the frequency response properties of the gear rotor bearing systems; the interaction between the mesh stiffness with the bearing nonlinearity was found to be weak based on their model, but only the primary resonance was considered.Özgüven H. [5] proposed a semi-definite model for the dynamic analysis of spur gear system and the results highlighted that the effects from the bearing compliances and shaft on the system dynamic behaviors were distinct when the mesh pair torsional mode was coupled with other modes governed by bearing compliances and shaft.A finite element model for the gear bearing rotor system was built by Kahraman A. and his cooperators [6]; the unbalance mass, geometric eccentricity were added into their model and the results showed that the bearing elasticity dramatically influenced the vibration behaviors.A number of experiments for the gear system with clearances and time-varying parameters were performed by Kahraman A. and Blankenship G. [7] and abundantly nonlinear phenomena such as amplitude jump, sub and super harmonic resonances, and coexisting stable , sub-harmonic and chaotic motions were observed.Jianjun W. and his colleagues [8] comprehensively reviewed the investigations on the nonlinear vibration for the gear transmission systems for the past two decades since that time; the mainstream on gear-bearing system dynamic modeling from those investigations was taking the multiple clearances and the time-varying mesh stiffness into account; for the bearing part, most of studies only considered the isotropic constant bearing stiffness and clearances through a gap function.
Since the recent ten years, the effects of load-dependent and time-varying characteristics of bearing on gear-bearing system dynamic behaviors drew attention.Baguet S. and Velex P. [9] studied the nonlinear behaviors of the gear-shaft-bearing transmission system; the vibration response of bearing forces was found to be connected with the mesh excitations.In [10] Baguet S. and his cooperator confirmed the importance of considerations of bearing nonlinear behaviors for the gear-bearing system, especially for the cases with larger dynamic forces.Cai-Wan C. [11] presented abundantly nonlinear phenomena of gear-bearing system considering geometric eccentricity, nonlinear mesh forces and nonlinear bearing forces, and with that he [12] emphasized the influences of operation conditions, especially for rotating speed, on the strongly nonlinear vibration behaviors of gear-bearing systems.Zhenxing L. and his/her colleagues [13] proposed a multiple-vibration-parameters-coupling model for learning the nonlinear behaviors of highly nonlinear gear-bearing systems.However, those strong nonlinear vibration behaviors were all from the journal-bearing-gear system.For the rolling bearings, the researches aiming at the nonlinearity of bearing force and time-dependent bearing properties were very limited.The piecewise functions for presenting the relationship between the bearing force with displacement were still introduced by [14,15], but for more complicated configurations such as planetary gear set or two-stage gear transmission.Miroslav B. and Vladimír Z. [16] built a very detailed bearing model for the gear-bearing system dynamic analysis; their research showed the nonlinearity of bearing stiffness caused by varied number of rolling elements in contact; however, the nonlinear characteristics of bearing was presented by complex formulas with a assumptive rolling elements in contact.The comprehensive bearing stiffness matrix techniques computed by the particle differential method were introduced in [17,18], thus the anisotropy property of bearing supporting performance was taking into account.Junyi Y. and Teik L. [17] studied the influences of the interaction between the time varying bearing stiffness with the time-varying bearing stiffness on the dynamic behaviors of the gear-bearing system; they discovered that the bearing loads were significantly affected by the bearing stiffness, while dynamic mesh force was not very sensitive to the fluctuation of bearing stiffness.More parametric studies on bearing stiffness from a coupled gear-bearing system were showed in [18] and they found that the coupling effects between the system nonlinearities deeply influenced the sensitivity of gear dynamic behaviors to the bearing clearance.Tristan E. and Robert P. [19] confirmed the anisotropy of bearing stiffness through experiment results and finite-element/contact mechanics combined method; additionally, the bearing stiffness was found to be load-dependent and could have impact on gear mode.Zehua H., Jinyuan T. and Siyu C. [20] brought the comprehensive bearing stiffness model into a lateral-torsional coupled gear transmission system to study the influence of gyroscopic effect on the vibration behaviors of the system.Recently, Guanghui L., Jun H. and Robert P. [21] investigated the influence of time-varying rolling element bearing stiffness on the system vibration very carefully; strong nonlinear vibration behaviors caused by fluctuation of bearing stiffness over ball pass cycle was observed and theoretical explanation based on perturbation for the system instability caused by the mesh and bearing stiffness excitations was given.
Based on the above descriptions, the effects from the bearing nonlinearity on the vibration behaviors of gear-bearing transmission system are becoming more prominent because of the rapid development of high-speed and heavy-duty rotary machine.The rolling bearing nonlinearity contains both time-dependent nature and load-dependency, while all the researches presented recently only can consider one of them in their model.Also ,the gear mesh models described above neglect the influence of the variation of the contact pressure angle and center distance during the dynamic rotation process, which is proven to be influential in system dynamic behaviors [22,23].Additionally, the geometric eccentricity is verified to be no negligible [9,10,13,22], and the dynamic eccentricity caused by the translational displacements due to coupling effects of gear elasticity and bearing compliance is always neglected.Overall, it is necessary to establish an advanced gear-bearing model, for comprehensively study the coupling effects of multiple dynamic parameters due to the instinct characteristics on the vibration behaviors of the system.
This work proposes a gear-bearing dynamic analysis model, which can consider the time-varying mesh stiffness, varied parameters such as contact pressure angle, center distance, unbalance force caused by static and dynamic eccentricity.More importantly, the fully nonlinear characteristics of bearing are built and captured by two kinds of bearing stiffnesses with definite physical meanings, which is collectively called as synthetic bearing stiffness model.The iterative dynamic process for the system dynamic model is obtained by an advanced numerical scheme, which can timely and accurately update the excitations caused by the change of instantaneous equilibrium position and the variation of displacement disturbance at each time step.The conventional constant bearing stiffness model and time-varying bearing alternating stiffness models are applied for comparison.Abundant nonlinear vibration behaviors results are obtained and discussed.

Coupled torsional-translational model for gear mesh pair
A six degree-of-freedom model is proposed for showing the translational and torsional motion for gear mesh pair in the X Y  plane.Both the pinion and gear are supported by rolling element bearings and the inner races are attached to the bearing outer raceways through interference fit.Fig. 1 shows that the pinion is applied with torque p T and the resisting torque on the gear is denoted by g T , in which the anticlockwise direction is defined as the positive in the generalized coordinate.are the original positions for the gear mesh pair, defined in the generalized coordinate; is the displacement vector introduced for describing the vibration motion of gear mesh pair in the corresponding reference coordinates.The designed center distance and contact pressure angle of gear mesh pair is denoted by 0 d and 0  , which is defined in the generalized coordinate: The center distance d and working pressure angle  are varied during dynamic process, which can be calculated through: Defining the relative mesh deformation  along the line of action is positive when the contact area is under pressure, we can obtain the calculation of  :  is defined as the position angle, obtained by:

Equations of motion
The gear mesh stiffness is defined as pg k , thus the dynamic mesh force along LOA is written as: Depending on the Newton's second law, the equations of motion for the gear-bearing system are written as: We can rewrite Eq.( 8) in matrix form as: Where: is the vibration velocity vector  is the damping matrix and the variables inside C are the viscous damping coefficients in the corresponding degree of freedom; is the bearing force vector, which is transferred from rolling element bearings and considers the asymmetric feature, varying compliance and nonlinear bearing load-deflection curves; is the mesh force vector; is the matrix of unbalance force, which is caused by the static mass eccentricity and the translational shift in X Y  plane.The detailed calculation and descriptions of excitations are presented in the following section 2.3.b F , pg F and u F are excitations for the dynamic system.For the rolling element bearing with internal clearance, the load-deflection functions are strongly sensitive to the variation of applied loads on the bearings due to the change of rolling elements in contact.In this case, the traditional dynamic analysis technique for the gear-bearing coupled system, which takes the bearing static elastic deformation characteristics by load-independent functions, is insufficient to capture the nonlinear behaviors by the complicated factors coupled with the geometric and conditional parameters.In this work, two kinds of bearing stiffnesses are introduced for comprehensively representing the transmissibility and disturbance resisting, respectively, which are called as mean stiffness and alternating stiffness.The total bearing stiffness models containing the mean stiffness and alternating stiffness is called as synthetic bearing stiffness model in the dynamic analysis part.

System excitations 2.3.1 Bearing force description
The mean stiffness denoted by b k , is based on the definition of SKF official website [24], that shows the ratio of the applied load to the elastic deformation on the corresponding direction, written as: The mean stiffness can be activated only when the displacements b u exceeds the no-load displacement 0 b u , which is examined or calculated when the external load is tending toward to nil; 0 b u at the position over one ball pass cycle showed in Fig. 3 can be written as: , where bc is the bearing internal clearance.Physically, mean stiffness of bearing presents the capacity for transferring loads on the corresponding directions.In static state, the load transferred from bearing equals to the loads applied on the bearing, thus Eq.( 10) can be rewritten as: (11) For the dynamic state, the alternating stiffness b k % of bearing is introduced for presenting the disturbance force caused by displacement perturbation around the instantaneous equilibrium position.Assuming the load-displacement curves of the bearing are smooth enough to bring in the finite difference method for calculating b k % , that is widely used as the comprehensive bearing stiffness matrix calculation [17,18,21,[25][26][27][28][29]., the calculation regime is showing in Eq.( 12) The bearing disturbance force caused by small displacement variation b u is gained by: The dynamic bearing forces b F are composed of the bearing supporting force at the dynamic equilibrium position and bearing disturbance force around the equilibrium position, written as: Thus, taking the pinion bearing as an example, the mean stiffness matrix is obtained by Eq.( 14).And the mean stiffness bp k is activated when the pinion bearing displacement bp u is larger than no-load displacement Eq.( 15) presents the alternating stiffness matrix for the pinion bearing.The non-diagonal terms bpxy k % , bpyx k % are orders of magnitude smaller than the diagonal terms bpxx k % , bpyy k % , which will be neglected in the following discussion.

Gear time-varying mesh stiffness simulation
Based on the description in in the introduction, there are plentiful publications paying close attention to the modeling of gear mesh excitation modeling and the system dynamic behaviors caused by gear mesh excitations.Considering the principle of grasping the main idea and neglecting the secondary, the meticulous geometric feature of gear mesh pair is ignored and replaced by introducing the rectangular wave mesh stiffness, which still concludes the time-varying characteristics.The gear rectangular wave mesh stiffness is expressed as: Where mpg k is the mean value of mesh stiffness; pg  is the decimal part of gear pair contact ratio; k  is the fluctuation ratio of mesh stiffness; pg T is the gear mesh cycle.

Eccentricity description
Because of static eccentricity, the mass center of pinion and gear cannot be matched with the geometric centers in real world.In addition, during the operation process the translational motion in the X Y  plane causes the dynamic eccentricity with respect to the center of rotation.For the high-speed condition, the unbalance force due to the combination of the static and dynamic eccentricity is prominent and excites the system vibration inevitably.The unbalance force matrix T 0 0 can be calculated by: In which p  , g  are the angular rotation speed of pinion and gear, respectively; p  and g  are the pinion and gear phase angle, which are presented by in Fig. 4; p e and g e are the total eccentricity of pinion and gear.The magnitude of p e and g e equal to the length of , thus the calculation process of eccentricity can be obtained geometrically through Eq.( 19): The static eccentricity is introduced and the dynamic eccentricity can be calculated based on the vibration displacements, showed as Eq.( 20): The total eccentricity is presented by Eq.( 21):

Determination of modal damping matrix
For the strong nonlinear system, the damping is extremely important for the system vibration response.It is difficult to accurately determine the damping for the practical engineering problems because that the mechanism of damping is complicated.The effect from damping on the system vibration response is beyond the scope of this work.But a method to obtain the rational value of the damping coefficients is given in the following content.The modal damping matrix will be introduced in this work, which is based on the damping ratio and corresponding eigenvalues of the system.
The bearing forces and gear mesh forces in the Eq.( 9), which are strong-nonlinear and dependent-on the variation of time and displacement, In order to get the determined bearing supporting performance and the mesh characteristics, some mathematical linearization is brought into.The equivalent bearing stiffness ˆb K and mesh stiffness ˆpg K are introduced, which are expressed in Eq.( 22) and Eq.( 23): In which , ˆˆ, , bpx bpy bgx bgy k k k k are equivalent bearing stiffnesses on the corresponding directions for a certain applied load and are only decided by the bearing alternating stiffnesses.Taking ˆbpy k as an example, the calculation process of bearing equivalent stiffness is showed as Eq.( 24 Thus, the global stiffness matrix is written as: The unbalance force vector is integrated with the external applied force vector to the equivalent external force vector, read as: Thus, the equation of motion is rewritten as: The simplified non-damping free vibration equation based on Eq.( 29) is obtained by: The eigenvalues i  ( 1, 2,..6 i  represents the degree of freedom number) are calculated by the characteristic equation, showed in Eq.( 31): (31) The mode shape matrix Thus the regular mode matrix is defined by Eq.(33): In which   , , . The Eq.(34) are obtained by the orthogonality and regularization of main modes: Where ψ is the eigenvalue matrix calculated by Eq.(35): The solution of the free equations of motion by introducing the complex vector (36) Combining with Eqs.(34)(35)(36), the Eq.( 29) can be rewritten as: The equivalent damping matrix is written by Ĉ , showing as: Ĉ is also called as regular mode damping matrix and can be calculated by: Thus the damping coefficients ( , , , , , ) can be decided by the regular modal damping matrix and the regular mode matrix.

Numerical solutions
Based on the description above, the system contains multiple strong and unpredictable nonlinear characteristics, which are intricately coupling with each other.The principal of superposition is not applicable any more, that makes it unprocurable to obtain the solution of Eq.( 8) through analytical approaches.An advanced numerical technique based on the central difference method is proposed in this work, which can comprehensively and rapidly gain the dynamic response caused by the multiple nonlinear and time-dependent excitations.The operation time t is discretized by time interval t  , which is recorded as Omitting the high order terms, Eq.( 41) is obtained by adding and subtracting Eq.( 40) Thus the Eq.( 8) at time step 1 s  is rewritten as: (42) The vectors of unbalance force, bearing supporting force and meshing force from Eq.( 9) are highly dependent on the system dynamic parameters, which can be calculated by Eqs.(43 is the bearing mean stiffness matrix from gear-bearing system, expressed as Eq.( 47): k % is the system bearing alternating stiffness matrix, is obtained by Eq(48): q is called as no-load displacement and related to the static parameters no-load displacement of pinion bearing and gear bearing 0 bp u and 0 bg u , is presented by: are the system dynamic parameters, which are decided by the magnitude of displacement vector s q .Defining the effective mass matrix and the effective force vector at time step s as: Thus the displacement vector at step 1 s  can be gotten by: Then the force vectors and the displacement disturbance . And the displacement 2 s q is obtained by repeating the similar process.The situation that the acceleration vector 0 q & & and the displacement vector which is only existing in the mathematical concepts, can be computed by Eqs.(53) (54): This numerical method is based on the central difference method, which is conditionally stable.The time interval is important for achieving the convergent solutions, which is decided by: max Where max  is the largest natural frequency of the system.k k s q 0 0 , q q s q 1 s x Fig. 5 Calculation flow chart Fig. 5 shows the calculation flow chart of this improved numerical method.We want to highlight here that other explicit iteration numerical method is also applicable as long as the excitations from the dynamic parameters inside of the equations of motion and the change of instantaneous equilibrium position and disturbance of bearing can be renewed at each time step.

Gear-bearing system parametric vibration analysis
The 6 DOFs model described above is simplified into a model with 4 DOFs for gaining the numerical results efficiently and isolating the intrinsic law between the complicated nonlinear characteristics with the nonlinear vibration behaviors of the system.For the simplified model, the pinion is assumed to have no displacements the translational direction, which is fixed to the ground.Thus the Eq.( 8) can be transformed as: Tab. 1 and Tab. 2 show the geometric parameters, mass and moment of Inertia of gear pair used in this part.The geometric parameters for the pinion supporting bearing are showed in Tab. 3 Three cases are designed to test the effect from bearing inherent nonlinearity on the system, of which the details of the nonlinear characteristics of bearings forces of different cases are described as below:  Case 1: the nonlinearity of bearing is fully captured and presented by the bearing mean stiffness and bearing alternating stiffness which is defined by Eq.( 10) and Eq.( 12), which is called as enhanced bearing force excitation model (called as EBFEM).The bearing forces are represented through Eqs.(11)(13) and computed by Eq.( 43) based on the improved numerical integration method, which are influenced not only by the time-varying compliance but also by the variation of displacements. Case2: From [25][26][27] y is the system no-load displacements, which is related to Eq.( 49) For analyzing the system vibration stability, there are two directions, one is to study the parametric resonance problem that will lead to excessive vibration strength; another one is the stability of the periodic motions with different orbitals for strong nonlinear system.For direction one-the parametric resonance problem, the vibration amplification factor y  is introduced as the system vibration strength stability evaluation reference standard, which is defined by Eq.( 57): is the peak-peak Y direction displacement when the static state at the equilibrium position over one whole cycle.Set 10 as the threshold.When the system vibration amplification factor y  is larger than 10, then the system will be judged as an unstable region with excessive vibration strength.
In this work, the load-deflection distribution is calculated from the FE/CM model [32] and it can also be obtained by other advanced analytical, computational approaches, or experiments techniques once it can fully capture and preserve the nonlinear properties of rolling element bearings.Fig. 6 shows the full distributions of bearing mean stiffness term bgy k and bearing alternating stiffness terms bgxx k % , bgyy k % .The system main internal excitations will change with the increasing of the system rotation speed, thus the system vibration responses will vary through those system internal parametric excitations.Fig. 7 is the contour map of vibration amplification factor y  comparisons from different bearing excitation models under changing rotation speed.The first order modal frequency 1 /2   calculated based on Eq.( 35) is denoted as agy f , for presenting the translational vibration modal frequency.The second order modal frequency 2 /2    is the system torsional vibration modal frequency, recorded as apg f .To facilitate analysis for the changing law of system parametric vibration, the torsional vibration modal frequency apg f is introduced as the base frequency for nondimensionalizing the frequency domain parameters, showed as:  Where m f is the system mesh cycle, gb f is the ball pass cycle for the gear bearing.

a. System vibration strength instability analysis:
The input torque p T is set as 100 N•m, and the mean value of mesh stiffness ) and 1.0, for both the three bearing excitations models the system translational and torsional vibration mode resonances are excited, the unbounded solutions that cannot converge obtained; with the change of mesh stiffness fluctuation ratio k  , the dimension of instability region is increasing, the results from and CBSM and TVBASM have larger dimension of instability region compared with the result from EBFEM.When the system locates the non-resonant area, the vibration strength from EBFEM is little larger than the results from other two models.Thus, when the excitation is within a relatively small range, the internal excitation of the system is mainly driven by gear mesh excitation.However, the nonlinear bearing support characteristics from EBFM, related to the transient dynamic equilibrium position in rolling element bearing can reduce the dimension of resonance instability region and improve the stability of the system to a certain extent.is increasing from 1.8 to 3.2, the system vibration amplification factor y  distributions in steady-state responses are showed in Fig. 7 b).With the increase of excitation frequency, the results obtained by the three models show great differences.There is a brown area for both three different bearing excitation models when ND m f  2.0, which is the super-harmonic resonant response of system torsional vibration mode, and the resonant area from EBFM is much smaller than the results from CBSM and TVBASM.Continuing to increase the system excitation frequency, there is no translational vibration mode resonant response excited by ball pass frequency For the vibration stability analysis of strong nonlinear system, in addition to the stability analysis of vibration intensity which focuses on the parametric resonance problem described as above, it also includes the stability of motion state.This work builds an enhanced bearing force excitation model for the analysis of nonlinear vibration behaviors of the gear-bearing system.The system contains the effect of multi-source excitations, which are strongly nonlinear and unpredictable.Combined with the above analysis, the vibration behavior of the system is found to have multiple periodic motions with different fundamental frequencies.Studying the vibration periodic orbits of the system under certain operation parameters and predicting the stability of the motion states with different coexisting orbits are another important direction for the system vibration instability analysis.

b. System nonlinear vibration behaviors analysis:
The system dynamic model based on EBFEM built in this work contains strong nonlinearity with coupling of multiple dynamic excitation parameters and considers the nonlinear geometric kinematic relationship of rolling element bearings and gear mesh pair, the analytical method with decided periodic solutions is not applicable.In the following analysis, the vibration distribution diagram in time and frequency domain, phase diagram and Poincare mapping diagram of the dynamic steady-state response of the system are drawn and compared for the three different bearing excitation models.According to the laws of shape and distribution, the states and regularities of the coexistence of multiple periodic motions of the system under the given operating parameters are discussed qualitatively.a) ND m f ≈0.6: Fig. 8 a) shows the different forms of presentation for the system translational vibration behaviors based on three different bearing excitation models, in which the input torque is p T =100 N•m，the mesh frequency is m f =4090 Hz ( ND m f =0.6) while the ball frequency gb f equals to 1011 Hz.Based on the results showed in Fig. 7 a), when ND m f = 0.6, the system is in the translational vibration modal resonant region.In order to avoid the excessive vibration strength, the mesh stiffness fluctuation ratio k  is chosen a relatively small value as 0.1. For CBSM, the time history diagram is a stable single cycle motion based on gear mesh cycle ; the phase diagram shows an elliptic curve and there is only one mapping point in the Poincare map; only the mesh frequency component shows in the frequency domain plot; according to these graphical features, the system vibration response based on CBSM is a harmonic response with meshing frequency as the fundamental frequency. For TVBASM, the time history result shows that the periodic motions with both meshing cycle pg T and ball pass cycle gb T can be observed obviously; the phase plane diagram shows two closed loops; Poincare mapping points are concentrated on two concentrated area; the frequency domain of FFT is mainly governed by meshing frequency m f , and there also exists some weak components of fractional harmonics of m f ( (3 / 4) m f , (5 / 4) m f ) and 2 ~ 5 harmonics of ball pass frequency gb f ; it can be seen from those results that when the bearing stiffness time-varying characteristic is considered, the system response shows multi-periodic coexisting motions based on mesh frequency m f and ball pass frequency gb f , but the stability of harmonic motion based on mesh frequency m f is the strongest. Compared with another two models, the results from EBFEM show big difference in motion morphology; the amplitude of vibration displacement peak-peak value in the time history curve is far less than the results of other two models; the periodic motion based on ball pass cycle gb T is very obvious, and the periodic motion with meshing period pg T is also visible, but the vibration amplitude changes weakly; the phase diagram and Poincare map point diagram both show the coexistence of multiple periodic orbits; the amplitude of ball pass frequency gb f component in FFT spectrum is the largest, and there are 1-6 harmonics of m f and the fractional harmonics, 2-order harmonics of meshing frequency m f ; it can be found from the results based on EBFEM, the periodic motions exists simultaneously, based on the ball pass frequency gb f and its higher-order terms and the fractional and second harmonics of meshing frequency m f , and the periodic motion stability with the fundamental frequency of ball pass frequency gb f is the strongest. For CBSM, the peak-peak value of vibration displacement in the time domain diagram decreases from 2.61 to 0.80 with the increase of system rotation speed, and the phase plane diagram is still elliptical and the area size of the ellipse in phase plane is decreasing compared to the case of 0.6. For TVBASM, the time history diagram shows complex periodic motion behaviors, and there is also a small decrease in the peak-peak value of the vibration displacement compared to the results of ND m f =0.6; the phase plane curve is non-circular and non-elliptic; Poincare mapping points are distributed in two regions on the phase plane curve; there are components of mesh frequency m f component, ball pass frequency gb f and its second harmonic in the FFT spectrum, where a continuous spectrum with a certain width is formed between the ball pass frequency gb f and mesh frequency m f ; in summary, the system vibration response based on TVBASM is in the form of periodic motions of multi periodic motions, in which the periodic motions with pg T and gb T are stable while the motion with cycle of gb T locates in the center of continuous spectrum with a tendency to instability. For EBFEM, the peak-peak value of vibration displacement is obviously larger than that of the two models; the coexistence of periodic motions with pg T and gb T can be observed; the phase plane is two groups of similar rings, and Poincare points are distributed on these two groups of rings; In FFT frequency domain diagram, the ball pass frequency gb f component is obvious, and there are weak components ofor mesh frequency fractional harmonics , .f is 4045 Hz, the system is in the translational vibration mode resonant area, which is excited by the ball pass frequency. For CBSM, the peak-peak value of vibration displacement still decreases to 0.13 μm with the increase of system rotation speed.Due to the small amplitude of vibration, the phase diagram and Poincare map are too concentrated in the interval between 35.77 and 35.90 μm.In FFT spectrum plot, there are continuous spectrum components near the fractional harmonic of mesh frequency (1 / 2) m f and the mesh frequency m f component also dominates the whole spectrum. For TVBASM, only the periodic motion with ball pass cycle gb f can be observed; therein only clear ball pass frequency gb f component in FFT spectrum.Poincare mapping points are distributed in four positions on the irregular ring, this is because the Poincare mapping surfaces are selected by the interval of mesh cycle pg T , which is a quarter of ball pass cycle gb T .In fact, when the Poincare surface is selected for the interval of ball pass cycle gb T , the Poincare mapping points of the system only exist at one concentrated point. For EBFEM, the peak-peak value of vibration displacement is obviously larger than that of the other two models, and the periodic motion with mesh cycle pg T is the most obvious in motion trajectories.The phase plane curve is an elliptic ring, and the Poincare mapping points are concentrated in a point on the ring.In the FFT spectrum, the mesh frequency m f component is very significant and the ball pass frequency gb f component is weak.
Based on three different bearing excitation models, the vibration response of the system under internal excitation is obtained.Through the highly efficient and refined numerical solution method in Chapter 3, combined with the nonlinear dynamic vibration stability analysis method, the vibration strength stability and the system motion stability of are studied by changing the internal excitations, respectively.Through three different bearing models, it can be found that:  The system based on constant bearing stiffness model (CBSM) is only excited by gear meshing, and the instability of its vibration intensity is mainly caused by the resonant response due to the change of gear parameters; when the system speed increases and passes through the resonance region, the system goes into stable state of single harmonic motion, and the vibration amplitude in the translational direction is too small, which does not meet to the reality. For the gear-bearing transmission system based on (TVBASM), the translational resonance mode caused by the bearing varying compliance excitation can be observed, but the vibration intensity in the translational resonant region is far less than that in the torsional mode resonant region; With the increase of the system rotation speed, the periodic motion of the system is mainly dominated by the stable mesh cycle motion, and the multi-periodic motions based on the ball pass frequency coexists with the vibration amplitude varying in a small range. The vibration behavior of the system based on three-dimensional nonlinear enhanced bearing force excitation model (EBFEM) is more sensitive to the change of system parameters.In the judgment of the resonance region of the system, the unstable region of vibration intensity caused by multi-source excitation can be accurately identified; In periodic motion stability analysis, when the system rotation speed is relatively low, the system periodic motion amplitude is given priority to ball pass frequency component.Along with the increase of rotation speed, system forms coexistence state of stable multiple-periodic motions, but in the resonance area, there may be the phenomenon that some periodic motion orbits are close to each other, resulting in excessive vibration intensity The fluctuated toque ap T is introduced as Eq.(59) to explore the dynamic characteristics of the system under the action of torque fluctuation excitation, gear meshing time-varying excitation and three-dimensional nonlinear excitation of rolling element bearing.Based on the above analysis, the vibration behavior is found to be more complex when the torque fluctuation excitation is considered, compared with the results when only the internal excitation is applied.Among them, the resonance region of translational vibration mode caused by bearing flexible and compliant support is still the region with the largest vibration intensity, which should be avoided firstly in the parameter design and selection of gear-bearing system.Considering the torque fluctuation excitation, the system is easier to excite the combined mode resonance region.Therefore, it is necessary to fully consider the external working conditions and judge the resonance region of the system based on the refined control model.

Conclusions
This work develops a new model for the nonlinear vibration analysis of gear-bearing system.The rolling element bearing nonlinearities are presented by mean stiffness for presenting the bearing transmissibility and alternating stiffness for showing the bearing disturbance resisting, respectively, that can fully preserve and describe the load-dependence and time-varying property of rolling bearing.For precisely and quickly solving the equations of motion of the gear-bearing system, that containing the dynamic parameters directly related to the instantaneous vibration displacements such as working pressure angle, center distance, dynamic eccentricity and bearing supporting force and disturbing force, an advanced numerical method is developed in detail.The conventional constant bearing stiffness model and time-varying bearing alternating stiffness models are applied in chapter 4 for comparison, abundant results are obtained and conclusions are summarized from the cases used in the chapter 4 are showed as: 1) The vibration behavior of the system based on three-dimensional nonlinear enhanced bearing force excitation model (EBFEM) is more sensitive to the change of system internal excitations.The unstable region of vibration intensity caused by internal excitations can be accurately identified based on EBFEM.In periodic motion stability analysis, the ball pass frequency component dominates the system periodic motion when the system rotation speed is relatively low.With the increase of rotation speed, there is coexistence state of stable multiple-periodic motions; for the resonant region, some periodic motion orbits which are close to each other, may lead to excessive vibration intensity.
2) Based on the EBFEM, the system combined mode resonant response is easier to be excited when both consider the external torque excitation and the internal excitations.

Fig. 1
Fig. 1 Model of gear-bearing rotating system

Fig. 2 Fig. 2
Fig. 2 presents the configuration and geometric relationship between the different reference coordinate.  , p p p O x y

1 YFig. 3
Fig. 3 Model of rolling element bearing ˆpg k is the mean value of mesh stiffness, is obtained through:


value of constant bearing stiffness bgy k is computed by the mean value of time-varying bearing mean stiffness bgy k over one whole ball pass cycle, showed as: is the gap function, for judging whether there is contact for the rolling elements with raceways: the peak-peak value for the system vibration displacement on the Y direction in the steady-state response region. /CM mixed method, equals to 4.40×105 N•mm -1 .Fig. 7 a) shows the results of vibration amplification factor y  when ND m f runs in range 0.4~1.8, in which the region in brown is unstable with excessive vibration strength; in this region, when the mesh frequency ND m

Fig. 7 .
Fig. 7. a) Vibration amplification factor distributions based on different bearing excitation models

Fig. 7 b
Fig. 7 b) Vibration amplification factor distributions based on different bearing excitation models

V i b r a t i o n d i s p l a c e m e n t / μ m g y 1 V i b r a t i o n v e l o c iFig. 8 af  1 2 :
Fig. 8 a) Translational vibration responses based on different bearing supporting force models

Fig. 9
Fig. 9 Translational vibration responses under both internal and external excitations .
, the bearing force excitation are only decided by the time-varying bearing alternating stiffnesses (called as TVBASM).The equations of motion based on the TVBASM, derived from Eq.(56), can obtained as: Case3: the time-varying and load-dependent characteristics are ignored for most of the researches in the past.For the comparison of enhanced bearing force excitation force model built in Chapter 3, the constant bearing stiffness model (called as CBSM) is introduced in this case.The system equations of motion based on the CBSM is got by: % 