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

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.


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 quickstart 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.

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.Figure 1 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.
The geometrical center S and gravitational center have the following relations: x c x e cos ψ, y c y e sin ψ.

2.2
The reaction force and damping force generated by the shaft and applied on the disk can be expressed as follows: where k is the stiffness of shaft and c 1 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
The equations of motion of the geometrical center S can be written as follows: m ẍ c 1 ẋ − e ψ sin ψ kx me ψ2 cos ψ me ψ sin ψ, m ÿ c 1 ẏ e ψ cos ψ ky me ψ2 sin ψ − me ψ cos ψ − mg.

2.7
The dimensionless form of system equations can be obtained after substituting ω 2 n k/m, ρ r/e, D * c/ 2mω n :

2.8
For convenience of calculation, system equations can be expressed as follows:

2.9
The solving process of these system equations of motion is shown as follows.a Using the relationships κ ω 0 /ω n , τ ω n t, λ ψ/2ω 2 n 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: where F 0 meω 2 0 .The steady response amplitude and phase angle of this single rotor system can be obtained 20 :
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.

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.
Figure 2 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 R b and R j , the external boundary center O b , the internal boundary center O j , 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 : For the turbulent flow, the following general equation refers to the turbulent theories 3, 15, 16 : x and G z are the turbulent factors and expressions given in 4.1 -4.3 .
By using the dimensionless expression as follows, the dimensionless form of 2.15 can be obtained refer to document 21 : where ω ω j ω b and B is the length of the annular seal channel.Add h 1 ε cos ϕ to 2.16 : 2.17 Suppose then 2.17 can be further simplified: then the final form of the turbulent equation can be obtained: Equation 2.20 has a form similar to the dimensionless equation of the laminar model 21 .However, there is no turbulent factor k x1 and k z1 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 k x1 , k z1 presented here are different from the turbulent factors k x , k z in 3.3 -4.3 and those presented by Roberts k x k z 6c f R em .The relationships between them are k x1 k x /12, k z1 k z /12.

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.Figure 3 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.
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 forces f x and f y can then be obtained from the From 3.2 we obtain

3.3
Equation 3.3 can then be solved as 2.9 .c Eccentricity of the disk e, 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 f x and f y are obtained.
d By repeating b and c , the curves of the numerical solutions are obtained.

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.The experimental results of the response were compared with the numerical solutions in Figure 6.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 Mathematical Problems in Engineering experiment apparatus.There is a 4π phase angle difference between the numerical results and the experiment results.

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.
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 Figures 10,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.
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 G x G z 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: 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 factor c f k x k z 6c f R em , to simulate the transition flow.The c f was then plotted against the Reynolds number based on the experimental results 19 .The turbulent factors k x , k z 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   at the hump position in either positive or negative region.By introducing the friction factor c f 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 c f presented by Roberts may need to be revised before it is used here.

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 Section 2. 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.
Figure 15 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 fluidinduced 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.
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.Figures 17 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 fluidinduced 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.

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 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.

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.

Figure 1 :
Figure 1: Position of the disk.

Figure 2 :
Figure 2: Schematic of the small clearance flow cross-section.

Figure 4 :− e ψ sin θ e ψ2 cos θ e ψ sin θ f x m , r φ 2 ṙ φ c 1 Figure 5 :
Figure 4: Comparison of the numerical results with the experimental results in the transition regime.

Figure 6 :Figure 7 :
Figure 6: Comparison between the numerical and experimental results.

Figure 8 :Figure 9 :
Figure 8: Comparison of the phase delay angle between the numerical and experimental results.

Figure 10 :
Figure 10: Comparison of the laminar results with the experimental results in the laminar state.

Figure 11 :
Figure 11: Comparison of the laminar results with the experimental results in the transition regime.

Figure 12 :
Figure 12: Comparison of the laminar results with the experimental results in the turbulent state.

Figure 13 :
Figure 13: Comparison of the different turbulent theory results with the experimental results.

Figure 14 :Figure 15 :
Figure 14: Comparison of the numerical results with the experimental results in the transition regime.

Figure 16 :
Figure 16: Fluid-induced forces in the x and y directions.

Figure 17 :
Figure 17: Influence of the annular seal flow on the angular velocity of the whirling motion.