Self-Excited Vibration Analysis of Gear-Bearing System with Multipoint Mesh and Variable Bearing Dynamic Coefficients

The gear-bearing system is the most important part of integrally centrifugal compressors. According to statistics, the majority of integrally geared compressor accidents are caused by excessive vibration of the geared rotor. However, its complicated dynamic characteristics and inevitable vibration faults in actual operation present signiﬁcant challenges throughout the analysis and design stages. In this paper, the coupled self-excited vibration of the gear system characterized by multipoint meshing and oil ﬁlm bearing supporting is investigated. Firstly, the structure of the gear system in an integrally geared compressor is used as a research object. The modeling approach of meshing excitation, including time-varying mesh stiﬀness, gear meshing error, and tooth backlash are introduced. However, the variable stiﬀness and damping coeﬃcient equations of journal bearing and oil ﬁlm thrust bearing are modeled and utilized to approximate the variable bearing force and simplify the vibration computation under the assumption of Newtonian ﬂuid. Then, a dimensionless modeling method of the gear system considering gyroscopic moment of gear disk, variable meshing force, as well as variable stiﬀness and damping coeﬃcient is proposed. Based on the dynamic model, the inﬂuence of the bearing dynamic coeﬃcients and load on the vibration of the entire gear system is studied. Among which, the vibration displacement and meshing force are examined using frequency-domain and time-domain analysis methods. The results suggest that the ﬂexible support can restrain the system’s nonlinear motion, whilst increasing load on the gear system can improve gear operation stability and reduce load ﬂuctuation.


Introduction
Integrally geared compressor is one of the most representative assembling units among big rotating machines, which are widely utilized in the domains of natural gas, petroleum, and coal chemical processing and meets the requirements, such as higher flux and longer running time. It is characterized by higher parameter, better performance and stability under various condition, extreme condition and multiinterference.
e dynamic stability of compressor geared rotor system has always attracted much attention. For example, it happened that self-excitation leads to the large vibration of compressor, as reported by Wachel (1975), Fulton (1984), Kirk(1985), Kuzdzal (1994), and Memmott (2000) [1]. e gear-bearing system is the key component in integrally centrifugal compressor which is often characterized by one helical gear meshing with many other gears at the same time and coupled in different directions, and also the nonlinear oil film support of rotors, which cause the vibration characteristics of the rotor systems in integrally geared compressors are different from those of the general gear system. According to the statistics, most accidents of integrally geared compressors are caused by excessive vibration of geared rotor [2], which is resulted by rotor unbalance, shafting coupling resonance, bearing instability, and so on. e rotor system in integrally centrifugal compressor always runs at high speed and high load. Due to the alternative engagement of teeth and the nonlinear film force, the meshing state changes during the meshing process, while the impact among the gear teeth caused by the mesh clearance and transmission error during the meshing process would excite the vibration of the whole rotor system. e excessive vibration will not only lead to the wear between rotors and stators, but also serious damage of the whole system. Although researchers have formed a relatively perfect theory and technology of nonlinear dynamic analysis and dynamic design for the main problems of the rotor system of a single shaft centrifugal compressor, and have a good means of prediction and control of its vibration problems, it is still impossible to make a reasonable explanation for the new vibration phenomenon of the geared rotor system in integrally centrifugal compressor and put forward an effective solution. In many times, the safe operation of the integrally centrifugal compressor can only be ensured by sacrificing the production capacity.
us, the meshing characteristics and performance of the gear system have important influences on the whole system even the whole devices as well as the system security.
Researchers have proposed many modeling methods on gear meshing force and oil film force. e main modeling methods of gear meshing force include linear meshing model and nonlinear meshing model. In linear meshing model, the meshing stiffness is regarded as a constant or harmonic value, which is usually used to solve the resonance frequency of the rotor or the vibration response of the rotor with less influence of meshing excitation. For example, Doan [3] presents a method for determining the resonance regions of the gear system under different design parameters according to the linear meshing assumption. Kang [4] considers the dynamics of a gear system with viscoelastic supports, in which, the gear pairs are simplified as two rigid discs connected by springs along meshing lines. In addition, in order to express the meshing excitation of gears, Shi [5] proposed a dynamics model for hypoid gears, considered the interaction between mesh stiffness and dynamic response, and the simulations show evident impact of dynamic mesh stiffness on hypoid gear dynamic response. Han [6] introduced time-dependent mesh stiffness to realize steady response analysis of rotor system, where meshing stiffness fluctuates in simple harmonic form, while the influence of various nonlinear factors on vibration is ignored, such as tooth backlash and meshing error. In many cases, the linear meshing model cannot describe the meshing force well, especially when the mesh clearance and gear transmission error to be taken into account. Feng [7] studied the increased vibration of geared rotor system caused by gear wear, and propose a gear wear model to describe wear process. With the continually updated model coefficient according to the real-time test data, the wear process of the gear mesh can be well monitored. en, Feng [8] demonstrated the ability and effectiveness of the proposed vibration-based methodology in monitoring and predicting gear wear, Inalpolat [9] analyzed the influence of different gear meshing error on dynamic meshing force of gear pair by experiment. eodossiadas [10] applied the non-linear meshing model to identify the periodic steady response of gear system involving backlash and time-dependent mesh stiffness under torsional moments. Eritenel [11] proposed nonlinear dynamic model to investigate the gear loads and bearing forces of the gear system with backlash and time-varying stiffness. Cho [12] studied the dynamics of a two-stage differential wind generator gearbox, where the non-linear meshing model was adopted. Guerine [13] investigated the influence of random uncertainty of mass, damping coefficient, bending stiffness, and torsion stiffness on the dynamic response of single-stage gear system, as well as the synthetically impact of these random parameters. Luo [14] proposed an improved analytical model to calculate the time-varying mesh stiffness of a planetary gear set, where the effect of sliding friction and spalling defects are considered. Zheng [15] considered centrifugal effect and developed an analytical-FEM framework to integrate the centrifugal field into the mesh stiffness and nonlinear dynamics, where the reasonable accuracy is demonstrated between the simulation and experiment results. Wang [16] studied the bendingtorsion coupling response of spur gear system, taking into account mesh stiffness variations, backlash, transmission errors, and loads. e results show that, due to the coupling effects of bending and torsional vibrations, the system exhibits a diverse range of periodic, subharmonic, and chaotic motion. Sánchez [17] investigated the contact conditions of modified teeth influenced by the profile modification on the load sharing and transmission error under load. Wang [18] suggested two approaches for determining time-varying meshing stiffness in internal gear pairs with tiny tooth differences, both of which were validated using the finite element method. Ma [19] established a fractal contact model suited for gear pair contact, where the influence of roughness on the normal contact stiffness of gears is considered. On the basis of modified fractal contact model, the asperity contact stiffness and the fractal contact compliance are calculated.
Also included are two methods that are commonly used in the numerical calculation of bearing force. e first one is to directly calculate the oil film force using a numerical or analytic method, while the other is an approximation of the dynamic coefficients of the oil film, such as stiffness and damping coefficients. It always takes a long time to solve the dynamic oil film force by finite element method or variation method in the first method. us, the Reynolds equation is often simplified to analyzing and solving practical problems, where the oil film force model of long bearing proposed by Sommerfeld and the short bearing proposed by Capone is the most representative. For example, Chang-Jian [20] and Amamou [21] use the long bearing hypothesis to investigate the non-linear dynamic features of a disk rotor system supported by a circular tile bearing and the dynamic stability of a circular tile bearing. However, Lin [22] and Soni [23] used the short bearing hypothesis to investigate the lubrication features of bearings in non-Newton ferromagnetic fluid and the nonlinear dynamic characteristics of circular tile bearings in thin film lubrication. Due to the complexity of the first technique's calculation process, the second method is extensively utilized in engineering applications to solve the rotor response problem due to its simple calculation procedure. Srikanth [24] ignores the degree of freedom of the tile thrust bearing block and solves the stiffness and damping parameters of tile thrust bearings of various 2 Shock and Vibration sizes. Furthermore, several researchers investigated the nonlinear vibration of a rotor system with oil film support from the perspectives of external stimulation and rotor fault. Zhang [25] investigates the change in bearing force of rotor system due to the change of submarine position. Liang [26] analyzes the influence of the film force in the squeeze film damper on the nonlinear vibration suppression of the rotor system produced by the misalignment. Cable [27] studies the static angular misalignment of the thrust collars and the main bull gear in integrally geared compressor induced by manufacturing inaccuracies and poor assembly process, which has major impact on the dynamics of the rotor system. Despite the fact that many studies have been done on the dynamic of gear system, they have always focused on gear systems with classical configurations. ere are few works concentrate on the time-varying features and non-linear coupled vibration of gear systems with multipoint mesh. Meanwhile, in some vibration analyses of complicated gear systems, bearings are always treated as invariable stiffness and damper coefficients, which would make the analysis results lose some important information. Many compressor manufacturing companies still lack a dynamic analysis process of the gear system during the design phase, resulting in many vibration problems of the gear system that cannot be effectively solved, and the safe operation of the equipment can only be ensured by sacrificing production capacity, such as reducing the working speed and load. e modeling approach of the gear system with multipoint mesh according to the structure of the gear system in integrally geared compressor is proposed in this study, where the variable meshing force and variable bearing dynamic coefficients are taken into account. In addition, the non-linear coupled vibrations of the gear system and meshing properties at various speeds and loads are investigated. e findings of this work could provide theoretical support for the dynamic design and fault diagnostics of multipoint mesh gear systems aiming for high reliability and minimal vibration.

Structure of Gear System in Integral Centrifugal
Compressor. Figure 1 depicts an entire centrifugal compressor rotor system comprised of five shafts engaged by helical gear, including one input shaft, one middle gear driving shaft, and three high speed output shafts, all supported by journal bearings and oil film trust bearings. e maximum speed of the input shaft of the rotor system is 4000r/min, and the maximum output speed of the output shaft is 12800r/min. e gears are treated as rigid body and the dynamic model of gear system is shown in Figure 2. Gears 1-5 refer to the gears on the input shaft, center shaft, and three output shafts. k ij (t), 2b ij , and e ij (t) are the time-variant stiffness parameter, gear clearance parameter, and transfer error parameter between Gears i and Gear j, respectively.
As shown in Figure 2, each gear is located in the center of the two supporting, and Gear i rotates around the rotation center O i at rotation speed ω i . k ij (t)， c ij ， e ij (t), and b ij are the meshing stiffness, meshing damping, transfer error, and unilateral backlash between Gears i and j, respectively. k ix is axial stiffness of the thrust bearing on Gear i, k iyL , k izL and k iyR , k izR are the vertical and horizontal stiffness of left and right journal bearings on Gear i, respectively.

Time-Variant Meshing Stiffness.
Due to the alternation of meshing teeth, the meshing stiffness changes periodically, which result in the dynamic stiffness excitation during the meshing process, the Ishikawa Formula shown as the following equation is used to describe the deformation of the tooth under mesh force F: where δ Br and δ Bt are the bending deformations of the rectangular and trapezoid section, respectively, and δ s and δ G are the shear and elasticity deformations, respectively. en the stiffness of the tooth can be described as e meshing stiffness appears obvious periodicity during the meshing process and can be expressed as a Fourier series: where k 0 is the average meshing stiffness; a i , b i , and k i are the Fourier coefficients; φ j is the phase angle; and ω Λ is the meshing frequency. If only the first order of Fourier series is preserved, Eq. (3) can be written as where ε is the fluctuation coefficient of stiffness.

Gear Transfer Error.
e transfer error of the gear pair is mainly caused by manufacture error and gear wear. Transfer errors make the meshing tooth profile deviate from the ideal meshing position, which destroy the correct meshing method of gear pairs, and resulting in tooth-totooth collisions and impacts during the meshing process.
Because the transfer error is of periodicity during the transfer process, the transfer error could also be expressed in the form of Fourier series: Shock and Vibration 3  Similarly, if only the first order of this Fourier series is preserved, and average transfer error is equaled to zero, Equation (3) can be written as where e s is the fluctuant amplitude of transfer error.

Backlash.
Backlash between gear pair changes the contact state and causes continuous impact between the gear pairs, which has a significant impact on the dynamic characteristics of the gear system. Assuming the teeth backlash of a gear pair is 2b n , and the relative deformation between the meshing teeth is Δ. erefore, the deformation functions of the tooth along the meshing line can be expressed as Equation (7) is a step function, and when the meshing point is within the scope of teeth backlash, there is no deformation. Otherwise, the teeth begin to deform.

Modelling Approach of Bearing Dynamic Coefficients.
e journal bearings on shafts can be modeled as Figure 3. Here, the journal rotates at the speed of Ω with stable load F. O and Ο b are the center of journal and bearing, respectively. r and R b are the radius of journal and bearing bore. φ is the displacement angle. E is the eccentric distance.
After coordinate translation, the Reynolds equation for oil film force analysis in journal bearing and its oil film boundary can be shown as where η is the dynamic viscosity of the oil; p is the oil film pressure; x is the position along the width direction; θ 1 and θ 2 are the starting angle boundary and end angle boundary; and L is the width of the bearing. Nondimensional parameters are defined as To simplify the vibration calculations, stiffness and damping coefficients can be used to approximate the bearing force of the lubricant film, which can be calculated by linearization of the unsteady load capacity in the vicinity of the static position of equilibrium and developed to the first derivative in a Taylor series as Equation (11), where the force caused by dip angle of the journal is neglected: where k yy k yz k zy k zz c yy c yz c zy c zz e oil film thrust bearing on shafts can be modeled as Figure 4, which consists of many bearing bushes. θ 0 , θ 1 , and θ 2 are the angle of bush, side oil seal lip, and thrust face, respectively; b is the width of the outer oil seal lip; r i and r o are the inner and outer radius of thrust face; φ is the angle of thrust face; θ and r are the coordinates in radial and circumferential directions; h(θ) is the film thickness with radial coordinate θ; C is the bearing clearance; and e is the relative displacement between trust bearing and thrust runner collar. e Reynolds equation for each bush in oil film thrust bearing and its oil film boundary can be shown as where η is the dynamic viscosity of the oil; p is the oil film pressure; and θ 1 and θ 2 are the starting angle boundary and end angle boundary. Nondimensional parameters are defined as e/(CΩ), τ � Ωt, and κ � 1/r. Equations (14)- (15) can be transferred to Shock and Vibration Similarly, the linearization result of nonlinear bearing force can be shown as Equation (18), where the force caused by dip angle of the thrust runner collar is neglected: e turnover stiffness of geared shaft to gear is expressed as Define k iyL and k izL as the left supporting bearing stiffness of Gear i along axis-x, y, and define k iyR and k izR as the right supporting bearing stiffness of Gear i along axis-x, y. e lateral supporting stiffness of Gear is expressed as e tilting stiffness of Gear i is expressed as Define c ix , c iy , c iz , c iθy , and c iθz as the shaft damping coefficient of Gear i along axis-x, y, z and around axis-y, z, respectively. e dynamic functions of gear system shown in Figure 2 are shown in Equations (22)-(26). e dynamic functions of Gear 1 are shown as follows: e dynamic functions of Gear 2 are shown as follows: Shock and Vibration e dynamic functions of Gear 3 are shown as follows: e dynamic functions of Gear 4 are shown as follows: e dynamic functions of Gear 5 are shown as follows: 8 Shock and Vibration where T i is the load on the Gear i; m i , I pi , and I di are the mass, polar moment of inertia, and diameter moment of inertia of Gear i, respectively; p ij is the length of meshing line between Gear i and Gear j; ϕ ij is the position angle of Gear j respect to Gear i; α is the pressure angle of gear; and g ij (t) is the clearance function expressed as where 2b ij is the clearance length between Gear i and j. e space vector of the system totally has 30 DOFs and can be expressed as T . (28) Equations (22)- (26) can be written as where M, C s , K s , C m , and K m are the mass matrix, supporting damping matrix, supporting stiffness matrix, meshing damping matrix, and meshing stiffness matrix of gear system, respectively; H is the gyro matrix of the gear system; G(p) is a tooth gap function related to p; and F is load vector. rough defining that x θi � r i θ xi , y θi � r i θ yi , and z θi � r i θ zi , the angle variable of space vector can be transformed into the arc length turning around the basic circle. Assume the characteristic frequency and length are ω c and b c , and define that S � diag 1, 1, 1, r 1 , r 1 , r 1 , · · · , 1, 1, 1, r 5 , r 5 , r 5 .
(30) (25) can be converted into proper dimensionless indexes and expressed as c F。 e backlash function of each gear pair can be defined as

Vibration Behaviors of the Gear System by Self-Excited
Vibration. e main parameters of gear pairs in Figure 2 are shown in Table 1. e backlash and meshing error of the gear pairs are 100 μm and 20 μm, respectively.

Condition 2. Rated load and variable bearing damping and stiffness calculated by Equations
e Newmark-β method is used to analyze the dynamics of the above three working conditions according to (31), in which the Newmark constants are 0.25 and 0.5. In order to describe the vibration and meshing behavior of the gear system in the whole process of increasing speed, the vibration bifurcation diagrams of Gear 4 along y-axis at different input speeds ranging from 500r/min to 4000r/min are obtained and presented in Figure 5 based on the dynamic analysis of the vibration of gear system under the three conditions. e nonlinear vibration characteristics of the gear system under the impact of speed are shown in Figure 5(a). Before the input speed reaches 1600r/min, the gear vibration Shock and Vibration 9 is harmonic.
en the vibration signals begin to period doubling bifurcate, and gradually evolve into chaos near 2280 r/min. With the increasing of the input speed, the vibration signal of the gear begins to return to the multi period motion, and eventually returning to harmonic motion at 3500r/min.
When bearing variable stiffness is taken into account, nonlinear vibration of gears is reduced. Although there is no further evolution after the vibration signal evolved to eight periods' motion, the non-linear speed region, which is also from 1600r/min to 3500r/min, is not greatly reduced in Figure 5(b). e stiffness and the damping of the bearing are positively correlated with the speed of the bearing, according to Equations (13) and (17). e dynamic stiffness of each oil film bearing before 4000 r/min is relatively low in comparison to Condition 1, so the stiffness of the entire gear system with Condition 2 is lower than Condition 1. ese show that the flexible support has inhibition to the nonlinear    motion of the system, whereas the support with greater stiffness could aggravate the nonlinear vibration of the system. e calculated result of the gear system under 1.5 times rated load differs significantly from that shown in Figure 5(c) under rated load. With a heavy load, the non-linear speed interval decreases noticeably. e vibration of the system begins to bifurcate after 1800r/min, and there is no further evolution after two periods' motion and back to harmonic motion after 3000r/min. e gear vibration signal alternates between two periods' motion and on period motion throughout the entire input speed range. e comparison results show that the load on the gear system has a significant impact on the gear system's nonlinear vibration. Figure 6 depicts the frequency-domain characteristics of the vibration displacement of Gear 4 along y-direction under the three conditions. ere is no discernible difference between Condition 1 and Condition 2 in terms of frequencydomain features. Except for the meshing frequency f m and its double frequency of 2f m , the subfrequency f m /2 begins to appear from 1500r/min to 3500r/min. F m /4 and f m /8 can also be found from 2000 r/min to 3000 r/min, and the peaks of the subfrequencies are even higher than the peak of meshing frequency at some speeds. e peak value of meshing frequency gradually decreases as the speed increases. It also shows that if nonlinear vibration in the nonlinear speed region is ignored, the higher the speed is, the relatively less the effect of self-excited vibration on the system vibration is. e frequency-domain features of Condition 3 are quite different from those of the first two conditions. Except for the meshing frequency f m and its double frequency of 2f m , only the subfrequency f m /2 can be found between 1700r/min to 3000r/min, and no other subfrequencies are presented. e peak value of subfrequency f m /2 with heavy load is much lower than in Conditions 1 and 2, due to the load effect, but there is no discernible difference in the peak value of f m compared to Conditions 1 and 2. e influence of variable bearing coefficients on the analysis results is compared using a correlation analysis between Condition 1 and Condition 2. Figure 7 depicts the results of the vibration signals in the time domain and frequency domain at various speeds. When the input speed is below 1000 r/min or above 3200 r/min, the correlation coefficients between time-domain and frequency-domain data are close to 1. e correlation coefficient of frequencydomain data decreases slightly in the nonlinear region, but it remains above 0.95. However, time-domain data has lower correlation coefficients than frequency-domain data, though there is still some correlation (above 0.6). It shows that variable stiffness and damper of bearing have a minor impact on the simulation results at frequency-domain features analysis, but have a significant impact at time-domain features analysis. Variable bearing coefficients should be considered in vibration analysis and non-linear characteristic analysis. e vibration signals of y4 under Condition 1 at different rotational speeds are shown in Figure 8 to compare and analyze the change of gear vibration with input speed. e variation of the vibration period with rotational speed can be clearly seen from the diagram. In order to compare the vibration information under the three conditions more clearly. Time-domain analysis methods are adopted in Figure 9, where the peak to peak value and standard deviation of vibration data in these three conditions are analyzed.
It can be seen that the peak to peak and standard deviation values of the vibration data in the non-linear speed region are much larger than the other speeds. e vibration will act directly on the supporting bearing, producing disturbed vibration force. e peak to peak value and standard deviation of the central displacement of Gear 4 are nearly identical under Conditions 1 and 2, but the result calculated using the variable dynamic parameter of bearing is slightly lower and smoother than constant one. Due to the influence of heavy load, the peak to peak value and standard deviation values are clearly reduced, while the fluctuated force of the bearing in the non-linear speed region is also relatively small.

Meshing Force of the Gear System by Self-Excited
Vibration. Figure 10 depicts the frequency-domain features of meshing force between Gear 2 and Gear 4 at various speeds.
ere is also no discernible difference between the frequency-domain features of meshing force under Conditions 1 and 2. In comparison to the results in Figure 6, the peaks of meshing frequency are unaffected by changing speed, and they are significantly greater than the peaks of subfrequency f m /2 even in the non-linear speed region. f m /4 and f m /8 can also be found from 1500 r/min to 3500 r/min, but their peaks of these subfrequencies are far smaller than f m /2. e frequency-domain features of meshing force differ significantly from those shown in Figure 6(c). Except for the meshing frequency f m and its double frequencies, the subfrequencies are quite small. Condition 3 has significantly higher peaks of meshing frequency than the other conditions due to the heavy load. Figure 11 compares and analyzes the mean, peak to peak, and RMS of the meshing force in these three conditions. e mean value of the three groups of meshing forces varies little as the rotational speed increases. e values under Condition 3 are roughly 1.5 times those of Conditions 1 and 2, which are proportional to the load. e peak to peak values and RMS values in the nonlinear speed region are obviously higher than the other speeds at rated load. e result calculated using a variable dynamic parameter of bearing changes slightly smoother than constant one. When the load have increased to 1.5 times rated load, the peak to peak values and RMS values in the nonlinear speed region do not increase, and are far lower than the results under rated loads. Meanwhile, the RMS values under high load are clearly higher than that under low load. e main cause of non-linear vibration in compressor gear systems is the impact between gears. With the increasing of the speed, the meshing state between the tooth surfaces would change from continuous meshing to collision, and different non-linear meshing states could be evolved. In the nonlinear speed region, the vibration and    14 Shock and Vibration meshing force will change. In general, increasing the load and decreasing the stiffness of the bearing appropriately can improve the stability of meshing.

Conclusion
e modeling approach of the geared-rotor system in integrally centrifugal compressors considering changing multipoint meshing force and variable bearing stiffness and damper is proposed in this research, and the non-linear coupled vibrations of the gear system and meshing properties are explored at different speeds and loads. e nonlinear vibration might appear as the meshing frequency increases due to the effect of nonlinear parameters such as meshing clearance and transmission error. e flexible support inhibits the system's nonlinear motion, whereas the bearing with increased stiffness may aggravate the system's nonlinear vibration. Considering the influence of the rotating shaft's dynamic stiffening effect, increasing the driven speed will not change the peak of meshing frequency of meshing force, but decrease the peak of meshing frequency of vibration in the gear system. e peak to peak and standard deviation values of the vibration data are substantially bigger in the nonlinear speed area than in the other speeds. is means that the disturbance force acting on the supporting bearing would be raised as well. e disrupted force in the nonlinear speed area would plainly diminish as the load increased. e peak to peak and RMS values of the meshing force are greater in the nonlinear speed area than in the other speeds at rated load. However, when the load is increased to a certain level, this phenomenon disappears, and the changes in peak to peak and RMS values in the nonlinear speed area are not significant.
e large load of the gear-bearing system can effectively reduce nonlinear vibration and increase meshing stability.

Data Availability
e data used to support the findings of this study are included within the article.

Conflicts of Interest
e authors declare that there are no conflicts of interest regarding the publication of this paper.