A coupling vibration model of multi-stage pump rotor system based on FEM

Leqin Wang*, Wenjie Zhou**, Xuesong Wei***, Lulu Zhai****, Guangkuan Wu***** *Institute of Process Equipment, Zhejiang University, Hangzhou 310027, China, E-mail: hj_wlq2@ zju.edu.cn **Institute of Process Equipment, Zhejiang University, Hangzhou 310027, China, E-mail: zhouwenjiezwj@zju.edu.cn ***Institute of Process Equipment, Zhejiang University, Hangzhou 310027, China, Faculty of Engineering, Kyushu Institute of Technology, Kitakyushu 804-8550, Japan, E-mail: wxs0773@126.com ****Institute of Process Equipment, Zhejiang University, Hangzhou 310027, China, E-mail: zju_zhailulu@foxmail.com *****Institute of Process Equipment, Zhejiang University, Hangzhou 310027, China, Institute of Water Resources and Hydro-electric Engineering, Xi’an University of Technology, Xi’an 710048, China, E-mail: wuguangkuan@163.com


Introduction
The vibration issue of multi-stage pump rotor system has been one of the main problems in multi-stage pump units.Compared with the single-stage pumps or turbines, the effect of annular seals on the vibration behavior of multi-stage pump rotor system needs additional consideration [1].The multi-stage pump rotor system can be simplified into a rotor-bearing-seal system.As an important factor influencing the vibration status of the rotor system, the dynamic characteristics of bearings and the dynamic vibration of bearing-rotor system have been the target of many researchers [2,3].
Fluid flows in the annular seals can generate reaction force, which increases the stiffness of the rotating shaft system [4].In fact, the seal fluid forces are usually expressed in the form of dynamic coefficients.Childs [5,6] proposed a finite-length method to calculate the dynamic coefficients of annular seals based on the Hirs' turbulent lubrication equations and perturbation method, and this method has been widely applied to the analysis of the stable motion of rotor-bearing-seal system.Kanemori and Iwatsubo [7,8] conducted groups of experiments to research the dynamic coefficients of long annular seals at various rotating speeds, whirling speeds, and pressure drops, they found the experimental results coincided fairly well with Childs' theory.However, the Childs' linear model will be no longer applicable in the case of large disturbance.Muszynska and Bently [9,10] proposed a classic nonlinear seal fluid dynamic force model on the basis of experiments, which could get satisfactory results even if the whirl motion was in large disturbance.The Muszynska nonlinear seal force model was commonly used to study the nonlinear vibration of rotor-bearing-seal system [11,12].
Although some FSI models of the rotor system were proposed in the previous researches, most of them were used to study the transient vibration.In this paper, a novel coupling vibration model of multi-stage pump rotorbearing-seal system for stable motion was proposed.The dynamic coefficients of short annular seals and bearings were calculated by Childs' finite-length method, Hertz contact theory and Elasto-Hydrodynamic Lubrication (EHL), respectively.The calculated results were consistent with the experimental results conducted by a model rotor system.The influence of rotating speed on the dynamic coefficients of short annular seals and bearings was also analyzed.

Ball bearing model
The rollers and raceways of ball bearing will deform when the ball bearing suffers from radial or axial loads.Considering that the double-row self-aligning ball bearings in this paper are not under conditions of high rotating speed, gyroscopic moments and centrifugal forces on the rolling elements are ignored.Similarly, the axial deformations are also ignored due to the small axial force and only the radial load is taken into consideration.Fig. 1 shows the contact deformations between the raceways and rollers under radial force Fr.

Compression deformations between raceways and rollers
The values of the contact deformations can be calculated by Hertz contact theory, which can be calculated by the following form [13]: where ; ∑ρ is curvature sum; Q is contact load; E1 and E2 are elastic moduli of rolling balls and races, respectively; μ1 and μ2 are Poisson's ratio for rolling balls and races, respectively, k is ellipticity; K and L are first and second kind of complete elliptic integrals, respectively.The curvature sum ∑ρ is determined by raceway groove curvature radius factor and geometric structure of the ball bearing and it can be got by Zhou et al [14], thus the sum contact deformations can be calculated as: where  is contact deformation and the subscript 'i' and 'o' denote inner race and outer race, respectively.
From Fig. 2, the load equilibrium equation of the ball bearing can be described as: is the angle of the adjacent rolling balls.Furthermore, due to the symmetrical structure of double-row self-aligning ball bearing and the radial clearance, the largest ball bearing load can be summarized based on massive computation: where z is the number of rolling balls,  is contact angle.
According to the definition of contact stiffness, the Hertz contact stiffness of ball radial bearing can be obtained as: where superscript 'b' denotes bearing.The rolling balls and the raceways are separated by a layer of oil film.Based on EHL, the systemic numerical simulations were carried out by Hamrock and Dowson to research the isothermal elliptical contact issue and propose the film thickness formula [15].The minimum oil film thickness can be expressed as: where ; Rx is equivalent radius of curvature in ball rolling direction, η0 is dynamic viscosity, u is surface velocity, λ is pressure viscosity coefficient.
The total minimum film thickness can be calculated as: where The calculation of lubrication oil film stiffness is similar to Hertz contact stiffness and the lubrication oil film stiffness of double-row self-aligning ball bearing can be expressed according to the definition of oil film stiffness:   Finally, an overall radial stiffness of double-row self-aligning ball bearing can be calculated as:

Annular seal model
The finite-length theory was propose by Childs [5,6], which has been widely applied to calculate the dynamic coefficients of annular seals and can be expressed as: where superscript 's' denotes short annular seal.The values of dynamic coefficients are determined by geometry size and working condition of the short annular seals.

System motion equations
For a rotor-bearing-seal system with n degrees of freedom, the motion state can be described by n independent generalized coordinates uj (j = 1, 2,…, n), which can be expressed by non-conservative Lagrange equation: where T is kinectic energy, V is strain energy, qj are generalized forces.
The motion equations of rigid disc can be obtained by Eq. ( 11) and they have been described in [3]: The motion equations of elastic shaft can be derived from FEM and each element can be expressed as: where superscript 'e' denotes elastic shaft element, e M , e  J and e K are 4 × 4 matrices.Combining the motion equations of the rigid discs and elastic shaft elements, the motion equations of rotor system can be expressed as: where [M1], [J1] and [K1] have the similar structure, and dimensions of these three coefficient matrixes are 2n × 2n; n is the total numbers of nodes.
Eq. ( 14) can be described in the following form: In order to introduce the dynamic coefficients of bearings and short annular seals, a FSI method [1] is used to establish the final vibration model of rotor-bearing-seal system.Assuming the seal is on node j, the displacement of node j can be expressed as:    16) and (10) into Eq.( 15) and the motion equation of the rotor system including seal force can be expressed as: For the general case of multi-stage annular seals and considering the effect of ball bearings, the final motion equation of rotor-bearing-seal system can be written as: and ns, nb are numbers of annular seals and bearings, respectively.
The solving process of coupling vibration mode for multi-stage pump is shown in Fig. 3.

Experimental test facilities
A model rotor system was established to verify the coupling vibration model described above.Fig. 4 shows the experimental model rotor system test facilities.The model rotor system is driven by an AC motor and its rotating speed is controlled by a variable-frequency electric cabinet.The maximum rotating speed of the drive motor and the model rotor are 3,000 r/min and 6,000 r/min, respectively, and the torque is transmitted by a synchronous belt.
In order to observe and record the real-time vibration data for further data processing, a data acquisition system including data acquisition and analysis instrument with multi-channel, notebook computer, dis-placement sensors and power supply is applied.The data acquisition and analysis instrument combines the multi-DSP parallel processing technology, low-noise design technology and high-speed data transmission, which ensures the real-time and accuracy of collected vibration data.Furthermore, the amplitude and phase of the fundamental frequency can be calculated by correlation analysis from the collected data.The model rotor system is the key component in the experiment facilities.As shown in Fig. 5, the model rotor system is designed in segmental structure.There are nine chambers in the system, including five low-pressure chambers and four high-pressure chambers.There are five disks instead of impellers in the corresponding five lowpressure chambers.The only fluid flow passage between high-pressure chamber and low-pressure chamber is the seal channel.The model shaft is supported by two double-row self-aligning ball bearings.In order to produce the pressure drop between the two ends of each annular seal channel, a booster pump is used to supply the highpressure fluid, which is connected with the high-pressure chambers.The displacement sensor is used to the vibration amplitude of the model shaft.
Fig. 6 shows the sealing ring and its inner seal channels.The inlets of the sealing ring are placed in highpressure chamber and the outlets are connected with lowpressure chamber.Each sealing ring has two seal channels between the inlets and outlets.For each sealing ring, the high-pressure fluid from the booster pump is first pumped into the high-pressure chamber and then flows into the low-pressure chambers through the two seal channels.There are totally four sealing rings in the rotor system, which means that there are eight sealing channels, i.

Results and discussion
In order to verify the accuracy of the coupling vibration model mentioned above, the experimental critical speed of model rotor system was compared with the simulation results calculated by the coupling vibration model.The damping coefficients and cross stiffness of ball bearings used in the paper were set to zero.The calculation parameters of bearings and each sealing channel used in the paper are listed in Table 1 and Table 2, respectively (the number of each sealing channel increases successively from the drive-side).Fig. 7 shows the dynamic relationships of Ω with b f k and b t k , it can be seen that the oil film stiffness of right bearing is larger than those of left bearing and both of them decrease rapidly when rotating speed is less than 500 r/min.As rotating speed increases, the oil film stiffness difference of the two bearing decreases and tends to reach the same constant.However, the overall radial stiffness of two bearings approximately decreases linearly as rotating speed increases, this is because the contact stiffness is the main factor influencing the overall radial stiffness.The results imply that the stiffness of the bearings could not be regarded as a constant under different rotating speed and it relates with the rotating speed.According to Fig. 10, It can be obviously seen that the first critical speeds under different pressure drops when Rc = 0.3 mm were calculated.Unlike the calculated results of radius clearances, the first critical speeds increase from 1969.4 r/min to 2605.3 r/min as the pressure drop increases from 0.1 MPa to 0.5 MPa.The calculated results illustrate that greater pressure drop will cause higher critical speed, which is consistent with the calculated results of Jiang et al [1].
Table 3 shows the experimental and calculated first critical speeds under different drop pressures and radius clearances.Comparing the experimental results and the calculated results, it is evident through the Table 3 that the radius clearance and pressure drop play a significant role in the dynamic vibration characteristics of the rotorbearing-seal system.Based on the coupling vibration model developed in the paper, the maximum error percentage and minimum error percentage are 5.5% and 0.2%, respectively and most error percentages are 1%~3%.The reasons of calculated error are mainly due to the ignorance of the eccentricity of shaft and the gyroscopic moments and centrifugal forces of rolling elements, which cause the increase in calculated error of dynamic coefficients of annular seals and ball bearing, respectively.Even so, the error percentages are smaller than those calculated by Jiang et al [1], in which they ignored the effect of rotating speed on the bearings and regarded the dynamic coefficients of the bearings as constants.Even if some other effects on the shaft are ignored in the paper, such as the rotating damping of the fluids and mechanical seals, the calculated results under different operating conditions are accurate enough to meet the engineering needs, which means the coupling vibration model proposed in the paper is feasible.

Conclusions
In order to predict the dynamic vibration characteristics of multi-stage pumps, a coupling vibration model of multi-stage pump rotor system was proposed in this paper.The multi-stage pump rotor system was simplified to a rotor-bearing-seal system and the dynamic coefficients of bearings and short annular seals were considered changing with the rotating speed.Based on FEM and Lagrange equation, the motion equations of the coupling vibration model were established.The dynamic relationships of Ω with the stiffness of bearings and annular seals were calculated and the results shown that the rotating speed was an important factor influencing the dynamic coefficients of bearing and annular seals.In order to verify the numerical results of the coupling vibration model, the first critical speed of a model rotor system was measured.The calculated results implied that smaller radius clearance and greater pressure drop would lead to remarkable Lomakin effect and the coefficients of bearings and annular seals had an important impact on the system.The experimental results of the model rotor showed good agreement with the simulated results.Most error percentages were between 1%~3%.In fact, the stable vibration of rotor system supported by journal bearings and other structure form bearings can also be calculated by the coupling vibration model.

Fr
Fig. 1 Compression deformations between raceways and rollers

Fig. 2
Fig. 2 Radial load distribution where Q0, Q1, Q2,…, Qj are loads on the roller, respectively, r is contact deformation in maximum load direction,  lumped mass,  is rotating speed, Jd, Jp are diametric moment of inertia and polar moment of inertia, respectively, u1 and u2 are the generalized coordinates vectors, q1 and q2 are the generalized forces vectors, superscript 'd' denotes rigid disc.

Fig. 3
Fig. 3 Flow chart of solving process

Fig. 5
Fig. 4 Experimental model rotor system test facilities

Fig. 7 Fig. 8
Fig. 7 The change of rotating speed with: a -oil film stiffness; b -overall radial stiffnessThe short annular sealing channels are divided into three forms according to the length-diameter ratio.The first form includes channel 1, channel 2 and channel 6 to channel 8.The second form includes channel 3. Channel 4 and channel 5 belong to the third form.In Fig.8, we can see that the sealing channel with larger length-diameter ratio has bigger dynamic characteristics.The principal stiffness decreases gradually as the rotating speed increases, but the dynamic characteristics of cross stiffness increases linearly when rotating speed increases.Similarly, The principal damping and cross damping have the same varying tendency.The results also show that the rotating

Table 1
Ball bearing parameters

Table 3
Experimental and calculated results of first critical speed