A Small Perturbation CFD Method for Calculation of Seal Rotordynamic Coefficients

Seal rotordynamic coefficients link the fluid reaction forces to the rotor motion, and hence are needed in the stability calculations for the overall rotating systems. Presented in this paper is a numerical method for calculations of rotordynamic coefficients of turbomachinery seals with rotors nominally at centered, eccentric and/or misaligned position. The rotor of the seal is assumed to undergo a prescribed small whirling motion about its nominal position. The resulting flow variable perturbations are expressed as Fourier functions in time. The N-S equations are used to generate the governing equations for the perturbation variables. Use of complex variables for the perturbations renders the problem quasi-steady. The fluid reaction forces are integrated on the rotor surface to obtain the fluid reaction forces at several different whirl frequencies. The rotordynamic coefficients are calculated using appropriate curve fitting. Details of the model are presented, and sample results for concentric and eccentric annular incompressible flow seals are included to demonstrate the capability and accuracy of the proposed method.


INTRODUCTION
Seal Rotordynamics urbomachinery seals interface between rotating parts such as rotors, blade tips, and stationary parts such as housings.The seals, used to isolate regions of different pressures, are non-contacting, and allow a leakage flow across.The fluid flow that exists in seals generates reaction forces on the rotor.As a rotor moves away from its nominal operating position, the fluid flow in the seal is altered and reaction forces are generated on the rotor.Knowledge of the type and magnitude of these forces is important when calculating the stability char- acteristics of the overall rotating system.The reaction forces can be destabilizing e.g. in a labyrinth seal, Alford  1965], or with proper designs, the reaction forces can be used to stabilize the rotor and provide load bearing capabilities in a seal, Von Pragenau [1982, 1987].
As a seal rotor moves away from its nominal position, the reaction forces generated by the fluid flow can be linked to the rotor displacement, velocity and accelera- tions using a spring-damper-mass system in 2 directions.
For the seal configuration shown in Figure 1, the reaction forces Fy and F can be linked to the rotor motion using the relation IFl + -Mzy Mzz (1) where y,/9, Y and z, , % are the displacements, velocities, and accelerations in the Y and Z directions.Kyy and Kzz are the direct and Ky and Kzy the cross-coupled stiffness coefficients, Cyy, Czz and Cyz, Czy direct and cross- coupled damping coefficients and Myy, Mz and My, My the direct and cross-coupled inertia (added mass) coeffi- cients.Linked in this fashion, the rotordynamic coeffi- cients can directly be used in the rotor stability calcula- tions.When the nominal position of the rotor is concentric, the coefficient matrices become simpler and assume a skew-symmetric fo, given as _Fy= Kk y+ -k K z -c C + (2) -mM with a total of only six distinct coefficients as against twelve for the general case.

Rotordynamic Coefficient Calculation Methods
The influence of annular seal forces on the dynamics of centrifugal pumps was first investigated by Black [1969]  and Black and Jensen [1970].A bulk flow model was developed to compute the rotordynamic coefficients of annular seals.The theory treats the seal as a single lumped flow domain with no variations across the volume with friction factor correlations used to calculate wall shear.A perturbation in the rotor location is then used to generate the fluid reaction forces and subse- quently rotordynamic coefficients.This method can be used to calculate only the skew-symmetric coefficient set.Bulk flow models based on Hirs' lubrication equation [1973] were developed by Childs for short [1983]  and finite length [1983] annular seals.The bulk flow model was also extended to compressible flow seals by Nelson [1985] for annular and tapered seals.
Recently, more detailed solution procedures, based on 2-D solution procedures have been reported.The models account for flow variations in the circumferential as well as the axial direction but the variations across the film are integrated to give mean flow parameters.San Andres [1991] used the Hirs' lubrication equations with a 2-D solution procedure similar to the SIMPLE method for flow solutions in eccentric and misaligned seals and bearings under constant and variable properties.Simon  and Frene [1989, 1991] used N-S equations integrated across the film thickness to treat eccentric annular seals.
In both the methods, the flow equations are perturbed to yield equations for the static flow and the perturbations in the flow.With a prescribed, small motion of the rotor, the perturbation variables are calculated and the pressures integrated on the rotor surface to generate fluid reaction forces and the rotordynamic coefficients.
The above described methods use integrated, mean values of the flow variables, averaged over either the overall flow domain, or across the fluid film in the seal.
These models, then, can be expected to be adequate when the film thicknesses are small, and the flow is well behaved.When the seal has flow regions where there are steps, and/or regions with large film thicknesses e.g grooved or labyrinth seals, and/or when strong secondary motion is present e.g. in labyrinth seals, or as shown by Tam et.al. 1988], in eccentric annular seals, the average property models can be expected to deteriorate in accu- racy.For such cases a full, 2-or 3-D, detailed analysis the flow becomes necessary for accurate flow and rotor- dynamic predictions.
With the advent of powerful computers, and sophisti- cated numerical techniques, the computational fluid dynamics (CFD) analyses of seal flows have become economical.Such a study was conducted by Tam et.al.   [1988].A 3-D Navier-Stokes analysis was done on an eccentric, whirling annular seal with incompressible flow.The fluid reaction forces were used to compute the rotordynamic characteristics of the seal.In addition, flow fields at several seal configurations and boundary condi- tions were computed.Computed flow fields showed the existence of recirculation bubbles in the seal flows under certain flow conditions.Nordmann et.al. [1987, 1989]  have reported a finite-difference method based on the N-S equations for seal rotordynamics, where they have used a small-perturbation analysis on the N-S equations in cylindrical frame.The method was used to compute rotordynamics of nominally concentric, annular and grooved seals with incompressible and compressible flows.They also reported a study for eccentric, incom- pressible flow seals, Nordmann and Dietzen [1988].More recently, Baskharone and Hensel [1991] have reported a 3-D finite element solution method, coupled with a perturbation analysis to compute flows and rotordynamics of nominally centered seals.They have also used this method to calculate the fluid reaction forces on rotating parts such as impellers.
Solution procedures based of the full CFD solutions for seal flow and rotordynamics have been reported by Athavale et.al. [1992].One of the methods is based on a whirling rotor with a circular whirl orbit.A grid transformation is done to render the problem quasi- steady, and solutions at several whirl frequencies are calculated.The rotor surface pressures are integrated to yield fluid reaction forces and appropriate curves are fit to the force-whirl frequency plot to obtain the rotordy- namic.coefficientsThis method is similar to that used by Tam et.al. [1988].This method is fast, but is applicable only for calculation of the skew-symmetric coefficients of a nominally concentric seal.
For nominally eccentric and/or misaligned seals, the full CFD solution method is based on the "shaker" method used in experimental determination of the rotor- dynamic coefficients (e.g.Childs et al [1996]).The rotor is moved along a radial direction from a nominally concentric or eccentric position, and the motion is sinusoidal in time.The resulting time-dependent flow is calculated using a moving grid algorithm.The pressures on the rotor surface are then integrated to provide time-accurate fluid reaction forces, which then are used to calculate the rotordynamic coefficients.This, method, based on the full CFD solutions, is flexible and can be used for all nominal rotor positions, however, the com- putational times involved can be expected to be higher than the whirling rotor method, due to the time-accurate solutions that have to be obtained.
The method presented in this paper is developed to treat seals with rotors that are nominally eccentric and/or misaligned.To avoid the high computational times asso- ciated with the full, time accurate CFD solutions, the method is based on a small-perturbation analysis of the 3-D flow field in the eccentric seals.The methodology is similar to that described by Nordmann and Dietzen  [1988], but the present approach uses the generalized Body Fitted Coordinate (BFC) formulation of the flow equations.For this reason, the present formulation is expected to be applicable to cylindrical seal problems as well as a variety of other problems with complicated shapes, e.g.blade tip flows and impeller flows.
Y Yo + ero cos 12t z Zo + ero sin Ot (3)   where Yo and Zo correspond to the nominal center location, e is a small number, and r o is a suitable length scale, e.g.nominal seal clearance, Co.The resultant time-varying flow variables are now expressed as, e.g.tt U 0 + eft iv v 0-+-ev w--w 0 q-ew p Po + ePl (4) where the subscript 0 corresponds to the time-indepen- dent values that describe the steady flow in the seal with the rotor at its nominal position, and the time-dependent perturbations are denoted by the subscript 1.Since the rotor center motion has sine and cosine functions of time, we assume that the time-dependent part of the resultant perturbations can also be expressed using Fourier series in time.Thus, we have, e.g.,

Ig
Ulc COS Ot + Uls sinlt V1 Vlc COS lit + Vls sinl2t (5 where the values Ulc Uls etc. now are functions of space only.When solutions of turbulent flows in the seals are considered, the variables connected with turbulence: turbulent kinetic energy (k), dissipation rate (), and the turbulent viscosity (t) are assumed to be unaffected by the rotor center motion, and hence are held constant at their steady-state, 0th order values.
To derive the equations for perturbation variables, we start with the N-S equations, and introduce a coordinate transformation to map the time-varying grid locations (x,y,z,t) to a time-independent grid frame (I, -q, , 'r).The resulting equations are: Continuity:

DESCRIPTION OF THE SOLUTION METHODOLOGY Governing Equations
The seal flow is described using the 3-D Navier-Stokes (N-S) equations in the generalized body-fitted-coordinate (BFC) system.The perturbations in the flow are intro- duced by assuming that the rotor center undergoes a whirl about its nominal, steady-state position.The rotor spin and whirl velocities are assumed to be o and 1) respectively (see Figure 1).The rotor center position, in general, can be described as (6) Momentum: (7) where + can be any of the Cartesian velocities u, v, or w, k are the coordinate directions, J is the Jacobian of transformation, e k are the contravariant base vectors, the elements of metric tensor and F is the diffusion coefficient.V is the grid velocity due to the rotor motion and is of order e.The term S,I contains pressure gradient terms and additional body force terms.
The definition of the dependent variables, Eqns. 4,are now used in the above equations.All terms containing powers of e higher than are ignored, and then the terms with and without the factor e are separated to yield two sets of equations.The 0th order equations now have no time-dependent terms or terms with grid velocity, and describe the steady-state seal flow.The 1st order equa- tions still have time-dependent terms.To render them quasi-steady, the Fourier expansions of the 1st order variables, Eqns. 5,are used in the 1st order equations, and the terms containing sine and cosine functions of time are separated out to yield independent flow equations for each of the Ulc, Uls etc.Since these variables are functions of space only, a steady-state solution treatment can be used to solve the equations for these 1st order variables.To avoid computations of, e.g.Ulc and Uls independently, we combine them to yield new variables: 1 Ulc + iUls "1 Vlc -]-iVls etc (8) where -1 and/1, now are complex variables.The flow equations for the complex variables are obtained in a similar manner, i.e.
(Eqn for/1) (EqFt.for Ulc + i(Eqn for Uls etc. (9) The flow equations for the 0th order, as stated before, are the steady-state N-S equations and are not repeated here.
The equations for the 1st order variables and for incom- pressible flows are: Continuity: where Vg is the complex grid velocity and V1 is the complex velocity vector.Momentum: grid motion terms contain complex grid metric coeffi- cients, and the 0th order flow variables and hence are independent of the 1st order flow variables.

Solution Methodology
The 0th and 1st order equations have a very similar structure, and can be solved using similar methodologies.
The starting point for the solutions of the 0th order was an advanced code, SCISEAL, which has a finite volume integration scheme, uses Cartesian velocities are velocity variables in a colocated arrangement, and a pressure based-solution algorithm.The flow equations are solved sequentially, and the velocity-pressure coupling is provided using a modified SIMPLEC method.The code has capabilities of high spatial (up to 3rd order) and temporal (2nd order) accuracy, and has a variety of turbulence models (including the standard k-G, low-Re k-G, multiple scale k-e and the 2-layer model) and treatment for isotropic surface roughness.
The 1st order flow module was added to SCISEAL to provide a seamless solution procedure for rotordynamics, which involved the following steps: 1. Solution of the steady-state equation for nominal rotor position. 2. Solution of 1st order equations at 5 different whirl frequencies.
3. A least-square curve fitting algorithm for rotordynamic coefficients.The 1st order equations are solved in a manner very similar to the 0th order equations.A sequential procedure is used to solve the equations, and pressure-velocity coupling is provided by the modified SIMPLEC algorithm.The major differences between the 1st and the 0th order equations are: a) The 1st order equations are constant coefficient equations, and hence linear in nature.This affords faster convergence of the 1st order set, and b) The 1st order variables are complex quantities and separate modules are required to treat the equations.The complex arithmetic also tends to increase the computa- tional times.
where the source term now contains the 1st order pressure gradient source terms, e.g. e OlOfo as well as the source terms that arise due to the grid motion.The

Boundary Conditions
Boundary condition specifications for the 1st order equations is straightforward.A Neumann boundary con- dition on any 0th order flow variable translates to a Neumann condition on the 1 st order variable.A Dirichlet condition on a 0th order flow variable does not allow a change in that quantity, so that the corresponding 1st order variable value becomes zero on that boundary.A seal-specific inlet boundary condition where a stagnation pressure and entrance loss factor is specified translates as: where u is the Cartesian velocity along the axis of the seal, Po is the specified stagnation pressure, is the loss coefficient, and subscript b refers to the value to be applied at the inlet boundary.

Calculation of Rotordynamic Coefficients
To calculate the rotordynamic coefficients, the fluid reaction forces on the rotor wall are needed.These are obtained by integrating the 1st order pressure solutions on the rotor wall, and then the total force is resolved along the y and z directions to yield the complex reaction forces Fy and Fzl.Using the definition of the rotor center location, Equation 3, expressions for the rotor center velocity and accelerations are found, and then substituted in the definition of the rotordynamic coefficients, Eqn. 4.

Czy . Mz -2 (16)
To evaluate the individual coefficients, the variation of reaction forces with whirl frequency is needed.This is done by calculating the 1st order solutions at several different whirl frequencies.Appropriate curves are then fit (2nd order polynomials) to the force Vs. frequency relations to obtain the complete set of rotordynamic coefficients.

NUMERICAL RESULTS
Two test cases are presented here to demonstrate the capability and accuracy of the perturbation method.In both cases experimental and/or published numerical data is used to validate the method.

Long Incompressible Flow Annular Seal
Experimental data for a long annular seal with incom- pressible flow was reported by subo [1992].The seal dimensions were: length 240 mm., rotor radius 39.656 mm., clearance 0.394 mm., entrance loss factor 0.5.The stator is nominally concentric, and results in a skew-symmetric set of coefficients.The boundary conditions were (Figure 2): specified static pressure at seal exit, specified pressure and loss coefficient at inlet.Experimental measurements of inlet swirl were imposed at the inlet boundary.
Solutions were obtained at four different rotor spin speeds: 600, 1080, 1800 and 3000 rpm, and at seal pressure differentials (Ap) of 20, 50, 100, 250, 500 and 900 KPa.The low-Re k-e model was used and the computational grids used had 15 cells in the axial direction, 30 in the circumferential direction and 12 or 16 cells in the radial direction, depending on the pressure differential applied.Central differencing with 0.1 damPing was used for spatial discretization of convective terms.The results of the calculations are shown in Figures 3 through 7. The experimental data are also plotted on these plots, and very good correlation between the two sets is obtained.The discrepancy in K value at 1980 rpm and 900KPa pressure is probably due to the Cross coupled damping coefficients, long annular seal low-Re kturbulence model.The computed stiffness coefficients K, k, inertia coefficients M, m, and cross- coupled damping coefficients C, c, were found to be within 10% to 15% of the experimental values for a majority of the data points.The direct damping coeffi- cients, however showed somewhat higher discrepancies (up to about 40%) in the lower Ap values and at higher rpm.This long seal exhibits negative direct stiffness coefficients at high rotor rpm over the entire Ap range, as well as at low pressure differentials at the lower rpm.The numerical results predict the negative stiffness as well as the crossover to positive values at higher seal pressure differentials Ap and low rpms.
The analysis given by Tam et.al. 1988]was applied to a few data points from the present computations to check the ability of the perturbation model to predict flow physics as outlined in their paper.Seal pressure differ- entials of 50 and 500 KPa at a rotor speed of 1080 rpm were selected as the two sample points.In the present set of calculations, the real parts of the reaction forces F y and Fz correspond to the direct and quadrature dynamic stiffness as described in Tam et.al.[1988].
Plot of the quadrature dynamic stiffness, K Vs. the whirl frequency [1 for the two pressure differentials is shown in Figure 8.The curves are straight lines that indicate that the rotordynamic coefficients are indepen- 1.0E35 AP I: 1.0,08 dent of the whirl frequency in these runs.Values of the average circumferential velocity ratio, h, for the 50 and 500 KPa cases are 0.478 and 0.412 respectively.This corresponds to a spin-up of the fluid in the seal from specified inlet h values of 0.46 and 0.32 (from experi- mental data).This inlet swirl generates positive values for the so called cross-coupled stiffness coefficient k.
The direct dynamic stiffness is a function of the modified whirl speed (f to), and Figure 9 shows the plots.Slopes of the lines give the direct inertia coefficient M. In the model described by Tam et.al.[1988], the direct stiffness coefficient is defined as the intercept on the vertical axis in Fig. 9; and is positive for both Ap values.However, as defined in Eqn. 2 in the present calculations, the direct stiffness coefficient K is related to the intercept Kr,o as K Kr,o M h 2 to2. (17) The negative K for Ap 50KPa in the present calcula- tions results due to the inertia term, and a small intercept Kr,o, while the larger intercept at 500KPa results in a net positive K.
Annular Eccentric Seal This case involves an annular seal with incompressible flow, with nominal static eccentricity ratios varying from 0 to 0.7.The seal dimensions were: seal length 40 mm., seal radius 80 mm., clearance 0.36 mm., entrance loss factor 0.5.The boundary conditions were similar to those before.Experimental data for the stiff- ness coefficients was reported by Falco et.al 1986], and a compilation of various numerical data by Simon and  Frene [1991].Solutions were obtained at a rotor speed of 4000 rpm, average preswirl of 0.3, a pressure differential of 1.0 Mpa, and at several static eccentricity ratios from 0.0 to 0.7.The standard k-e model with wall functions was used and the grid used had 12 cells in the axial, 6 in radial and 30 in the circumferential direction.Central differencing with 0.1 damping was used for convective terms.Results of these calculations are shown in Figures 10 through 19.The experimental as well asnumerical data is also plotted for comparison, and very good 1.EOS-Fy lq,-K M(Q-X= --to ",,, (C) .0.0"1'% "",, K. Nrn re= 0 4,Sx 105 , =_1' " , ., ., , .correlation between the present results and published experimental and numerical data is seen.

CONCLUDING REMARKS
A small-perturbation method for rotordynamic coeffi- cients of seals has been developed.The method is based on the perturbations in the 3-D, Navier Stokes equations for BFC grids.The method has the capability of treating nominally eccentric and/or misaligned rotors.Solutions presented here included incompressible flow concentric and eccentric annular seals.In the concentric seal, numerical results correlate well with experimental data, including the negative direct stiffness that result in this long annular seal.For the eccentric seal case, the present results were compared with experimental and other numerical data, and again good comparison between thetwo was obtained for all coefficients except one cross-coupled damping coefficient.The use of 3-D BFC grids and N-S equations allow this model to treat a wide range of seals, including seals such as labyrinth, grooved, tip and stepped.This method can also be easily extended to compute interactions between fluid flow and larger rotating components, e.g.shrouded impellers.Current plans include extension to compressible flows and more complicated geometry seals.direct damping coefficients, N-s/m cross coupled damping coefficients; N-s/m small number, perturbation parameter Force, N elements of metric tensor complex axis Jacobian of transformation direct stiffness coefficients, N/m cross-coupled stiffness coefficients, N/m direct inertia (mass) coefficients, Kg cross-coupled inertia coefficients, Kg pressure, N/m length scale for whirl radius, m source term in momentum equation, N time, Cartesian velocities in x, y and z directions, fla/id velocity vector, m/s complex grid velocity

FIGURE
FIGURE