# Fluid-Structure Interaction Analysis on Turbulent Annular Seals of Centrifugal Pumps during Transient Process

Authors: Qinglei Jiang; Lulu Zhai; Leqin Wang; Dazhuan Wu
Journal: Mathematical Problems in Engineering
 (2011)

---

## Abstract

The current paper studies the influence of annular seal flow on the transient response of centrifugal pump rotors during the start-up period. A single rotor system and three states of annular seal flow were modeled. These models were solved using numerical integration and finite difference methods. A fluid-structure interaction method was developed. In each time step one of the three annular seal models was chosen to simulate the annular seal flow according to the state of rotor systems. The objective was to obtain a transient response of rotor systems under the influence of fluid-induced forces generated by annular seal flow. This method overcomes some shortcomings of the traditional FSI method by improving the data transfer process between two domains. Calculated results were in good agreement with the experimental results. The annular seal was shown to have a supportive effect on rotor systems. Furthermore, decreasing the seal clearance would enhance this supportive effect. In the transient process, vibration amplitude and critical speed largely changed when the acceleration of the rotor system increased.


---

## Body

## 1. Introduction

In multistage centrifugal pumps, annular seal passages, such as sealing rings and balance drums, can provide fluid-induced forces, which play an important role in determining the rotor dynamic stability of a centrifugal pump [1, 2]. The flow in these annular seal passages is generally called the annular seal flow. The annular seal flow in pumps can generate reaction forces to support rotating shafts similar to journal bearings. A completely developed turbulent flow usually exists in annular seals than a laminar flow in journal bearings. In the past, the influence of the annular seal flow on rotor systems was usually studied indirectly in terms of dynamic coefficients. In the 1970s, Hirs presented a turbulent lubrication theory called the “bulk flow theory” in his doctoral thesis [3]. This result was published in the Journal of the American Society of Mechanical Engineers [4]. Childs [5] developed a method to obtain dynamic coefficients (K, k, M, C, c) based on the “bulk flow theory” using perturbation analysis. This method has been widely used in numerical and experimental studies of the pump shaft stability analysis [6, 7].The start-up process of the centrifugal pump has gained increasing attention recently. In many cases, centrifugal pumps must experience a quick-start process due to its function demand [8, 9]. In this process, the rotor response involves simultaneous free and forced vibrations, a typical transient response very much different from the normal steady response [10]. The dynamic coefficients method mentioned above is not applicable in such situations. First, this method disregards the impeller difference eccentricity, whereas, in the quick-start process, the eccentricity rapidly increases. Second, the amplitude and phase angle of the impellers change rapidly with the increase in rotating speed. The dynamic coefficients method cannot consider nonlinear characteristics of fluid-induced forces and inertia effects of rotor systems due to its linear simplification [5]. Given these disadvantages of the traditional methods, studying the influence of the annular seal on rotor systems in the quick-start process using other methods is beneficial. The fluid-structure interaction method is a good choice because it can directly transfer data between two domains.  Richardet and Lornage [11, 12] studied shaft-fluid coupling problems using two similar FSI methods. The finite element method and finite difference method were used to solve structural analysis and fluid film analysis. Afterwards, the data of the structural and fluid meshes were transferred between each other based on an interfacing grid concept.To conduct the FSI analysis, a proper annular seal model must be selected first. Although Richardet and Lornage [11, 12] used the classical laminar model to simulate a small clearance flow in the FSI analysis, the laminar model is not applicable in the present study. The reason is that in centrifugal pumps, the annular seal flow is not maintained in the laminar state but grows from a laminar state to transition regime state and finally a turbulent state with the increase of the rotor rotating speed. A similar phenomenon was presented by Taylor in 1922 [13]. He found that the axial vertex (Taylor vertex) appears between two concentric cylinders when the internal cylinder starts to rotate. This vertex will break if the angular velocity continues to increase, indicating that turbulent flow occurs in the angular clearance. Laminar flow can be easily simulated by classical Reynolds equations [14]; thus, much attention has been given to the transition regime and turbulent flow. Con et al. [3, 15, 16] proposed three turbulent theories on the angular clearance flow in the 1970s. These three theories found a way to express turbulent shear stress using the Prandtl mixing length theory, law of wall, and experimental results. These expressions of turbulent stress were then submitted to the momentum equations of clearance flow, leading to turbulent lubrication equations, including turbulent factors. The angular clearance flow is usually in the transition regime state between the laminar and turbulent states when the Reynolds number is between 700 and 2000. The transition regime flow cannot be directly simulated by classical Reynolds equations or turbulent theories. In the late 1980s [17, 18], numerical and experimental studies were carried out for the transition regime flow. Roberts and Mason studied these characteristics in 1986 [19] and then proposed a transition flow model. The transition flow model is similar in form to the turbulent model except for the difference in turbulent factor expressions.In traditional FSI methods, fluid mesh is usually much finer than the structure mesh. Thus data transfer between two domains is usually difficult and inaccurate due to missing data. In the current paper, the rotor system was modeled based on the node-element method instead of the traditional FE method. The deformation of the impeller itself is ignored in the node-element method. Because only the radial displacement and the deflection angle of the impeller are useful in the rotor dynamic analysis and the deformation of the impeller is usually not necessary in FSI calculation, this simplification is accurate enough. Moreover, this simplification can overcome problems on missing data and remarkably reduce computing time by replacing the impeller with a node. In the current paper, a typical single rotor system was modeled and coupled with an annular seal flow model. The transient response in the start-up process was then obtained using the FSI method presented. Fluid-induced forces showed great influence on the rotor dynamic analysis. The results can be helpful in the improvement of the vibration characteristics of the centrifugal pump during the start-up period.


## 2. Models of Rotors and Turbulent Flow

### 2.1. Dynamic Model of a Rotor System

The dynamic model of a single rotor system was first established. System equations of motion were obtained under the following assumptions. (a) A rotating shaft is supported simply at both ends and carries one disk at the middle of span. The mass of the shaft is negligible. (b) The angular acceleration of the shaft is constant with time. Figure1 shows the cross-section of the shaft at the disk position, including the geometrical center S, the gravitational center of the disk, eccentricity of the disk e, static steady position of geometrical center O′, and deflection of shaft r.Figure 1
Position of the disk.The geometrical centerS and gravitational center have the following relations:(2.1)xc=x+ecos⁡ψ,yc=y+esin⁡ψ.
After derivation, (2.2)ẋc=ẋ-eψ̇sin⁡ψ,ẏc=ẏ+eψ̇cos⁡ψ,ẍc=ẍ-eψ̈sin⁡ψ-eψ̇2cos⁡ψ,ÿc=ÿ+eψ̈cos⁡ψ-eψ̇2sin⁡ψ.The reaction force and damping force generated by the shaft and applied on the disk can be expressed as follows:(2.3)Fx=-kx,Fy=-ky,Rx=-c1ẋc,Ry=-c1ẏc,
where k is the stiffness of shaft and c1 is the lateral motion damping coefficient of the disk, which is usually approved by the experiments.Based on the theorem of motion of centre of mass(2.4)mẍc=Fx+Rx,mÿc=Fy+Ry-mg.
The equations of motion of the geometrical center S can be written as follows: (2.5)mẍ+c1(ẋ-eψ̇sin⁡ψ)+kx=meψ̇2cos⁡ψ+meψ̈sin⁡ψ,mÿ+c1(ẏ+eψ̇cos⁡ψ)+ky=meψ̇2sin⁡ψ-meψ̈cos⁡ψ-mg.
The system equations can be simplified as follows after substituting x=X, y=Y-mg/k:(2.6)mẌ+c1(Ẋ-eψ̇sin⁡ψ)+kX=meψ̇2cos⁡ψ+meψ̈sin⁡ψ,mŸ+c1(Ẏ+eψ̇cos⁡ψ)+kY=meψ̇2sin⁡ψ-meψ̈cos⁡ψ.
Suppose X=rcosφ, Y=rsinφ, ψ=φ+θ; then(2.7)r̈-rφ̇2+kmr+c1m(ṙ-eψ̇sin⁡θ)=eψ̇2cos⁡θ+eψ̈sin⁡θ,rφ̈+2ṙφ̇+c1m(rφ̇+eψ̇cos⁡θ)=eψ̇2sin⁡θ-eψ̈cos⁡θ.
The dimensionless form of system equations can be obtained after substituting ωn2=k/m, ρ=r/e, D*=c/(2mωn): (2.8)ρ̈-ρφ̇2+2D*(ρ̇-ψ̇sin⁡θωn)+ρ=ωn2(ψ̇2cos⁡θ+ψ̈sin⁡θ),ρφ̈+2ρ̇φ̇+2D*(ρφ̇+ψ̇cos⁡θωn)=ωn2(ψ̇2sin⁡θ-ψ̈cos⁡θ).
For convenience of calculation, system equations can be expressed as follows:(2.9)ρ̈=-(1-φ̇2)ρ-2D*(ρ̇-ψ̇sin⁡θ)+(ψ̇2cos⁡θ+ψ̈sin⁡θ),φ̈=-2ρ̇φ̇ρ-2D*(ρφ̇+ψ̇cos⁡θ)ρ+(ψ̇2sin⁡θ-ψ̈cos⁡θ)ρ.
The solving process of these system equations of motion is shown as follows.(a) Using the relationshipsκ=ω0/ωn, τ=ωnt, λ=ψ̈/2ωn2 and supposing the single rotor system is initially rotating at a uniform speed when τ=0, ψ̇=ω0, the equations of motion of the disk can be expressed by the following equation:(2.10)mr̈+c1ṙ+kr=F0eiω0t,
where F0=meω02.The steady response (amplitude and phase angle) of this single rotor system can be obtained [20]: (2.11)r0=F0/k(1-s2)2+(2ξs)2,θ0=arctan2ξs1-s2,
where ξ=c/2km, s=ω0/ωn.The dimensionless form of response is(2.12)ρ0=κ2(1-κ2)2+(2κD*)2,θ0=arctan2κD*1-k2.
Initially, ψ̇0=κ, φ̇0=κ, ρ̇0=0, ψ0=0, φ0=-θ0.(b) Equation (2.9) is solved after substituting the initial values. Then, ρ̈0 and φ̈0 are obtained.(c) After all initial conditions become available, the next time step starts by calculating the variables as follows:(2.13)ρ̇1=ρ̇0+Δτρ̈0,ρ1=ρ0+Δτ(ρ̇0+ρ̇1)2,φ̇1=φ̇0+Δτφ̈0,φ1=φ0+Δτ(φ̇0+φ̇1)2.(d)  Solving (2.9) by substituting ψ̇i=κ+2λτ, ψ̈i=κτ+λτ2, and θi=ψi-φi, then ρ̈1 and φ̈1 are obtained.(e) Repeat steps (b)–(d).Equation (2.9) is solved using the numerical time-marching method. For linear systems, the time step verifies <1/10f, where f is the highest frequency of the shaft rotation. However, due to the highly nonlinear property of the fluid-induced force, the time step used in the numerical computation is =10–6 s, which is significantly lower than the critical limit (10–3 s). This indicates that about 2000 calculations are made per rotation.


### 2.2. Turbulent Flow Model of the Annual Seal

With the increase of the rotor angular velocity, the small clearance flow in the centrifugal pump will grow from a laminar state to a transition regime and finally to a turbulent state. Proper models that simulate all three states flow should be selected to study the influence of the annular seal flow on the transient response of rotor systems through the start-up process. Based on the three turbulent theories and transition regime model presented by con and Roberts et al. [3, 15, 16, 19] numerical results were obtained and then compared with the experimental results.Figure2 is the schematic of the annular seal flow channel cross-section. The internal boundary has eccentricity e. Figure 2 shows the angular velocity of the internal and external boundaries ωj, ωb; the radius of the internal and external boundaries Rb and Rj, the external boundary center  Ob, the internal boundary center Oj, the eccentricity of internal boundary e′, and the local clearance of flow channel h. Suppose there is a laminar flow in this channel; then the Reynolds laminar equation can be used to simulate the flow [14]: (2.14)∂∂x(h3μ∂p∂x)+∂∂z(h3μ∂p∂z)=6(Uj+Ub)∂h∂x+12∂h∂t.
For the turbulent flow, the following general equation refers to the turbulent theories [3, 15, 16]:(2.15)∂∂x(Gxh3μ∂p∂x)+∂∂z(Gzh3μ∂p∂z)=Uj+Ub2∂h∂x+∂h∂t,
where h=(Rb-Rj)(1+εcosφ), Uj=Rjωj≃Rbωj, Ub=Rbωb. Gx and Gz are the turbulent factors and expressions given in (4.1)–(4.3).Figure 2
Schematic of the small clearance flow cross-section.By using the dimensionless expression as follows, the dimensionless form of (2.15) can be obtained (refer to document) [21]: (2.16)ψ=Rb-RjRb,ε=eRbψ,φ=xRb,z̅=2zB,h̅=1+εcos⁡φ,p̅′=ψ26μωp,λ=B2Rb,∂∂φ(6Gxh̅3∂p̅′∂φ)+6Gzh̅3λ2∂2p̅′∂z̅2=∂h̅2∂φ+1ω∂h̅∂t,
where ω=ωj+ωb and B is the length of the annular seal channel.Addh̅=1+εcosφ to (2.16):(2.17)∂∂φ(Gxh̅3∂p̅′∂φ)+Gzh̅3λ2∂2p̅′∂z̅2=-εsin⁡φ12+ε6ωsin⁡φdδdt+cos⁡φ6ωdεdt.
Suppose (2.18)ω*=(ω-2dδdt),p̅=ψ26μω*p;
then (2.17) can be further simplified: (2.19)∂∂φ(Gxh̅3∂p̅∂φ)+Gyh̅3λ2∂2p̅∂z2=-εsin⁡φ12+cos⁡φ6ω*dεdt.
Add q=(2/ω*)(dε/dt), kx1=1/12Gx, kz1=1/12Gz; then the final form of the turbulent equation can be obtained:(2.20)∂∂φ(h̅3kx1∂p̅∂φ)+h̅3kz1λ2∂2p̅∂z2=-εsin⁡φ+qcos⁡φ.Equation (2.20) has a form similar to the dimensionless equation of the laminar model [21]. However, there is no turbulent factor kx1 and kz1 exists in laminar model. Thus, solving the laminar and turbulent models using one program is convenient. Furthermore, this program can also be used in solving the transition regime flow, except for the difference between the expressions of these models. The turbulent factors (kx1, kz1) presented here are different from the turbulent factors (kx, kz) in (3.3)–(4.3) and those presented by Roberts (kx=kz=6cfRem). The relationships between them are kx1=kx/12, kz1=kz/12.


## 2.1. Dynamic Model of a Rotor System

The dynamic model of a single rotor system was first established. System equations of motion were obtained under the following assumptions. (a) A rotating shaft is supported simply at both ends and carries one disk at the middle of span. The mass of the shaft is negligible. (b) The angular acceleration of the shaft is constant with time. Figure1 shows the cross-section of the shaft at the disk position, including the geometrical center S, the gravitational center of the disk, eccentricity of the disk e, static steady position of geometrical center O′, and deflection of shaft r.Figure 1
Position of the disk.The geometrical centerS and gravitational center have the following relations:(2.1)xc=x+ecos⁡ψ,yc=y+esin⁡ψ.
After derivation, (2.2)ẋc=ẋ-eψ̇sin⁡ψ,ẏc=ẏ+eψ̇cos⁡ψ,ẍc=ẍ-eψ̈sin⁡ψ-eψ̇2cos⁡ψ,ÿc=ÿ+eψ̈cos⁡ψ-eψ̇2sin⁡ψ.The reaction force and damping force generated by the shaft and applied on the disk can be expressed as follows:(2.3)Fx=-kx,Fy=-ky,Rx=-c1ẋc,Ry=-c1ẏc,
where k is the stiffness of shaft and c1 is the lateral motion damping coefficient of the disk, which is usually approved by the experiments.Based on the theorem of motion of centre of mass(2.4)mẍc=Fx+Rx,mÿc=Fy+Ry-mg.
The equations of motion of the geometrical center S can be written as follows: (2.5)mẍ+c1(ẋ-eψ̇sin⁡ψ)+kx=meψ̇2cos⁡ψ+meψ̈sin⁡ψ,mÿ+c1(ẏ+eψ̇cos⁡ψ)+ky=meψ̇2sin⁡ψ-meψ̈cos⁡ψ-mg.
The system equations can be simplified as follows after substituting x=X, y=Y-mg/k:(2.6)mẌ+c1(Ẋ-eψ̇sin⁡ψ)+kX=meψ̇2cos⁡ψ+meψ̈sin⁡ψ,mŸ+c1(Ẏ+eψ̇cos⁡ψ)+kY=meψ̇2sin⁡ψ-meψ̈cos⁡ψ.
Suppose X=rcosφ, Y=rsinφ, ψ=φ+θ; then(2.7)r̈-rφ̇2+kmr+c1m(ṙ-eψ̇sin⁡θ)=eψ̇2cos⁡θ+eψ̈sin⁡θ,rφ̈+2ṙφ̇+c1m(rφ̇+eψ̇cos⁡θ)=eψ̇2sin⁡θ-eψ̈cos⁡θ.
The dimensionless form of system equations can be obtained after substituting ωn2=k/m, ρ=r/e, D*=c/(2mωn): (2.8)ρ̈-ρφ̇2+2D*(ρ̇-ψ̇sin⁡θωn)+ρ=ωn2(ψ̇2cos⁡θ+ψ̈sin⁡θ),ρφ̈+2ρ̇φ̇+2D*(ρφ̇+ψ̇cos⁡θωn)=ωn2(ψ̇2sin⁡θ-ψ̈cos⁡θ).
For convenience of calculation, system equations can be expressed as follows:(2.9)ρ̈=-(1-φ̇2)ρ-2D*(ρ̇-ψ̇sin⁡θ)+(ψ̇2cos⁡θ+ψ̈sin⁡θ),φ̈=-2ρ̇φ̇ρ-2D*(ρφ̇+ψ̇cos⁡θ)ρ+(ψ̇2sin⁡θ-ψ̈cos⁡θ)ρ.
The solving process of these system equations of motion is shown as follows.(a) Using the relationshipsκ=ω0/ωn, τ=ωnt, λ=ψ̈/2ωn2 and supposing the single rotor system is initially rotating at a uniform speed when τ=0, ψ̇=ω0, the equations of motion of the disk can be expressed by the following equation:(2.10)mr̈+c1ṙ+kr=F0eiω0t,
where F0=meω02.The steady response (amplitude and phase angle) of this single rotor system can be obtained [20]: (2.11)r0=F0/k(1-s2)2+(2ξs)2,θ0=arctan2ξs1-s2,
where ξ=c/2km, s=ω0/ωn.The dimensionless form of response is(2.12)ρ0=κ2(1-κ2)2+(2κD*)2,θ0=arctan2κD*1-k2.
Initially, ψ̇0=κ, φ̇0=κ, ρ̇0=0, ψ0=0, φ0=-θ0.(b) Equation (2.9) is solved after substituting the initial values. Then, ρ̈0 and φ̈0 are obtained.(c) After all initial conditions become available, the next time step starts by calculating the variables as follows:(2.13)ρ̇1=ρ̇0+Δτρ̈0,ρ1=ρ0+Δτ(ρ̇0+ρ̇1)2,φ̇1=φ̇0+Δτφ̈0,φ1=φ0+Δτ(φ̇0+φ̇1)2.(d)  Solving (2.9) by substituting ψ̇i=κ+2λτ, ψ̈i=κτ+λτ2, and θi=ψi-φi, then ρ̈1 and φ̈1 are obtained.(e) Repeat steps (b)–(d).Equation (2.9) is solved using the numerical time-marching method. For linear systems, the time step verifies <1/10f, where f is the highest frequency of the shaft rotation. However, due to the highly nonlinear property of the fluid-induced force, the time step used in the numerical computation is =10–6 s, which is significantly lower than the critical limit (10–3 s). This indicates that about 2000 calculations are made per rotation.


## 2.2. Turbulent Flow Model of the Annual Seal

With the increase of the rotor angular velocity, the small clearance flow in the centrifugal pump will grow from a laminar state to a transition regime and finally to a turbulent state. Proper models that simulate all three states flow should be selected to study the influence of the annular seal flow on the transient response of rotor systems through the start-up process. Based on the three turbulent theories and transition regime model presented by con and Roberts et al. [3, 15, 16, 19] numerical results were obtained and then compared with the experimental results.Figure2 is the schematic of the annular seal flow channel cross-section. The internal boundary has eccentricity e. Figure 2 shows the angular velocity of the internal and external boundaries ωj, ωb; the radius of the internal and external boundaries Rb and Rj, the external boundary center  Ob, the internal boundary center Oj, the eccentricity of internal boundary e′, and the local clearance of flow channel h. Suppose there is a laminar flow in this channel; then the Reynolds laminar equation can be used to simulate the flow [14]: (2.14)∂∂x(h3μ∂p∂x)+∂∂z(h3μ∂p∂z)=6(Uj+Ub)∂h∂x+12∂h∂t.
For the turbulent flow, the following general equation refers to the turbulent theories [3, 15, 16]:(2.15)∂∂x(Gxh3μ∂p∂x)+∂∂z(Gzh3μ∂p∂z)=Uj+Ub2∂h∂x+∂h∂t,
where h=(Rb-Rj)(1+εcosφ), Uj=Rjωj≃Rbωj, Ub=Rbωb. Gx and Gz are the turbulent factors and expressions given in (4.1)–(4.3).Figure 2
Schematic of the small clearance flow cross-section.By using the dimensionless expression as follows, the dimensionless form of (2.15) can be obtained (refer to document) [21]: (2.16)ψ=Rb-RjRb,ε=eRbψ,φ=xRb,z̅=2zB,h̅=1+εcos⁡φ,p̅′=ψ26μωp,λ=B2Rb,∂∂φ(6Gxh̅3∂p̅′∂φ)+6Gzh̅3λ2∂2p̅′∂z̅2=∂h̅2∂φ+1ω∂h̅∂t,
where ω=ωj+ωb and B is the length of the annular seal channel.Addh̅=1+εcosφ to (2.16):(2.17)∂∂φ(Gxh̅3∂p̅′∂φ)+Gzh̅3λ2∂2p̅′∂z̅2=-εsin⁡φ12+ε6ωsin⁡φdδdt+cos⁡φ6ωdεdt.
Suppose (2.18)ω*=(ω-2dδdt),p̅=ψ26μω*p;
then (2.17) can be further simplified: (2.19)∂∂φ(Gxh̅3∂p̅∂φ)+Gyh̅3λ2∂2p̅∂z2=-εsin⁡φ12+cos⁡φ6ω*dεdt.
Add q=(2/ω*)(dε/dt), kx1=1/12Gx, kz1=1/12Gz; then the final form of the turbulent equation can be obtained:(2.20)∂∂φ(h̅3kx1∂p̅∂φ)+h̅3kz1λ2∂2p̅∂z2=-εsin⁡φ+qcos⁡φ.Equation (2.20) has a form similar to the dimensionless equation of the laminar model [21]. However, there is no turbulent factor kx1 and kz1 exists in laminar model. Thus, solving the laminar and turbulent models using one program is convenient. Furthermore, this program can also be used in solving the transition regime flow, except for the difference between the expressions of these models. The turbulent factors (kx1, kz1) presented here are different from the turbulent factors (kx, kz) in (3.3)–(4.3) and those presented by Roberts (kx=kz=6cfRem). The relationships between them are kx1=kx/12, kz1=kz/12.


## 3. FSI Analysis Method

To accomplish the FSI, motion equations of a single rotor system and the turbulent lubrication equation must be solved separately. The deflection of the shaft under the influence of the fluid force can then be obtained through data transfer in each time step. Figure3 shows the annular seal flow channel when considering fluid-induced forces. Fluid forces generated by the annular seal flow require the decomposition of the coordinate system to be fixed on the wall. Figure 4 shows the process of data transfer. Details of the calculation process are as follows.Figure 3
Schematic of the FSI theory.Figure 4
Comparison of the numerical results with the experimental results in the transition regime.(a) Solve (2.12) to obtain the initial state of rotor systems based on boundary conditions.(b) Select a proper model to simulate the annular seal flow based on the initial state of rotor systems. The fluid-induced forcesfx and fy can then be obtained from the integration of pressure applied on the annular seal channel. These forces are introduced in (2.11):(3.1)mẍ+c1(ẋ-eψ̇sin⁡ψ)+kx=meψ̇2cos⁡ψ+meψ̈sin⁡ψ+fxcos⁡φ-fysin⁡φ,mÿ+c1(ẏ+eψ̇cos⁡ψ)+ky=meψ̇2sin⁡ψ-meψ̈cos⁡ψ-mg+fxsin⁡φ+fycos⁡φ.
The dimensionless form of (3.1) is as follows: (3.2)r̈-rφ̇2+kmr+c1m(ṙ-eψ̇sin⁡θ)=eψ̇2cos⁡θ+eψ̈sin⁡θ+fxm,rφ̈+2ṙφ̇+c1m(rφ̇+eψ̇cos⁡θ)=eψ̇2sin⁡θ-eψ̈cos⁡θ+fym.
From (3.2) we obtain(3.3)ρ̈=-1(1-φ̇2)ρ-2D*(ρ̇-ψ̇sin⁡θ)+(ψ̇2cos⁡θ+ψ̈sin⁡θ)+fxmeωn2,φ̈=-2ρ̇φ̇ρ-2D*(ρφ̇+ψ̇cos⁡θ)ρ+(ψ̇2sin⁡θ-ψ̈cos⁡θ)ρ+fymeωn2.
Equation (3.3) can then be solved as (2.9).(c) Eccentricity of the diske, angular velocity of rotation dψ/dt, angular velocity of whirling motion dφ/dt, and deflection velocity of the whirling motion dr/dt can be obtained from the results ρ, ψ̇, φ̇, and ρ̇ calculated from (3.2). The four variables e,dψ/dt,dφ/dt, and dr/dt will be introduced into the turbulent lubrication equations. Hirs’ turbulent lubrication theory is used in this paper. By solving (2.6), fluid forces fx and fy are obtained.(d) By repeating (b) and (c), the curves of the numerical solutions are obtained.


## 4. Results and Discussion

### 4.1. Verification of the Rotor Dynamic Model

The dynamic model was verified based on the experimental results presented by Yanabe and Tamura [10] in 1971. Structural parameters of the experimental apparatus are shown in Figure 5. A rotating shaft supported on the ball bearings carries a single disk midway between them. The shaft with a diameter d1 of 17 mm has a span L1 of 520 mm. The disk has a diameter D1 of 120 mm, thickness T of 50 mm, and weight of 4.7 kg. This rotor is driven by a 0.5 sp motor with variable speed and is given a speed varying at 0–90 1/s.Figure 5
Schematic of the experimental apparatus.The experimental results of the response were compared with the numerical solutions in Figure6. There is a slight discrepancy between the experimental value and the calculated one. However, the periodicity of the beats was determined to be in good agreement. Figure 7 compares the rotating speed and the whirling angular velocity of the shaft. When acceleration is constant, the rotating speed increases linearly. Figure 7 also shows that the whirling angular velocity is divided into three regions. In the first region, the rotating speed is almost the same with the whirling angular velocity. However, in the second and third regions, the whirling angular velocity fluctuates around the rotating speed. In the second region, the whirling angular velocity is always smaller than the rotating speed, whereas it is always bigger in the third region at the peak value position. The second region starts from the rotating speed close to the critical speed of the rotor system. The peak value increases with the acceleration of the rotating speed, indicating that the backwards whirling motion may occur when rotor systems pass through the critical speed. Based on the results shown in Figure 7, the phase delay results in Figure 8 are also divided into three regions. In the second region, the phase delay starts to increase steadily because the whirling angular velocity is always smaller than the rotating speed and then fluctuates in the third region. The phase delay angle mainly depends on the damping ratio of the shaft materials. The damping ratio used in this paper was obtained from empirical data, which might not be identical to the damping ratio of the materials used in the experiment apparatus. There is a 4π phase angle difference between the numerical results and the experiment results.Figure 6
Comparison between the numerical and experimental results.Figure 7
Comparison between the angular velocity of whirling motion and rotation.Figure 8
Comparison of the phase delay angle between the numerical and experimental results.


### 4.2. Verification of the Turbulent Flow Model

The turbulent flow model was verified based on the experimental apparatus built in the 1980s by Roberts and Hinton [22]. The annular seal passages used in this investigation were formed by a shaft and a plain cylinder around it. There is a central circumferential groove in the annular seal through which the medium is supplied. The shaft is made of mild steel, has a fixed axis of rotation, and is mounted on two ball bearings. Using this apparatus, Roberts carried out two experiments, namely, experiment a and experiment b. The parameters of the annular seal channel used in these two experiments were different. The structural parameters of this apparatus are shown in Figure 9. The radius of shaft D, radial clearance C, and length of clearance channel L in experiment a are 100, 0.55, and 33 mm, respectively. In experiment b, these three parameters are 120, 0.186, and 30 mm, respectively.Figure 9
Structural parameters of the experimental apparatus.First, a laminar model was used to simulate the annular seal flow in four different conditions. The numerical results were compared with the experiment results in Figures10, 11, and 12. These four conditions correspond to the mean Reynolds numbers 17, 663, 1076, and 48600, respectively. Results were chosen from the pressure distribution at the center position of the annular seal channel in the circumferential direction. Boundary conditions corresponding to the numerical results shown in Figures 10–12 are provided in Table 1.Table 1
Boundary conditions used in Figures10–12.
Figure no.Experiments no.MediumsRotating speed  r/minEccentricity ratioεNormal Reynolds numberRem10btellus22320 r/min0.49/11bwater/0.52107612atellus23/0.7766313aWater/0.9148600Figure 10
Comparison of the laminar results with the experimental results in the laminar state.Comparison of the laminar results with the experimental results in the transition regime.
(a)(b)Figure 12
Comparison of the laminar results with the experimental results in the turbulent state.The numerical results are shown to be in good agreement with the experiment results when the annular seal flow is at the laminar state. However, in the transition and turbulent states, the numerical results are no longer correct. This indicates that the laminar model is only applicable to the laminar state.The turbulent theories presented by Con et al. [15] and the transition regime model presented by Roberts were then used to simulate the annular seal flow in the regime and turbulent states. Based on the finite difference method, the numerical results obtained by solving the equations of the three turbulent theories were compared with the experiment results in Figure 13. The numerical results have the same boundary conditions as those shown in Figure 12. Turbulent coefficients Gx  Gz of the three turbulent lubrication theories developed by Con, Ng-PAN, and Hirs, respectively, are presented in (4.1)–(4.3) [12]. The numerical results show better agreement with the experimental results compared with Figure 12, except that the absolute values of the numerical results are more or less still larger than those of the experimental results. Compared with the other two turbulent theories, Con’s theory is clearly the best in simulating the annular seal flow: (4.1)kx=1Gx=12+0.026Re0.8265,kz=1Gz=12+0.0198Re0.741,(4.2)kx=1Gx=12+0.0136Re0.9,kz=1Gz=12+0.0043Re0.96,(4.3)kx=1Gx=12+0.0687Re0.75,kz=1Gz=12+0.0392Re0.75Figure 13
Comparison of the different turbulent theory results with the experimental results.The annular seal flow is maintained in the transition regime flow state when the Reynolds number is between 500 and 2000. In this regime called the transition regime flow, Roberts presented the turbulent expression of the friction factorcf kx=kz=6cfRem, to simulate the transition flow. The cf was then plotted against the Reynolds number based on the experimental results [19]. The turbulent factors kx, kz of the transition regime used in this paper were obtained by consulting Roberts’ plots. The numerical results of the laminar theory, turbulent theory, and transition regime model were compared with the experimental results, respectively, in Figure 13. The experimental results and the boundary conditions of the numerical computation used in Figure 14 are the same as those used in Figure 11(a). The results of the laminar and turbulent theories have greater absolute values at the hump position in either positive or negative region. By introducing the friction factor cf in simulating the transition regime flow, the numerical results clearly improved compared with the other results based on the laminar or turbulent theories. However, there are still discrepancies between these results, indicating that the friction factor cf presented by Roberts may need to be revised before it is used here.Figure 14
Comparison of the numerical results with the experimental results in the transition regime.


### 4.3. Results of the FSI

Based on the FSI method and related equations presented above, a single rotor system and annular clearance flow were modeled and solved. The transient response of the rotor system in the start-up process was obtained when considering the fluid-induced forces. Considering the nonlinear characteristics of the fluid-induced forces, the time step was set to 1e-6s to guarantee the accuracy of the results.Unless specifically explained, the rotor system has the same structural parameters and boundary conditions as those in Section2. Length of the annular seal is 33 mm, inner radius is 50 mm, outer radius is 50.5 mm, and the medium is water with viscosity of 0.001 Pa.s.Figure15 provides the transient response of the single rotor system. The amplitude of the response decreases, whereas the critical speed increases under the influence of the fluid-induced forces. This proves the existence of the supporting effects of the annular seal flow. Fluid-induced forces change with the increase in rotating velocity. Figure 16 shows that forces in the x direction are smaller by one order of magnitude than those in the y direction. Fluid forces in the y direction clearly have a shape similar to the curve in Figure 15, illustrating that the forces are mainly affected by the radial displacement of disk. The larger shaft deformation leads to higher fluid-induced forces, thus changing the critical speed corresponding to the largest shaft deformation of rotor systems. Moreover, the magnitude of fluid-induced forces is about 1 N, which leads to small effects of the annular seal flow on the vibration characteristics of rotor systems.Figure 15
Transient response of the rotor systems under the influence of the annular seal flow.Figure 16
Fluid-induced forces in thex and y directions.According to the study of Childs et al. [5, 6], annular seals play a supporting role similar to journal bearing in rotor systems. This can explain the increase of critical speed and the decrease of response in Figure 15. However, supporting effect is directly bound up with the annular clearance size and pressure drop of the annular seals. Because pressure drop between the two ends of the annular seals is relatively small in this paper, the increase of critical speed and the decrease of response are not very obvious.Figures17 and 18 show the variation of the whirling angular velocity, the phase delay of the whirling motion with the increase in rotating velocity under the influence of fluid-induced forces, and the comparison between these results when considering the annular seal flow or not. Figure 18 clearly shows that the whirling angular velocity, which considers the fluid-induced forces, is slightly smaller in most regions but larger in several peak value positions than that which does not consider fluid-induced forces. This result leads to two groups of phase delay, as shown in Figure 18, which are almost the same, indicating that the annular seal flow does not affect the damping ratio of rotor systems.Figure 17
Influence of the annular seal flow on the angular velocity of the whirling motion.Figure 18
Influence of the annular seal flow on the phase delay angle.


### 4.4. Transient Effects of the Rotor System

To study the influence of the annular seal flow on a rotor system response in the transient process, three different accelerations of the rotating velocity and size of the annular clearance were solved.First, three different accelerations (λ=0.00026, 0.00052, 0.00104) were selected to determine the influence of starting time on the vibration characteristics of rotor systems. Figure 19 shows that the amplitude clearly decreases when acceleration increases. The amplitude decreases from 37 μm to 18 μm when starting time increases from 0.184 s to 1.84 s, showing a decrease of 51%. The critical speed increases from 268.83 1/s to 303.73 1/s when acceleration increases from 0.00026 to 0.00104. This proves that acceleration increase can clearly postpone the rotating speed corresponding to the critical speed. Here, we apply the term critical speed to the rotating speed at which the maximum deflection occurs in a transient motion.Figure 19
Comparison of the transient response under different accelerations.Subsequently, the influence of fluid-induced forces on the single rotor system was studied by selecting three different clearances while setting the accelerating time to 0.184 s. Clearances are 0 (regardless of the fluid-induced forces), 0.45, and 0.85 mm, respectively. The numerical results are shown in Figure20. The annular seal can clearly reduce the amplitude of a single rotor system. The smaller clearance results in a larger amplitude decrease. The amplitude decreases from 40 μm to 28 μm when the clearance is 0.45 mm compared with the results that disregard the fluid-induced forces. Figure 20 also shows that along with the decrease in clearance size, critical speed markedly increases as well. When the clearance was set to 0.45 mm, the critical speed increased to 314 1/s, which is a 10% growth rate compared with that not considering the annular seal flow.Figure 20
Comparison of the transient response of different seal clearances.


## 4.1. Verification of the Rotor Dynamic Model

The dynamic model was verified based on the experimental results presented by Yanabe and Tamura [10] in 1971. Structural parameters of the experimental apparatus are shown in Figure 5. A rotating shaft supported on the ball bearings carries a single disk midway between them. The shaft with a diameter d1 of 17 mm has a span L1 of 520 mm. The disk has a diameter D1 of 120 mm, thickness T of 50 mm, and weight of 4.7 kg. This rotor is driven by a 0.5 sp motor with variable speed and is given a speed varying at 0–90 1/s.Figure 5
Schematic of the experimental apparatus.The experimental results of the response were compared with the numerical solutions in Figure6. There is a slight discrepancy between the experimental value and the calculated one. However, the periodicity of the beats was determined to be in good agreement. Figure 7 compares the rotating speed and the whirling angular velocity of the shaft. When acceleration is constant, the rotating speed increases linearly. Figure 7 also shows that the whirling angular velocity is divided into three regions. In the first region, the rotating speed is almost the same with the whirling angular velocity. However, in the second and third regions, the whirling angular velocity fluctuates around the rotating speed. In the second region, the whirling angular velocity is always smaller than the rotating speed, whereas it is always bigger in the third region at the peak value position. The second region starts from the rotating speed close to the critical speed of the rotor system. The peak value increases with the acceleration of the rotating speed, indicating that the backwards whirling motion may occur when rotor systems pass through the critical speed. Based on the results shown in Figure 7, the phase delay results in Figure 8 are also divided into three regions. In the second region, the phase delay starts to increase steadily because the whirling angular velocity is always smaller than the rotating speed and then fluctuates in the third region. The phase delay angle mainly depends on the damping ratio of the shaft materials. The damping ratio used in this paper was obtained from empirical data, which might not be identical to the damping ratio of the materials used in the experiment apparatus. There is a 4π phase angle difference between the numerical results and the experiment results.Figure 6
Comparison between the numerical and experimental results.Figure 7
Comparison between the angular velocity of whirling motion and rotation.Figure 8
Comparison of the phase delay angle between the numerical and experimental results.


## 4.2. Verification of the Turbulent Flow Model

The turbulent flow model was verified based on the experimental apparatus built in the 1980s by Roberts and Hinton [22]. The annular seal passages used in this investigation were formed by a shaft and a plain cylinder around it. There is a central circumferential groove in the annular seal through which the medium is supplied. The shaft is made of mild steel, has a fixed axis of rotation, and is mounted on two ball bearings. Using this apparatus, Roberts carried out two experiments, namely, experiment a and experiment b. The parameters of the annular seal channel used in these two experiments were different. The structural parameters of this apparatus are shown in Figure 9. The radius of shaft D, radial clearance C, and length of clearance channel L in experiment a are 100, 0.55, and 33 mm, respectively. In experiment b, these three parameters are 120, 0.186, and 30 mm, respectively.Figure 9
Structural parameters of the experimental apparatus.First, a laminar model was used to simulate the annular seal flow in four different conditions. The numerical results were compared with the experiment results in Figures10, 11, and 12. These four conditions correspond to the mean Reynolds numbers 17, 663, 1076, and 48600, respectively. Results were chosen from the pressure distribution at the center position of the annular seal channel in the circumferential direction. Boundary conditions corresponding to the numerical results shown in Figures 10–12 are provided in Table 1.Table 1
Boundary conditions used in Figures10–12.
Figure no.Experiments no.MediumsRotating speed  r/minEccentricity ratioεNormal Reynolds numberRem10btellus22320 r/min0.49/11bwater/0.52107612atellus23/0.7766313aWater/0.9148600Figure 10
Comparison of the laminar results with the experimental results in the laminar state.Comparison of the laminar results with the experimental results in the transition regime.
(a)(b)Figure 12
Comparison of the laminar results with the experimental results in the turbulent state.The numerical results are shown to be in good agreement with the experiment results when the annular seal flow is at the laminar state. However, in the transition and turbulent states, the numerical results are no longer correct. This indicates that the laminar model is only applicable to the laminar state.The turbulent theories presented by Con et al. [15] and the transition regime model presented by Roberts were then used to simulate the annular seal flow in the regime and turbulent states. Based on the finite difference method, the numerical results obtained by solving the equations of the three turbulent theories were compared with the experiment results in Figure 13. The numerical results have the same boundary conditions as those shown in Figure 12. Turbulent coefficients Gx  Gz of the three turbulent lubrication theories developed by Con, Ng-PAN, and Hirs, respectively, are presented in (4.1)–(4.3) [12]. The numerical results show better agreement with the experimental results compared with Figure 12, except that the absolute values of the numerical results are more or less still larger than those of the experimental results. Compared with the other two turbulent theories, Con’s theory is clearly the best in simulating the annular seal flow: (4.1)kx=1Gx=12+0.026Re0.8265,kz=1Gz=12+0.0198Re0.741,(4.2)kx=1Gx=12+0.0136Re0.9,kz=1Gz=12+0.0043Re0.96,(4.3)kx=1Gx=12+0.0687Re0.75,kz=1Gz=12+0.0392Re0.75Figure 13
Comparison of the different turbulent theory results with the experimental results.The annular seal flow is maintained in the transition regime flow state when the Reynolds number is between 500 and 2000. In this regime called the transition regime flow, Roberts presented the turbulent expression of the friction factorcf kx=kz=6cfRem, to simulate the transition flow. The cf was then plotted against the Reynolds number based on the experimental results [19]. The turbulent factors kx, kz of the transition regime used in this paper were obtained by consulting Roberts’ plots. The numerical results of the laminar theory, turbulent theory, and transition regime model were compared with the experimental results, respectively, in Figure 13. The experimental results and the boundary conditions of the numerical computation used in Figure 14 are the same as those used in Figure 11(a). The results of the laminar and turbulent theories have greater absolute values at the hump position in either positive or negative region. By introducing the friction factor cf in simulating the transition regime flow, the numerical results clearly improved compared with the other results based on the laminar or turbulent theories. However, there are still discrepancies between these results, indicating that the friction factor cf presented by Roberts may need to be revised before it is used here.Figure 14
Comparison of the numerical results with the experimental results in the transition regime.


## 4.3. Results of the FSI

Based on the FSI method and related equations presented above, a single rotor system and annular clearance flow were modeled and solved. The transient response of the rotor system in the start-up process was obtained when considering the fluid-induced forces. Considering the nonlinear characteristics of the fluid-induced forces, the time step was set to 1e-6s to guarantee the accuracy of the results.Unless specifically explained, the rotor system has the same structural parameters and boundary conditions as those in Section2. Length of the annular seal is 33 mm, inner radius is 50 mm, outer radius is 50.5 mm, and the medium is water with viscosity of 0.001 Pa.s.Figure15 provides the transient response of the single rotor system. The amplitude of the response decreases, whereas the critical speed increases under the influence of the fluid-induced forces. This proves the existence of the supporting effects of the annular seal flow. Fluid-induced forces change with the increase in rotating velocity. Figure 16 shows that forces in the x direction are smaller by one order of magnitude than those in the y direction. Fluid forces in the y direction clearly have a shape similar to the curve in Figure 15, illustrating that the forces are mainly affected by the radial displacement of disk. The larger shaft deformation leads to higher fluid-induced forces, thus changing the critical speed corresponding to the largest shaft deformation of rotor systems. Moreover, the magnitude of fluid-induced forces is about 1 N, which leads to small effects of the annular seal flow on the vibration characteristics of rotor systems.Figure 15
Transient response of the rotor systems under the influence of the annular seal flow.Figure 16
Fluid-induced forces in thex and y directions.According to the study of Childs et al. [5, 6], annular seals play a supporting role similar to journal bearing in rotor systems. This can explain the increase of critical speed and the decrease of response in Figure 15. However, supporting effect is directly bound up with the annular clearance size and pressure drop of the annular seals. Because pressure drop between the two ends of the annular seals is relatively small in this paper, the increase of critical speed and the decrease of response are not very obvious.Figures17 and 18 show the variation of the whirling angular velocity, the phase delay of the whirling motion with the increase in rotating velocity under the influence of fluid-induced forces, and the comparison between these results when considering the annular seal flow or not. Figure 18 clearly shows that the whirling angular velocity, which considers the fluid-induced forces, is slightly smaller in most regions but larger in several peak value positions than that which does not consider fluid-induced forces. This result leads to two groups of phase delay, as shown in Figure 18, which are almost the same, indicating that the annular seal flow does not affect the damping ratio of rotor systems.Figure 17
Influence of the annular seal flow on the angular velocity of the whirling motion.Figure 18
Influence of the annular seal flow on the phase delay angle.


## 4.4. Transient Effects of the Rotor System

To study the influence of the annular seal flow on a rotor system response in the transient process, three different accelerations of the rotating velocity and size of the annular clearance were solved.First, three different accelerations (λ=0.00026, 0.00052, 0.00104) were selected to determine the influence of starting time on the vibration characteristics of rotor systems. Figure 19 shows that the amplitude clearly decreases when acceleration increases. The amplitude decreases from 37 μm to 18 μm when starting time increases from 0.184 s to 1.84 s, showing a decrease of 51%. The critical speed increases from 268.83 1/s to 303.73 1/s when acceleration increases from 0.00026 to 0.00104. This proves that acceleration increase can clearly postpone the rotating speed corresponding to the critical speed. Here, we apply the term critical speed to the rotating speed at which the maximum deflection occurs in a transient motion.Figure 19
Comparison of the transient response under different accelerations.Subsequently, the influence of fluid-induced forces on the single rotor system was studied by selecting three different clearances while setting the accelerating time to 0.184 s. Clearances are 0 (regardless of the fluid-induced forces), 0.45, and 0.85 mm, respectively. The numerical results are shown in Figure20. The annular seal can clearly reduce the amplitude of a single rotor system. The smaller clearance results in a larger amplitude decrease. The amplitude decreases from 40 μm to 28 μm when the clearance is 0.45 mm compared with the results that disregard the fluid-induced forces. Figure 20 also shows that along with the decrease in clearance size, critical speed markedly increases as well. When the clearance was set to 0.45 mm, the critical speed increased to 314 1/s, which is a 10% growth rate compared with that not considering the annular seal flow.Figure 20
Comparison of the transient response of different seal clearances.


## 5. Conclusions

A single rotor system and an annular seal flow were modeled and verified. Based on these models, an FSI method was developed in this paper. The influence of the annular seal flow on rotor systems during the start-up process was studied through an FSI analysis. Compared with the traditional rotor dynamic analysis, the amplitude of the transient response decreases, whereas the critical speed increases when considering the annular seal flow. This indicates a supportive effect of the annular seal similar with bearings on the rotors. Moreover, the decrease in clearance size enhances the supportive effects. The amplitude and critical speed of the transient response were greatly influenced by the acceleration of the rotor system in the start-up process. The following are the useful results obtained in this paper.(a)
Con’s turbulent theory and Roberts’ transition regime model can be used to simulate the turbulent and transition regime flows in the annular seal of the centrifugal pump. By considering three different flow states existing in the annular seal, a whole start-up process of the centrifugal pump can be simulated in the calculating process.(b)
The transient response of rotor systems when considering the fluid-induced forces can be obtained by solving two domain equations and then transferring the data between these two domains. In the FSI calculation, rotor systems can be reduced to the node-element model, which does not affect the accuracy of the results.(c)
The fluid-induced forces generated by the annular seal flow can lead to a decrease in the response amplitude and increase in the critical speed of the centrifugal pump. Furthermore, these results can be enhanced by reducing the clearance size. The vibration characteristics of rotor systems improve when acceleration of the rotating speed is increased.



---

*Source: 929574-2011-11-20.xml*