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

Authors: Hao Zhang; Shiheng Cao; Pengyu Li; Qingkai Han
Journal: Shock and Vibration
 (2022)

---

## Abstract

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 significant 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 film 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 stiffness, gear meshing error, and tooth backlash are introduced. However, the variable stiffness and damping coefficient equations of journal bearing and oil film thrust bearing are modeled and utilized to approximate the variable bearing force and simplify the vibration computation under the assumption of Newtonian fluid. Then, a dimensionless modeling method of the gear system considering gyroscopic moment of gear disk, variable meshing force, as well as variable stiffness and damping coefficient is proposed. Based on the dynamic model, the influence of the bearing dynamic coefficients 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 flexible support can restrain the system’s nonlinear motion, whilst increasing load on the gear system can improve gear operation stability and reduce load fluctuation.


---

## Body

## 1. 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 multi-interference. The 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].The 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. The 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. The 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. Thus, 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. The 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. Then, 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. Theodossiadas [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 bending-torsion coupling response of spur gear system, taking into account mesh stiffness variations, backlash, transmission errors, and loads. The 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. The 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. Thus, 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 sizes. Furthermore, several researchers investigated the non-linear 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. There 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.The 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. The 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.


## 2. Gear Bearing System with Multipoint Mesh and Its Modelling Method

### 2.1. Structure of Gear System in Integral Centrifugal Compressor

Figure1 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. The 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.Figure 1
Rotor system of integral centrifugal compressor.The gears are treated as rigid body and the dynamic model of gear system is shown in Figure2. Gears 1–5 refer to the gears on the input shaft, center shaft, and three output shafts. kijt, 2bij, and eijt are the time-variant stiffness parameter, gear clearance parameter, and transfer error parameter between Gears i and Gear j, respectively.Figure 2
Dynamic model of the gear system.As shown in Figure2, each gear is located in the center of the two supporting, and Gear i rotates around the rotation center Oi at rotation speed ωi. kijt， cij， eijt, and bij are the meshing stiffness, meshing damping, transfer error, and unilateral backlash between Gears i and j, respectively. kix is axial stiffness of the thrust bearing on Gear i, kiyL, kizL and kiyR, kizR are the vertical and horizontal stiffness of left and right journal bearings on Gear i, respectively.


### 2.2. Modelling Approach of Meshing Excitation

#### 2.2.1. 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 forceF:(1)δ=δBr+δBt+δs+δG,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. Then the stiffness of the tooth can be described as(2)k=Fδ.The meshing stiffness appears obvious periodicity during the meshing process and can be expressed as a Fourier series:(3)kt=k0+∑j=1∞ajcosjωΛt+bjsinjωΛt=k0+∑j=1∞kjcosjωΛt+φj,where k0 is the average meshing stiffness; ai, bi, and ki 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(4)kt=k01+εsinωΛt+φ,where ε is the fluctuation coefficient of stiffness.


#### 2.2.2. Gear Transfer Error

The 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-to-tooth 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:(5)et=e0+∑j=1∞cjcosjωt+djsinjωt=e0+∑j=1∞ejcosjωΛt+φj,where e0 is the average transfer error and ci and di are the Fourier coefficients.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(6)et=essinωΛt+φ,where es is the fluctuant amplitude of transfer error.


#### 2.2.3. 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 is2bn, and the relative deformation between the meshing teeth is Δ. Therefore, the deformation functions of the tooth along the meshing line can be expressed as(7)gΔ=Δ−bn,Δ>bn,0,−bn≤Δ≤bn,Δ+bn,Δ<−bn,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.


### 2.3. Modelling Approach of Bearing Dynamic Coefficients

The journal bearings on shafts can be modeled as Figure3. 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 Rb are the radius of journal and bearing bore. φ is the displacement angle. E is the eccentric distance.Figure 3
Dynamic model of the journal bearing.After coordinate translation, the Reynolds equation for oil film force analysis in journal bearing and its oil film boundary can be shown as(8)1Rb2∂∂θh312η∂p∂θ+∂∂xh312η∂p∂x=12Ω∂h∂θ+∂h∂t,(9)pθ1,x=pθ2,x=0,pθ,−L2=pθ,L2=0,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 asp¯=p/6ηΩRb/C2, x¯=x/L/2， λ=2Rb/L, e¯=e/C, e¯˙=e˙/CΩ, h¯=h/C=1+e¯cosθ−φ, and τ=Ωt. Equation (8)-(9) can be transferred to(10)∂∂θh¯3∂p¯∂θ+λ2∂∂x¯h¯3∂p¯∂x¯=−e¯sinθ−φ+2e¯˙cosθ−φde=fθ,pθ1,x¯=pθ2,x¯=0,pθ,−1=pθ,1=0.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:(11)FyFz=Fy0Fz0+kyykyzkzykzzyz+cyycyzczyczzy˙z˙,where(12)kyykyzkzykzz=∫−L/2L/2∫θ1θ2∂p∂ycosθ−∂p∂ysinθ∂p∂zcosθ−∂p∂zsinθRbdθdx,cyycyzczyczz=∫−L/2L/2∫θ1θ2∂p∂y˙cosθ−∂p∂y˙sinθ∂p∂z˙cosθ−∂p∂z˙sinθRbdθdx.The oil film thrust bearing on shafts can be modeled as Figure4, 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; ri and ro 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.Figure 4
Dynamic model of the oil film thrust bearing.The Reynolds equation for each bush in oil film thrust bearing and its oil film boundary can be shown as(13)1r2∂∂θh312η∂p∂θ+∂∂rh312η∂p∂r=12Ω∂h∂θ+∂h∂t,(14)pri,θ=pro,θ=0,pr,θ1=pr,θ2=0,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 asp¯=p/6ηΩri/C2, r¯=r/ri, b¯=b/ri, h¯=h/C, e¯=e/C, e¯˙=e˙/CΩ, τ=Ωt, and κ=1/r¯. Equations (14)-(15) can be transferred to(15)κ2∂∂θh¯3∂p¯∂θ+∂∂r¯h¯3∂p¯∂r¯=∂h¯∂θ+2∂h¯∂τ,p1,θ=pr¯o,θ=0,pr¯,θ1=pr¯,θ2=0.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:(16)Fx≈Fx0+kxxΔx+cxxΔx˙,where(17)kxx=−∫riro∫θ1θ2∂p∂xrdθdr,cxx=−∫riro∫θ1θ2∂p∂x˙rdθdr.


### 2.4. Dynamic Model of the Gear-Bearing System

For the shaft of Geari, if total length is li, average diameter is di, and Young’s modulus is E, the supporting stiffness of the shaft to Gear i is expressed as(18)ks=3πdi4E4li3.The turnover stiffness of geared shaft to gear is expressed as(19)ko=πdi4E2li.DefinekiyL and kizL as the left supporting bearing stiffness of Gear i along axis-x, y, and define kiyR and kizR as the right supporting bearing stiffness of Gear i along axis-x, y. The lateral supporting stiffness of Gear is expressed as(20)kiy=kiyL+kiyRkskiyL+kiyR+ks,kiz=kizL+kizRkskizL+kizR+ks.The tilting stiffness of Geari is expressed as(21)kiθy=likizL+kizRkolikizL+kizR+2ko,kiθz=likiyL+kiyRkolikiyL+kiyR+2ko.Definecix, ciy, ciz, ciθy, and ciθz as the shaft damping coefficient of Gear i along axis-x, y, z and around axis-y, z, respectively. The dynamic functions of gear system shown in Figure 2 are shown in Equations (22)–(26).The dynamic functions of Gear 1 are shown as follows:(22)m1x¨1+c1xx˙1+k1xx1−c12p˙12+k12g12tsinβ=0,m1y1+c1yy˙1+k1yy1+c12p˙12+k12g12tcosβsinα−ϕ12=0,m1z1+c1zz˙1+k1zz1+c12p˙12+k12g12tcosβcosα−ϕ12=0,Ip1θ¨x1+c1θxθ˙x1+r1c12p˙12+k12g12tcosβ=T1,Id1θ¨y1+c1θyθ˙y1+ω1Ip1θ˙z1+k1θyθy1+r1c12p˙12+k12g12tsinβsinα−ϕ12=0,Id1θ¨z1+c1θzθ˙z1−ω1Ip1θ˙y1+k1θzθz1+r1c12p˙12+k12g12tsinβcosα−ϕ12=0.The dynamic functions of Gear 2 are shown as follows:(23)m2x¨2+c2xx˙2+k2xx2+c12p˙12+k12g12tsinβ−∑i=35c2ip˙2i+k2ig2itsinβ=0,m2y¨2+c2yy˙2+k2yy2−c12p˙12+k12g12tcosβsinα−ϕ12,+∑i=35c2ip˙2i+k2ig2itcosβsinα−ϕ12=0,m2z¨2+c2zz˙2+k2zz2−c12p˙12+k12g12tcosβcosα−ϕ12,+∑i=35c2ip˙2i+k2ig2itcosβcosα−ϕ12=0,Ip2θ¨2+c2θxθ˙x2−r2c12p˙12+k12g12tcosβ+r2∑i=35c2ip˙2i+k2ig2itcosβ=0,Id2θ¨y2+c2θyθ˙y2+ω2Ip2θ˙z2+k2θyθy2+r2c12p˙12+k12g12tsinβsinα−ϕ12,+r2∑i=35c2ip˙2i+k2ig2itsinβsinα−ϕ12=0,Id2θ¨z2+c2θzθ˙z2−ω2Ip2θ˙y2+k2θzθz2+r2c12p˙12+k12g12tsinβcosα−ϕ12,+r2∑i=35c2ip˙2i+k2ig2itsinβcosα−ϕ12=0.The dynamic functions of Gear 3 are shown as follows:(24)m3x¨3+c3xx˙3+k3xx3+c23p˙23+k23g23tsinβ=0,m3y¨3+c3yy˙3+k3yy3−c23p˙23+k23g23tcosβsinα−ϕ23=0,m3z¨3+c3zz˙3+k3zz3−c23p˙23+k23g23tcosβcosα−ϕ23=0,Ip3θ¨x3+c3θxθ˙x3−r3c23p˙23+k23g23tcosβ=T3,Id3θ¨y3+c3θyθ˙y3+ω3Ip3θ˙z3+k3θyθy3+r3c23p˙23+k23g23tsinβsinα−ϕ23=0,Id3θ¨z3+c3θzθ˙z3−ω3Ip3θ˙y3+k3θzθz3+r3c23p˙23+k23g23tsinβcosα−ϕ23=0.The dynamic functions of Gear 4 are shown as follows:(25)m4x¨4+c4xx˙4+k4xx4+c24p˙24+k24g24tsinβ=0,m4y¨4+c4yy˙4+k4yy4−c24p˙24+k24g24tcosβsinα−ϕ24=0,m4z¨4+c4zz˙4+k4zz4−c24p˙24+k24g24tcosβcosα−ϕ24=0,Ip4θ¨x4+c4θxθ˙x4−r4c24p˙24+k24g24tcosβ=T4,Id4θ¨y4+c4θyθ˙y4+ω4Ip4θ˙z4+k4θyθy4+r4c24p˙24+k24g24tsinβsinα−ϕ24=0,Id4θ¨z4+c4θzθ˙z4−ω4Ip4θ˙y4+k4θzθz4+r4c24p˙24+k24g24tsinβcosα−ϕ24=0.The dynamic functions of Gear 5 are shown as follows:(26)m5x¨5+c5xx˙5+k5xx5+c25p˙25+k25g25tsinβ=0,m5y¨5+c5yy˙5+k5yy5−c25p˙25+k25g25tcosβsinα−ϕ25=0,m5z¨5+c5zz˙5+k5zz5−c25p˙25+k25g25tcosβcosα−ϕ25=0,Ip5θ¨x5+c5θxθ˙x5−r5c25p˙25+k25g25tcosβ=T,5Id5θ¨y5+c5θyθ˙y5+ω5Ip5θ˙z5+k5θyθy5+r5c25p˙25+k25g25tsinβsinα−ϕ25=0,Id5θ¨z5+c5θzθ˙z5−ω5Ip5θ˙y5+k5θzθz5+r5c25p˙25+k25g25tsinβcosα−ϕ25=0,where Ti is the load on the Gear i; mi, Ipi, and Idi are the mass, polar moment of inertia, and diameter moment of inertia of Gear i, respectively; pij 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 gijt is the clearance function expressed as(27)gijt=pij−bij,pij>bij,0,pij≤bij,pij+bij,pij<−bij,where 2bij is the clearance length between Gear i and j.The space vector of the system totally has 30 DOFs and can be expressed as(28)X=x1,y1,z1,θx1,θy1,θz1,⋯,x5,y5,z5,θx5,θy5,θz5T.Equations (22)-(26) can be written as(29)MX¨+Cs+Cm+HX˙+KsX+KmGp=F,where M, Cs, Ks, Cm, and Km 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; Gp is a tooth gap function related to p; and F is load vector.Through defining thatxθi=riθxi, yθi=riθyi, and zθi=riθ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 bc, and define that(30)τ=ωct,X¯ij=Xijbc,e¯ij=eijbc,b¯ij=bijbc,X¯=x¯1,y¯1,z¯1,x¯θ1,y¯θ1,z¯θ1,⋯,x¯5,y¯5,z¯5,x¯θ5,y¯θ5,z¯θ5T,dX¯dτ=X¯˙,d2X¯dτ2=X¯¨,S=diag1,1,1,r1,r1,r1,⋯,1,1,1,r5,r5,r5.(25) can be converted into proper dimensionless indexes and expressed as(31)X¯¨+C¯s+C¯m+H¯X¯˙+K¯sX¯+K¯mG¯p¯=F¯,where C¯s=SM−1CsS−1/ωc， C¯m=SM−1CmS−1/ωc， K¯s=SM−1KsS−1/ωc， K¯m=SM−1KmS−1/ωc2， H¯s=SM−1HS−1/ωc, and F¯=SM−1/bcωc2F。The backlash function of each gear pair can be defined as(32)G¯ijp¯ij=X¯ij−14SijBijb¯ij,p¯ij>b¯ij,0,p¯ij≤b¯ij,X¯ij+14SijBijb¯ij,p¯ij<−b¯ij.


## 2.1. Structure of Gear System in Integral Centrifugal Compressor

Figure1 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. The 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.Figure 1
Rotor system of integral centrifugal compressor.The gears are treated as rigid body and the dynamic model of gear system is shown in Figure2. Gears 1–5 refer to the gears on the input shaft, center shaft, and three output shafts. kijt, 2bij, and eijt are the time-variant stiffness parameter, gear clearance parameter, and transfer error parameter between Gears i and Gear j, respectively.Figure 2
Dynamic model of the gear system.As shown in Figure2, each gear is located in the center of the two supporting, and Gear i rotates around the rotation center Oi at rotation speed ωi. kijt， cij， eijt, and bij are the meshing stiffness, meshing damping, transfer error, and unilateral backlash between Gears i and j, respectively. kix is axial stiffness of the thrust bearing on Gear i, kiyL, kizL and kiyR, kizR are the vertical and horizontal stiffness of left and right journal bearings on Gear i, respectively.


## 2.2. Modelling Approach of Meshing Excitation

### 2.2.1. 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 forceF:(1)δ=δBr+δBt+δs+δG,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. Then the stiffness of the tooth can be described as(2)k=Fδ.The meshing stiffness appears obvious periodicity during the meshing process and can be expressed as a Fourier series:(3)kt=k0+∑j=1∞ajcosjωΛt+bjsinjωΛt=k0+∑j=1∞kjcosjωΛt+φj,where k0 is the average meshing stiffness; ai, bi, and ki 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(4)kt=k01+εsinωΛt+φ,where ε is the fluctuation coefficient of stiffness.


### 2.2.2. Gear Transfer Error

The 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-to-tooth 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:(5)et=e0+∑j=1∞cjcosjωt+djsinjωt=e0+∑j=1∞ejcosjωΛt+φj,where e0 is the average transfer error and ci and di are the Fourier coefficients.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(6)et=essinωΛt+φ,where es is the fluctuant amplitude of transfer error.


### 2.2.3. 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 is2bn, and the relative deformation between the meshing teeth is Δ. Therefore, the deformation functions of the tooth along the meshing line can be expressed as(7)gΔ=Δ−bn,Δ>bn,0,−bn≤Δ≤bn,Δ+bn,Δ<−bn,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.


## 2.2.1. 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 forceF:(1)δ=δBr+δBt+δs+δG,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. Then the stiffness of the tooth can be described as(2)k=Fδ.The meshing stiffness appears obvious periodicity during the meshing process and can be expressed as a Fourier series:(3)kt=k0+∑j=1∞ajcosjωΛt+bjsinjωΛt=k0+∑j=1∞kjcosjωΛt+φj,where k0 is the average meshing stiffness; ai, bi, and ki 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(4)kt=k01+εsinωΛt+φ,where ε is the fluctuation coefficient of stiffness.


## 2.2.2. Gear Transfer Error

The 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-to-tooth 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:(5)et=e0+∑j=1∞cjcosjωt+djsinjωt=e0+∑j=1∞ejcosjωΛt+φj,where e0 is the average transfer error and ci and di are the Fourier coefficients.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(6)et=essinωΛt+φ,where es is the fluctuant amplitude of transfer error.


## 2.2.3. 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 is2bn, and the relative deformation between the meshing teeth is Δ. Therefore, the deformation functions of the tooth along the meshing line can be expressed as(7)gΔ=Δ−bn,Δ>bn,0,−bn≤Δ≤bn,Δ+bn,Δ<−bn,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.


## 2.3. Modelling Approach of Bearing Dynamic Coefficients

The journal bearings on shafts can be modeled as Figure3. 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 Rb are the radius of journal and bearing bore. φ is the displacement angle. E is the eccentric distance.Figure 3
Dynamic model of the journal bearing.After coordinate translation, the Reynolds equation for oil film force analysis in journal bearing and its oil film boundary can be shown as(8)1Rb2∂∂θh312η∂p∂θ+∂∂xh312η∂p∂x=12Ω∂h∂θ+∂h∂t,(9)pθ1,x=pθ2,x=0,pθ,−L2=pθ,L2=0,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 asp¯=p/6ηΩRb/C2, x¯=x/L/2， λ=2Rb/L, e¯=e/C, e¯˙=e˙/CΩ, h¯=h/C=1+e¯cosθ−φ, and τ=Ωt. Equation (8)-(9) can be transferred to(10)∂∂θh¯3∂p¯∂θ+λ2∂∂x¯h¯3∂p¯∂x¯=−e¯sinθ−φ+2e¯˙cosθ−φde=fθ,pθ1,x¯=pθ2,x¯=0,pθ,−1=pθ,1=0.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:(11)FyFz=Fy0Fz0+kyykyzkzykzzyz+cyycyzczyczzy˙z˙,where(12)kyykyzkzykzz=∫−L/2L/2∫θ1θ2∂p∂ycosθ−∂p∂ysinθ∂p∂zcosθ−∂p∂zsinθRbdθdx,cyycyzczyczz=∫−L/2L/2∫θ1θ2∂p∂y˙cosθ−∂p∂y˙sinθ∂p∂z˙cosθ−∂p∂z˙sinθRbdθdx.The oil film thrust bearing on shafts can be modeled as Figure4, 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; ri and ro 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.Figure 4
Dynamic model of the oil film thrust bearing.The Reynolds equation for each bush in oil film thrust bearing and its oil film boundary can be shown as(13)1r2∂∂θh312η∂p∂θ+∂∂rh312η∂p∂r=12Ω∂h∂θ+∂h∂t,(14)pri,θ=pro,θ=0,pr,θ1=pr,θ2=0,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 asp¯=p/6ηΩri/C2, r¯=r/ri, b¯=b/ri, h¯=h/C, e¯=e/C, e¯˙=e˙/CΩ, τ=Ωt, and κ=1/r¯. Equations (14)-(15) can be transferred to(15)κ2∂∂θh¯3∂p¯∂θ+∂∂r¯h¯3∂p¯∂r¯=∂h¯∂θ+2∂h¯∂τ,p1,θ=pr¯o,θ=0,pr¯,θ1=pr¯,θ2=0.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:(16)Fx≈Fx0+kxxΔx+cxxΔx˙,where(17)kxx=−∫riro∫θ1θ2∂p∂xrdθdr,cxx=−∫riro∫θ1θ2∂p∂x˙rdθdr.


## 2.4. Dynamic Model of the Gear-Bearing System

For the shaft of Geari, if total length is li, average diameter is di, and Young’s modulus is E, the supporting stiffness of the shaft to Gear i is expressed as(18)ks=3πdi4E4li3.The turnover stiffness of geared shaft to gear is expressed as(19)ko=πdi4E2li.DefinekiyL and kizL as the left supporting bearing stiffness of Gear i along axis-x, y, and define kiyR and kizR as the right supporting bearing stiffness of Gear i along axis-x, y. The lateral supporting stiffness of Gear is expressed as(20)kiy=kiyL+kiyRkskiyL+kiyR+ks,kiz=kizL+kizRkskizL+kizR+ks.The tilting stiffness of Geari is expressed as(21)kiθy=likizL+kizRkolikizL+kizR+2ko,kiθz=likiyL+kiyRkolikiyL+kiyR+2ko.Definecix, ciy, ciz, ciθy, and ciθz as the shaft damping coefficient of Gear i along axis-x, y, z and around axis-y, z, respectively. The dynamic functions of gear system shown in Figure 2 are shown in Equations (22)–(26).The dynamic functions of Gear 1 are shown as follows:(22)m1x¨1+c1xx˙1+k1xx1−c12p˙12+k12g12tsinβ=0,m1y1+c1yy˙1+k1yy1+c12p˙12+k12g12tcosβsinα−ϕ12=0,m1z1+c1zz˙1+k1zz1+c12p˙12+k12g12tcosβcosα−ϕ12=0,Ip1θ¨x1+c1θxθ˙x1+r1c12p˙12+k12g12tcosβ=T1,Id1θ¨y1+c1θyθ˙y1+ω1Ip1θ˙z1+k1θyθy1+r1c12p˙12+k12g12tsinβsinα−ϕ12=0,Id1θ¨z1+c1θzθ˙z1−ω1Ip1θ˙y1+k1θzθz1+r1c12p˙12+k12g12tsinβcosα−ϕ12=0.The dynamic functions of Gear 2 are shown as follows:(23)m2x¨2+c2xx˙2+k2xx2+c12p˙12+k12g12tsinβ−∑i=35c2ip˙2i+k2ig2itsinβ=0,m2y¨2+c2yy˙2+k2yy2−c12p˙12+k12g12tcosβsinα−ϕ12,+∑i=35c2ip˙2i+k2ig2itcosβsinα−ϕ12=0,m2z¨2+c2zz˙2+k2zz2−c12p˙12+k12g12tcosβcosα−ϕ12,+∑i=35c2ip˙2i+k2ig2itcosβcosα−ϕ12=0,Ip2θ¨2+c2θxθ˙x2−r2c12p˙12+k12g12tcosβ+r2∑i=35c2ip˙2i+k2ig2itcosβ=0,Id2θ¨y2+c2θyθ˙y2+ω2Ip2θ˙z2+k2θyθy2+r2c12p˙12+k12g12tsinβsinα−ϕ12,+r2∑i=35c2ip˙2i+k2ig2itsinβsinα−ϕ12=0,Id2θ¨z2+c2θzθ˙z2−ω2Ip2θ˙y2+k2θzθz2+r2c12p˙12+k12g12tsinβcosα−ϕ12,+r2∑i=35c2ip˙2i+k2ig2itsinβcosα−ϕ12=0.The dynamic functions of Gear 3 are shown as follows:(24)m3x¨3+c3xx˙3+k3xx3+c23p˙23+k23g23tsinβ=0,m3y¨3+c3yy˙3+k3yy3−c23p˙23+k23g23tcosβsinα−ϕ23=0,m3z¨3+c3zz˙3+k3zz3−c23p˙23+k23g23tcosβcosα−ϕ23=0,Ip3θ¨x3+c3θxθ˙x3−r3c23p˙23+k23g23tcosβ=T3,Id3θ¨y3+c3θyθ˙y3+ω3Ip3θ˙z3+k3θyθy3+r3c23p˙23+k23g23tsinβsinα−ϕ23=0,Id3θ¨z3+c3θzθ˙z3−ω3Ip3θ˙y3+k3θzθz3+r3c23p˙23+k23g23tsinβcosα−ϕ23=0.The dynamic functions of Gear 4 are shown as follows:(25)m4x¨4+c4xx˙4+k4xx4+c24p˙24+k24g24tsinβ=0,m4y¨4+c4yy˙4+k4yy4−c24p˙24+k24g24tcosβsinα−ϕ24=0,m4z¨4+c4zz˙4+k4zz4−c24p˙24+k24g24tcosβcosα−ϕ24=0,Ip4θ¨x4+c4θxθ˙x4−r4c24p˙24+k24g24tcosβ=T4,Id4θ¨y4+c4θyθ˙y4+ω4Ip4θ˙z4+k4θyθy4+r4c24p˙24+k24g24tsinβsinα−ϕ24=0,Id4θ¨z4+c4θzθ˙z4−ω4Ip4θ˙y4+k4θzθz4+r4c24p˙24+k24g24tsinβcosα−ϕ24=0.The dynamic functions of Gear 5 are shown as follows:(26)m5x¨5+c5xx˙5+k5xx5+c25p˙25+k25g25tsinβ=0,m5y¨5+c5yy˙5+k5yy5−c25p˙25+k25g25tcosβsinα−ϕ25=0,m5z¨5+c5zz˙5+k5zz5−c25p˙25+k25g25tcosβcosα−ϕ25=0,Ip5θ¨x5+c5θxθ˙x5−r5c25p˙25+k25g25tcosβ=T,5Id5θ¨y5+c5θyθ˙y5+ω5Ip5θ˙z5+k5θyθy5+r5c25p˙25+k25g25tsinβsinα−ϕ25=0,Id5θ¨z5+c5θzθ˙z5−ω5Ip5θ˙y5+k5θzθz5+r5c25p˙25+k25g25tsinβcosα−ϕ25=0,where Ti is the load on the Gear i; mi, Ipi, and Idi are the mass, polar moment of inertia, and diameter moment of inertia of Gear i, respectively; pij 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 gijt is the clearance function expressed as(27)gijt=pij−bij,pij>bij,0,pij≤bij,pij+bij,pij<−bij,where 2bij is the clearance length between Gear i and j.The space vector of the system totally has 30 DOFs and can be expressed as(28)X=x1,y1,z1,θx1,θy1,θz1,⋯,x5,y5,z5,θx5,θy5,θz5T.Equations (22)-(26) can be written as(29)MX¨+Cs+Cm+HX˙+KsX+KmGp=F,where M, Cs, Ks, Cm, and Km 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; Gp is a tooth gap function related to p; and F is load vector.Through defining thatxθi=riθxi, yθi=riθyi, and zθi=riθ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 bc, and define that(30)τ=ωct,X¯ij=Xijbc,e¯ij=eijbc,b¯ij=bijbc,X¯=x¯1,y¯1,z¯1,x¯θ1,y¯θ1,z¯θ1,⋯,x¯5,y¯5,z¯5,x¯θ5,y¯θ5,z¯θ5T,dX¯dτ=X¯˙,d2X¯dτ2=X¯¨,S=diag1,1,1,r1,r1,r1,⋯,1,1,1,r5,r5,r5.(25) can be converted into proper dimensionless indexes and expressed as(31)X¯¨+C¯s+C¯m+H¯X¯˙+K¯sX¯+K¯mG¯p¯=F¯,where C¯s=SM−1CsS−1/ωc， C¯m=SM−1CmS−1/ωc， K¯s=SM−1KsS−1/ωc， K¯m=SM−1KmS−1/ωc2， H¯s=SM−1HS−1/ωc, and F¯=SM−1/bcωc2F。The backlash function of each gear pair can be defined as(32)G¯ijp¯ij=X¯ij−14SijBijb¯ij,p¯ij>b¯ij,0,p¯ij≤b¯ij,X¯ij+14SijBijb¯ij,p¯ij<−b¯ij.


## 3. Self-Excited Vibration Analysis of Gear-Bearing System

### 3.1. Vibration Behaviors of the Gear System by Self-Excited Vibration

The main parameters of gear pairs in Figure2 are shown in Table 1. The backlash and meshing error of the gear pairs are 100 μm and 20 μm, respectively.Table 1
Parameters of gears.
Gear 1Gear 2Gear 3Gear 4Gear 5Modulus (mm)2.5Face width (mm)45Pressure angle (deg)20Tooth number80193413325Radius of basic circle (mm)191.6462.398.279.159.9Mass (kg)13503.82.51.7Polar inertia moment (kgm2)0.071.60.0070.0040.003Inertia moment (kgm2)0.040.80.0050.0030.002Rated loads on the three output shafts areT3 = 24 Nm, T4 = 20 Nm, and T5 = 18 Nm. The following three simulated conditions are depicted.Condition 1.
Rated load with constant bearing damping and stiffness at 4000 r/min.Condition 2.
Rated load and variable bearing damping and stiffness calculated by Equations (13) and (17).Condition 3.
1.5 times the rated load with variable bearing damping and stiffness calculated by Equations (13) and (17).
The 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.
The nonlinear vibration characteristics of the gear system under the impact of speed are shown in Figure5(a). Before the input speed reaches 1600r/min, the gear vibration is harmonic. Then 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 Figure5(b). The stiffness and the damping of the bearing are positively correlated with the speed of the bearing, according to Equations (13) and (17). The 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. These 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.
The calculated result of the gear system under 1.5 times rated load differs significantly from that shown in Figure5(c) under rated load. With a heavy load, the non-linear speed interval decreases noticeably. The 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. The gear vibration signal alternates between two periods’ motion and on period motion throughout the entire input speed range. The comparison results show that the load on the gear system has a significant impact on the gear system’s nonlinear vibration.
Figure6 depicts the frequency-domain characteristics of the vibration displacement of Gear 4 along y-direction under the three conditions. There is no discernible difference between Condition 1 and Condition 2 in terms of frequency-domain features. Except for the meshing frequency fm and its double frequency of 2fm, the subfrequency fm/2 begins to appear from 1500r/min to 3500r/min. Fm/4 and fm/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. The 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. The frequency-domain features of Condition 3 are quite different from those of the first two conditions. Except for the meshing frequency fm and its double frequency of 2fm, only the subfrequency fm/2 can be found between 1700r/min to 3000r/min, and no other subfrequencies are presented. The peak value of subfrequency fm/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 fm compared to Conditions 1 and 2.
The influence of variable bearing coefficients on the analysis results is compared using a correlation analysis between Condition 1 and Condition 2. Figure7 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. The correlation coefficient of frequency-domain 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.
The vibration signals of y4 under Condition 1 at different rotational speeds are shown in Figure8 to compare and analyze the change of gear vibration with input speed. The 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. The vibration will act directly on the supporting bearing, producing disturbed vibration force. The 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.Figure 5
Bifurcation diagram of y4. (a) Condition 1. (b) Condition 2. (c) Condition 3.
(a)(b)(c)Figure 6
Frequency-domain analysis of y4. (a) Condition 1. (b) Condition 2. (c) Condition 3.
(a)(b)(c)Figure 7
Correlation coefficient of vibration results between Conditions 1 and 2.Figure 8
Time-domain signals of y4 under Condition 1. (a) 1000r/min. (b) 2000r/min. (c) 2500r/min. (d) 4000r/min.
(a)(b)(c)(d)Figure 9
Time-domain analysis of meshing force between Gear 2 and Gear 4. (a) Peak to peak value. (b) Standard deviation.
(a)(b)


### 3.2. Meshing Force of the Gear System by Self-Excited Vibration

Figure10 depicts the frequency-domain features of meshing force between Gear 2 and Gear 4 at various speeds.Figure 10
Frequency-domain analysis of meshing force between Gear 2 and Gear 4. (a) Condition 1. (b) Condition 2. (c) Condition 3.
(a)(b)(c)There is also no discernible difference between the frequency-domain features of meshing force under Conditions 1 and 2. In comparison to the results in Figure6, the peaks of meshing frequency are unaffected by changing speed, and they are significantly greater than the peaks of subfrequency fm/2 even in the non-linear speed region. fm/4 and fm/8 can also be found from 1500 r/min to 3500 r/min, but their peaks of these subfrequencies are far smaller than fm/2. The frequency-domain features of meshing force differ significantly from those shown in Figure 6(c). Except for the meshing frequency fm 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.Figure11 compares and analyzes the mean, peak to peak, and RMS of the meshing force in these three conditions. The mean value of the three groups of meshing forces varies little as the rotational speed increases. The values under Condition 3 are roughly 1.5 times those of Conditions 1 and 2, which are proportional to the load. The peak to peak values and RMS values in the nonlinear speed region are obviously higher than the other speeds at rated load. The 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.Figure 11
Time-domain analysis of meshing force between Gear 2 and 4. (a) Mean value. (b) Peak to peak value. (c) RMS.
(a)(b)(c)The 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 meshing force will change. In general, increasing the load and decreasing the stiffness of the bearing appropriately can improve the stability of meshing.


## 3.1. Vibration Behaviors of the Gear System by Self-Excited Vibration

The main parameters of gear pairs in Figure2 are shown in Table 1. The backlash and meshing error of the gear pairs are 100 μm and 20 μm, respectively.Table 1
Parameters of gears.
Gear 1Gear 2Gear 3Gear 4Gear 5Modulus (mm)2.5Face width (mm)45Pressure angle (deg)20Tooth number80193413325Radius of basic circle (mm)191.6462.398.279.159.9Mass (kg)13503.82.51.7Polar inertia moment (kgm2)0.071.60.0070.0040.003Inertia moment (kgm2)0.040.80.0050.0030.002Rated loads on the three output shafts areT3 = 24 Nm, T4 = 20 Nm, and T5 = 18 Nm. The following three simulated conditions are depicted.Condition 1.
Rated load with constant bearing damping and stiffness at 4000 r/min.Condition 2.
Rated load and variable bearing damping and stiffness calculated by Equations (13) and (17).Condition 3.
1.5 times the rated load with variable bearing damping and stiffness calculated by Equations (13) and (17).
The 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.
The nonlinear vibration characteristics of the gear system under the impact of speed are shown in Figure5(a). Before the input speed reaches 1600r/min, the gear vibration is harmonic. Then 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 Figure5(b). The stiffness and the damping of the bearing are positively correlated with the speed of the bearing, according to Equations (13) and (17). The 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. These 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.
The calculated result of the gear system under 1.5 times rated load differs significantly from that shown in Figure5(c) under rated load. With a heavy load, the non-linear speed interval decreases noticeably. The 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. The gear vibration signal alternates between two periods’ motion and on period motion throughout the entire input speed range. The comparison results show that the load on the gear system has a significant impact on the gear system’s nonlinear vibration.
Figure6 depicts the frequency-domain characteristics of the vibration displacement of Gear 4 along y-direction under the three conditions. There is no discernible difference between Condition 1 and Condition 2 in terms of frequency-domain features. Except for the meshing frequency fm and its double frequency of 2fm, the subfrequency fm/2 begins to appear from 1500r/min to 3500r/min. Fm/4 and fm/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. The 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. The frequency-domain features of Condition 3 are quite different from those of the first two conditions. Except for the meshing frequency fm and its double frequency of 2fm, only the subfrequency fm/2 can be found between 1700r/min to 3000r/min, and no other subfrequencies are presented. The peak value of subfrequency fm/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 fm compared to Conditions 1 and 2.
The influence of variable bearing coefficients on the analysis results is compared using a correlation analysis between Condition 1 and Condition 2. Figure7 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. The correlation coefficient of frequency-domain 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.
The vibration signals of y4 under Condition 1 at different rotational speeds are shown in Figure8 to compare and analyze the change of gear vibration with input speed. The 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. The vibration will act directly on the supporting bearing, producing disturbed vibration force. The 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.Figure 5
Bifurcation diagram of y4. (a) Condition 1. (b) Condition 2. (c) Condition 3.
(a)(b)(c)Figure 6
Frequency-domain analysis of y4. (a) Condition 1. (b) Condition 2. (c) Condition 3.
(a)(b)(c)Figure 7
Correlation coefficient of vibration results between Conditions 1 and 2.Figure 8
Time-domain signals of y4 under Condition 1. (a) 1000r/min. (b) 2000r/min. (c) 2500r/min. (d) 4000r/min.
(a)(b)(c)(d)Figure 9
Time-domain analysis of meshing force between Gear 2 and Gear 4. (a) Peak to peak value. (b) Standard deviation.
(a)(b)


## 3.2. Meshing Force of the Gear System by Self-Excited Vibration

Figure10 depicts the frequency-domain features of meshing force between Gear 2 and Gear 4 at various speeds.Figure 10
Frequency-domain analysis of meshing force between Gear 2 and Gear 4. (a) Condition 1. (b) Condition 2. (c) Condition 3.
(a)(b)(c)There is also no discernible difference between the frequency-domain features of meshing force under Conditions 1 and 2. In comparison to the results in Figure6, the peaks of meshing frequency are unaffected by changing speed, and they are significantly greater than the peaks of subfrequency fm/2 even in the non-linear speed region. fm/4 and fm/8 can also be found from 1500 r/min to 3500 r/min, but their peaks of these subfrequencies are far smaller than fm/2. The frequency-domain features of meshing force differ significantly from those shown in Figure 6(c). Except for the meshing frequency fm 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.Figure11 compares and analyzes the mean, peak to peak, and RMS of the meshing force in these three conditions. The mean value of the three groups of meshing forces varies little as the rotational speed increases. The values under Condition 3 are roughly 1.5 times those of Conditions 1 and 2, which are proportional to the load. The peak to peak values and RMS values in the nonlinear speed region are obviously higher than the other speeds at rated load. The 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.Figure 11
Time-domain analysis of meshing force between Gear 2 and 4. (a) Mean value. (b) Peak to peak value. (c) RMS.
(a)(b)(c)The 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 meshing force will change. In general, increasing the load and decreasing the stiffness of the bearing appropriately can improve the stability of meshing.


## 4. Conclusion

The 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.The nonlinear vibration might appear as the meshing frequency increases due to the effect of nonlinear parameters such as meshing clearance and transmission error. The 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.The peak to peak and standard deviation values of the vibration data are substantially bigger in the nonlinear speed area than in the other speeds. This means that the disturbance force acting on the supporting bearing would be raised as well. The disrupted force in the nonlinear speed area would plainly diminish as the load increased. The 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. The large load of the gear-bearing system can effectively reduce nonlinear vibration and increase meshing stability.



---

*Source: 4132673-2022-03-25.xml*