A transient bulk flow model with circular whirl motion for rotordynamic coefficients of annular seals

Peng XIA,Zhansheng LIU,Xiangyu YU,Jingming ZHAO

School of Energy Science and Engineering,Harbin Institute of Technology,Harbin 150001,China

1.Introduction

Rotordynamic characteristics of high-performance turbomachinery significantly depend on the hydrodynamic forces developed by annular pressure seals.1Accordingly,a large number of studies have been conducted to theoretically predict and to experimentally measure the hydrodynamic forces in a number of seal configurations under various operating conditions.

The bulk flow model introduced by Hirs2has been widely used to predict the dynamic coefficients for turbulent annular seals.In order to simplify the Navier-Stokes equations and reduce substantial computational costs,the bulk flow theory associated bulk mean flow in the seal clearance with averaged turbulence forces.Childs3firstly used 1D bulk flow model with fluid inertia effect to calculate the rotor dynamic coefficients for smooth annular seals concentric with rotors.San Andres4developed a Computational Fluid Dynamics(CFD)solution for eccentric smooth annular seals.A finite difference scheme was implemented to solve 2D bulk flow model and the numerical solution procedure was based on the Semi-Implicit Method for Pressure Linked Equation(SIMPLE)to accelerate convergence.To improve the accuracy of bulk flow model,simplified turbulence models for different surface textures were developed based on a large number of experimental data in the Refs.5,6Although the turbulent models were proved to be effective in those studies,it was difficult to capture the full nature of turbulence flow based on analogical experiments.In the Refs.7,8,the analytical results and experiment data had good agreement by manually adjusting the friction factors in the turbulence models.Saba et al.9evaluated the friction factors by minimizing leakage derivation between analytical results and experiment data,and the prediction of rotor dynamic coefficients had good agreement with experiment results.However,empirical boundary conditions used in the Refs.7–9were another important error source.As the boundary conditions of bulk flow model also had an influence on the leakage flowrate,the friction factors could not be determined independently in the Ref.9Actually,the boundary conditions such as pressure loss coefficient and pre-swirl ratio could be measured in experiments,but the measurements were very expensive and seldom have been conducted.In the Refs.3–9,the bulk flow model was a static model and a perturbation method was implemented to calculate rotordynamic coefficients with a mathematical simpli fication that the amplitude of rotor whirl motion was much smaller than the clearance thickness,but the amplitude in experiments and engineering hardly fulfilled the perturbation simplification.In addition,there were many other cases where the perturbation condition did not ful fill.In the Refs.10,11,the large-scale dynamic responses of floating ring seals were investigated and the hydrodynamic forces were calculated by the perturbed bulk flow model.

In recent years,CFD has been increasingly used to predict the performance of annular seals.CFD method solves the complete Navier-Stokes equations and is able to predict turbulence flow at a highly detailed level in the solution.Moore and Palazzolo12conducted CFD simulations to investigate the liquid mean velocity of a grooved liquid annular seal.Numerical predictions of the mean velocity at different positions across the seal showed good agreement with experimental results.With the development of computational technology,CFD method was used to calculate rotordynamic coefficients for annular seals.Moore13introduced a steady-state CFD method to calculate rotordynamic coefficients of a labyrinth seal.The method used a steady-state flow field and a rotor whirl motion with small eccentricities.The axisymmetric geometry of the labyrinth seal allowed frequency-independent rotordynamic coefficients to be evaluated with impedance calculation at three different whirl frequencies.Untaroiu et al.14calculated rotordynamic coefficients for a circumferentially grooved liquid seal with the steady-state method.To study pre-swirl effect,the upstream region for the seal was explicitly modeled in CFD study.Chochua and Soulas15developed a 3D transient,dynamic mesh method to calculate dynamic characteristics of a hole-pattern gas seal.Yan et al.16improved the dynamic mesh method with a small circular rotor whirl motion,and the computational costs for a simulation needed 15 days on an Intel Core2 Quad6600 2.4 GHz CPU.Generally,full 3D CFD method demonstrated superior capacity in predicting leakage and rotordynamic coefficients over bulk flow model while the CFD method required long solution time for calculating the dynamic coefficients.In addition,for wall-bounded flow in liquid seals,turbulence intensity in the clearance was much lower than requirements of the wall function in CFD theory and was inevitable to exponentially increase computational costs to solve the wall-bounded flow in 3D CFD analysis.

Several investigations were conducted to take advantage of the two methods.Arghir et al.17modified the turbulence model in bulk flow model with a great deal of CFD results for a typical cell extracted from seal textured surfaces.Migliorini et al.18substituted wall shear stress predicted by steady-state CFD analysis into bulk flow model to calculate rotordynamic coefficients.However,the wall shear stress was only part of the turbulence resistance in seal clearance.

2.Analytical model

2.1.Bulk flow model with arbitrary rotor motion

Under high pressure drop and rotational speed,the flow of low viscosity is characterized by high levels of turbulence.By neglecting velocity variation in the radial(y)direction and incorporating simpli fied turbulence models into the Navier-Stokes equations,the bulk flow continuity and momentum equations(Eqs.(1)–(3))in the axial(z)and circumferential(x)directions are given by the Ref.4within each control volume of 2D flow field in seal clearance.

wheret is the time,w is the axial averaged velocity,u is the circumferential averaged velocity,usis the velocity relative to the stator,uris the velocity relative to the rotor,h is the clearance thickness,p is the fluid pressure,ρ is the fluid density,R is the rotor radius,Ω is the rotor speed,frand fsare Moody’s turbulence formulas at rotor and stator surfaces,n is the friction factor.

To directly consider rotor whirl motion,the squeeze velocity Vsin the control volume is explicitly represented as the time derivative of film thickness.

where θ is the circumferential angle of a control volume,Vradialis the radial velocity and Vtangentialis the tangential velocity in the local coordinate system shown in Fig.1.

Within a trivial time interval Δt,it is assumed that flow variables in Eqs.(1)–(3)lineally vary.For incompressible fluid,the equations are discretized between t and t+Δt time step.

As shown in Fig.1,the 2D control volumes keep the same shape in the local coordinate system at each time step.The node variables of t time step(φt=pt,ut,wt)in the coordinate system of t+Δt time step are obtained by linear interpolation inside the linear discretized control volumes of the t time step.

2.2.Boundary conditions with CFD calibration

Boundary conditions for the velocity and pressure fields at the seal entrance and exit are:

(1)At the seal entrance,the pressure loss effect is modeled aswhere pinis the upstream pressure,penis the seal entrance pressure and δ is the pressure loss coefficient and wenis the mean flow axial velocity at seal entrance.

(2)Depending on upstream conditions,the entrance circumferential velocity is uen= βRΩ,where β is the pres wirl ratio.

(3)At the seal exit,no pressure recovery effect is assumed to exist,pex=pout,where poutis the downstream pressure and pexis the seal exit pressure.

To obtain the above parameters of the boundary conditions,CFD method is used to simulate the steady-state flow ifeld in and around the seal clearance.By extracting weighted-average CFD results at the seal entrance and exit,the parameters are determined,and the friction factor of bulk flow model is further evaluated by minimizing leakage flowrate deviation between CFD method and bulk flow model.

2.3.Numerical solution procedure

The flow field of each time step is discretized with rectangular elements shown in Fig.2,and the distributions of flow variables are obtained by assuming that the variables are linearly distributed inside each element domain Ωe.

where φ=u,w,p;Niis the 2D Lagrange interpolation function and i is the node number in an element.

To obtain the finite element matrix represented by Eq.(9),the Eqs.(5)–(7)are multiplied by the Galerkin test function and integrated over a master isopara metric element.Xu19detailed coordinate transformation and numerical integration procedure with the Gauss-Legendre quadrature formulas.

The sub-matrixes in the left side of the Eq.(9)are 4×4 2D matrixes and the sub-matrixes in the right side are 4×1 1D matrixes.The sub-coefficients in each sub-matrix are given in detail in Appendix A.ue,we,and peare the 4×1 1D matrixes of the node variable ui,wiand piin an element.

The finite element matrixes of all the elements in the flow field are assembled to form an equation system.The alternating linear iteration algorithm is used to solve the equation system and obtain the velocity and pressure distributions in seal clearance.

2.4.Fluid film forces and rotor dynamic coefficients

The leakage flowrate is given by integration of the axial velocity at the seal exit as

At each time step,the hydrodynamic forces in the local coordinate system are given by Eq.(11)and can be transferred into the fixed reference system to obtain the force components in X and Y directions.

where Ftand Frare the hydrodynamic forces in the tangential and radial directions of the local coordinate system,and L is the length of the seal.

When the rotor moves around the seal geometry center with a circular path,the relationship14between the hydrodynamic forces in the local coordinate system and the rotor dynamic coefficients is

where Rwis rotor whirl radius and ω is rotor whirl speed,K is the direct stiffness,k is the cross-coupled stiffness,C is the direct damping,c is the cross-coupled damping,M is the direct inertia coefficient and m is the cross-coupled inertia coefficient.

In order to identify the stiffness,damping,and inertia coefficients in the Eqs.(12)and(13),the hydrodynamic forces are calculated at multiple whirl speeds and a linear least square curve fit is conducted.

3.Results and discussion

3.1.CFD calibration

The smooth annular seal schematically illustrated in Fig.3 is chosen from an experimental seal tested at Texas A&M University.7The seal geometry and operating conditions are shown in Table 1,and the test operating speed is respectively 10200,17400 and 24600 r/min,much lower than the first critical speed(30000 r/min)of the rotor in the test rig.When the seal is concentric with the shaft,a 2D axisymmetric analysis of wall-bounded flow in the clearance with high resolution is conducted to reduce substantial computational costs.Thermal effect is not incorporated,since the specific heat of liquid water is large and leakage flow quickly takes energy dissipation away.Rotor surface is defined as a moving wall with tangential velocity and stator surface is specified as a stationary wall.The walls for the surfaces are set to be nonslip,and centrifugal growth20of the rotor is calculated by

The software ICEM 14.0 is used to generate structured mesh.Under low Reynold number,turbulence flow in the clearance is significantly influenced by the presence of the walls.The k-ε standard turbulence model is used with the enhanced-wall treatment,and the y+values for the first grids near the walls are below 1 with more than 20 nodes inside boundary layers.Total pressure and turbulence quantities are defined at the inlet boundary,and static pressure is specified at the outlet boundary.Discretization scheme is the second-order upwind to guarantee numerical precision.Convergence condition is that the weighted-average pressure,axial velocity and swirl velocity monitored at the seal entrance and exit do not fluctuate with the corresponding scaled residuals of continuous,momentum and turbulence equations reaching 10-8.

In the case of5.52 MPa and 10200 r/min,a gridindependence study is conducted.Fig.4 depicts the grid distributions at the seal entrance with respectively 1.31(Mesh A),2.73(Mesh B)and 5.32(Mesh C)million elements.The velocities and pressure at the seal entrance are selected to be the key variables for the grid-independence study.Fig.5 shows that the pro files of the static pressure have good agreement,and Fig.6 shows the pro files of the axial velocity and the tangential velocity are almost the same with each other.Consequently,Mesh B is convergent and used in the CFD analysis.

As there are no turbulence models which demonstrate superior for turbulent flow in annular seals,the comparison including the k-ε standard model,the k-ε renormalization model,the k-ω model and the SST(Shear Stress Transport)model is conducted.The enhanced wall treatment is used with the k-ε turbulence models,and the k-ω model and the SST model candirectly solve boundary layers.Fig.7 shows the leakage comparison of the turbulence models.The k-ω model underestimates the leakage flowrate, while the SST model overestimates.The leakage flowrates of the k-ε renormalization model and the k-ε standard model are between the k-ω model and the SST model.Compared with the test results,the overall leakage flowrate error for the k-ε standard model is smaller than that of other turbulent models with the maximum derivation below 2.68%.Thus,the k-ε standard model with enhanced wall treatment is used to calibrate bulk flow model.

Table 1 Seal geometry and operating conditions.7

As there is a sudden narrowing at the seal entrance, fluid inertia accelerates the fluid from an upstream stagnant condition to a flow with high axial speed.Fig.8 shows the distributions of pressure and velocity at the seal entrance in the case of Δp=4.14 MPa and Ω =10200 r/min.

In all the cases of this work,the violent variation of the flow field at the entrance is con fined to the entrance region,and the weighted-average variables are measured at the position with 4cl away from the seal entrance.Table 2 shows the parameters of bulk flow model calibrated with CFD results.The pressureloss coefficient decreases with increasing pressure drop.The preswirl ratio increases with increasing rotational speed and decreases with increasing pressure drop.

3.2.Fluid film forces

To minimize numerical errors of the linear curve fit,eleven whirl frequencies ranging from 0.5Ω/60 to 1.5Ω/60 are used to evaluate rotor dynamic coefficients.For each frequency,360 time steps are used during every period and the rotor whirl radius is 0.025cl.

Table 2 Bulk flow parameters with CFD calibration.

Fig.9 shows the reaction forces of four excitation frequencies on the rotor surface in the case of Δp=5.52 MPa and Ω=10200 r/min.The reaction forces FXand FYhave the same frequency with the rotor whirl motion and the amplitudes of the forces increase with increasing the rotor whirl frequency.Table 3 shows the radial force and the tangential force in thelocal coordinate system at the eleven whirl frequencies.With the whirl frequency increasing,the radial force decreases while the tangential force increases.It indicates that the rotor whirl frequency is important for the dynamic stability of the rotor and seal system.Fig.10 shows the data of the forces agree very well with the least square fitted curve.

Table 3 Radial forces and tangential forces on rotor surface with different whirl frequencies.

In addition,to reduce computational costs of calculating rotor dynamic coefficients,the convergence condition for each whirl frequency is that the scaled residuals for the radial force and the tangential force reach 10-6.The computational time for the eleven frequencies is about 1.5 h on an Intel Core i7-4770 3.40 GHz CPU.

3.3.Comparisons and validation

In order to validate the bulk flow model developed in present study,as shown from Fig.11,the rotor dynamic coefficients are compared with the experimental data and analytic results in Ref.7The experimental difference of the direct damping(CXXand CYY)and direct inertia coefficients(MXXand MYY)is due to experimental errors,and the results of the cross-coupled inertia coefficients are not shown in the reference.In addition,the results calculated by the perturbed bulk flow model4with calibrated parameters are included in the comparison.

Compared with the analytic results,the direct stiffness of the calibrated bulk flow model has better agreement with the test data.Other results excluding the direct damping and cross-coupled damping agree well with the experimental results.Although the overestimation of the direct damping is apparent and the cross-coupled damping is underestimated,the comparisons indicate that CFD calibration improves the accuracy of bulk flow model.In the reference,the pressure loss coefficient is empirically selected to be 0.1,and the preswirl ratio and the friction factor are not shown explicitly.However,the calibrated pressure loss coefficient shown in Table 2 is much higher than 0.1.The reason may be that the direct stiffness and direct damping of the calibrated model are much higher than the results of the reference.

On the other hand,the transient bulk flow model provides further improvement for the rotordynamic coefficients in comparison with the perturbed bulk flow model.Although the overestimation of the direct damping still exists,the direct damping of the transient model is smaller than the result of the perturbed model and is closer to the experiment data.Meanwhile,the cross-coupled damping of the transient model is higher than the result of the perturbed model and has better agreement with the test data.The direct stiffness,crosscoupled stiffness and direct inertia coefficients of the transient model show slight improvement.The cross-coupled inertia coefficient of the present model is much larger than the result of the perturbed model.

4.Conclusions

(1)In comparison with other turbulence models,the k-ε standard turbulence model with the enhanced wall treatment is accurate to calculate the wall-bounded flow in seal clearance under low Reynolds number.

(2)The pressure loss coefficient decreases with increasing pressure drop;the preswirl ratio increases with increasing the rotational speed and decreases with increasing the pressure drop.

(3)The reaction force on the rotor surface has the same frequency with the rotor whirl motion.The radial force decreases and the tangential force increases with increasing the whirl frequency.

(4)The accuracy of the perturbed bulk flow model with the CFD calibration is improved in comparison with the analytical results in the reference.Furthermore,compared with the perturbed model,the predictions of the transient bulk flow model are in better agreement with experimental data.

Acknowledgement

This study was supported by the National Natural Science Foundation of China(No.11176010).

Appendix A

These are the sub-coefficients in the sub-matrix of the Eq.(9).The coefficient of lower-case letter is corresponding to the submatrix of capital letter.The subscript m and n is corresponding to the dimension of the sub-matrix.

References

1.Childs DW.Turbomachinery rotordynamics.Hoboken:John Wiley&Sons;1993.p.8–10.

2.Hirs GG.A bulk- flow theory for turbulence in lubricant films.J Lubr Technol 1973;95(2):137–45.

3.Childs DW.Finite-length solutions for rotordynamic coefficients of turbulent annular seals.J Lubr Technol 1983;105(1):437–45.

4.San Andres L.Analysis of variable fluid properties,turbulent annular seals.J Tribol 1991;113(4):694–702.

5.Ha T,Childs DW.Annular honeycomb-stator turbulent gas seal analysis using new friction-factor model based on flat plate tests.J Tribol 1994;116(2):352–60.

6.Childs DW,Kheireddin B,Phillips S,Asirvatham TD.Friction factor behavior from flat-plate tests of smooth and hole-pattern roughened surfaces with supply pressures up to 84 bars.J Eng Gas Turbines Power 2011;133(9):092504.

7.Marquette OR,Childs DW,San Andres L.Eccentricity effects on the rotordynamic coefficients of plain annular seals:Theory versus experiment.J Tribol 1997;119(3):443–7.

8.Lindsey WT,Childs DW.The effects of converging and diverging axial taper on the rotordynamic coefficients of liquid annular pressure seals:theory versus experiment.J Vibr Acoust 2000;122(2):126–31.

9.Saba D,Forte P,Vannini G.Review and upgrade of a bulk flow model for the analysis of honeycomb gas seals based on new high pressure experimental data.J Mech Eng 2014;60(5):321–30.

10.Arghir M,Nguyen MH,Tonon D,Dehouve J.Analytic modeling of floating ring annular seals.J Eng Gas Turbines Power 2012;134(5):052507.

11.Mariot A,Arghir M,Helies P,Dehouve J.Experimental analysis of floating ring annular seals and comparisons with theoretical predictions.J Eng Gas Turbines Power 2016;138(4):042503.

12.Moore JJ,Palazzolo AB.CFD comparison to 3D laser anemometer and rotordynamic force measurements for grooved liquid annular seals.J Tribol 1999;121(2):306–14.

13.Moore JJ.Three-dimensional CFD rotordynamic analysis of gas labyrinth seals.J Vibr Acoust 2003;125(4):427–33.

14.Untaroiu A,Hayrapetian V,Untaroiu CD,Wood HG,Schiavello J,Mcguire J.On the dynamic properties of pump liquid seals.J Fluids Eng 2013;135(5):051104.

15.Chochua G,Soulas TA.Numerical modeling of rotordynamic coefficients for deliberately roughened stator gas annular seals.J Tribol 2007;129(2):424–9.

16.Yan X,Li J,Feng Z.Investigations on the rotordynamic characteristics of a hole-pattern seal using transient CFD and periodic circular orbit model.J Vibr Acoust 2011;133(4):041007.

17.Arghir M,Billy F,Pineau G,Frene J,Texier A.Theoretical analysis of textured damper annular seals.J Tribol 2007;129(3):669–78.

18.Migliorini PJ,Untaroiu A,Wood HG,Allaire PE.A computational fluid dynamics/bulk- flow hybrid method for determining rotordynamic coefficients of annular gas seals.J Tribol 2012;134(2):022202.

19.Xu CW.Finite element method.Beijing:Tsinghua University Press;2003.p.130–56[chinese].

20.Warren Y,Richard B.Roark’s formulas for stress and strain.7th ed.New York:McGraw-Hill;2002.p.683.