Review and Upgrade of a Bulk Flow Model for the Analysis of Honeycomb Gas Seals Based on New High Pressure Experimental Data

Since the 1960s, smooth-rotor/honeycomb-stator annular seals have been used in process centrifugal compressors where they were employed to replace aluminium labyrinth seals consumed by the process fluid. As it turned out, honeycomb seals had significantly less leakage compared to conventional see-through labyrinth seals for the same clearances. In a smooth-rotor/honeycomb-stator annular seal, a honeycomb pattern of hexagonal cavities is present on the stator, see Fig. 1. Hole-pattern seals are very similar. The same pattern of cavities is present, but the their shape is cylindrical, see Fig. 2. It has long been known that internal seals can have a strong impact on the dynamics of a turbomachine. A destabilizing effect due to the seals was reported for the Kaybob compressor in 1975 [1] and for the Ekofisk compressor in 1976 [2].


INTRODUCTION
Since the 1960s, smooth-rotor/honeycomb-stator annular seals have been used in process centrifugal compressors where they were employed to replace aluminium labyrinth seals consumed by the process fluid.As it turned out, honeycomb seals had significantly less leakage compared to conventional see-through labyrinth seals for the same clearances.
In a smooth-rotor/honeycomb-stator annular seal, a honeycomb pattern of hexagonal cavities is present on the stator, see Fig. 1.Hole-pattern seals are very similar.The same pattern of cavities is present, but the their shape is cylindrical, see Fig. 2.
It has long been known that internal seals can have a strong impact on the dynamics of a turbomachine.A destabilizing effect due to the seals was reported for the Kaybob compressor in 1975 [1] and for the Ekofisk compressor in 1976 [2].The damping capabilities of honeycomb seals were noticed and exploited later.In 1985, honeycomb seals were used to stabilize the high-pressure oxygen turbopump of the space shuttle main engine [3].Currently, honeycomb seals are often used for the balance piston seal of a high performance compressor, especially in a back-to-back arrangement.When a seal is positioned half way between the bearings, where the first mode of vibration of the shaft has its maximum amplitude, its dynamic behavior becomes of greater importance.
The industry would benefit from a reliable way of predicting the behavior of honeycomb and holepattern gas seals.However, despite the combined analytical and experimental effort, the predictive ability of the presently available tools is not sufficient, especially at very high pressures.
In this paper, new experimental data for honeycomb gas seals are presented and compared to

Review and Upgrade of a Bulk Flow Model for the Analysis of Honeycomb Gas Seals Based on New High Pressure Experimental Data
Saba, D. -Forte, P. -Vannini, G. Diego Saba1,* -Paola Forte 1 -Giuseppe Vannini 2the solutions of a simplified fluid dynamic model.In the current program the seals have been tested at higher pressures than in previous works, to the authors' knowledge.The fluid dynamic model is based on an isothermal bulk flow model developed by Kleynhans and Childs [4] and Kleynhans [5].Since the simulations with the isothermal model did not agree with the authors' experimental data, a new simulation tool has been developed which adds the possibility to make different assumptions on the transport of heat and momentum.
Section 1 introduces the dynamical coefficients, as a way to characterize the dynamical behavior of seals.Section 2 briefly describes the testing apparatus used to measure the dynamic coefficients.Section 3 gives a detailed description of the analytical model used in this paper.Since the aim of this work was the comparison of different variants of the bulk-flow formulation, the model is described as the composition of different submodels, each one governed by its own parameters and options.Sections 4 and 5 present the results and give some concluding remarks.

DYNAMIC COEFFICIENTS
The dynamical behavior of a rotordynamic component is characterized by means of dynamic coefficients that describe the linear relationship between the displacements of the rotor and the forces it exerts on the surrounding bodies, in this case the gas in the clearance.The relevant displacements for gas seals are the lateral displacements of the rotor.The relationship can be expressed as: where x  , y  are the Fourier transforms of the lateral displacements in two fixed orthogonal directions and F x  , F y  are the Fourier transforms of the corresponding forces.The frequency-response function H ij is a dynamic stiffness, but is called impedance in this paper, following the use of other authors in the same context [4], [6] and [7].This choice has some advantage because the term stiffness can be used for the real part, i.e. the in-phase response, without further specification.Frequency dependent stiffness and damping coefficients are defined as: ( The terms on the diagonal H xx , H yy are called direct impedances, while the terms off the diagonal H xy , H yx are called cross-coupled impedances. Forces and displacements in Eq. ( 1) are represented through their Cartesian components.In other words, the displacements are decomposed into harmonic oscillations in two fixed directions.
It is also possible to represent the lateral displacements of the rotor as the superposition of forward and backward circular orbits, using a new basis.The transformation rules for the change of basis are: where x f  , xb  are the components in the new basis.The forces are decomposed accordingly into forward and backward rotating forces.The transformation rules are the same: In the new basis, Eq. ( 1) becomes: and defines forward and backward impedances.The dynamic behavior is called isotropic if it is independent of the orientation of the reference system x, y.For moderate eccentricities, the behavior of a seal is usually close to isotropic.It can be shown that, for an isotropic system, the following properties hold: The effective stiffness and damping are defined as: and are important indicators of the performance of a seal.When the rotor axis moves along a forward precessional orbit, a positive effective stiffness indicates a centripetal reaction force, and a positive effective damping indicates a tangential stabilizing force, see Fig. 3.For an isotropic system the following relations hold: Since forces and displacements are real variables, their Fourier transforms have Hermitian symmetry.The four impedance coefficients, in x / y notation, also have Hermitian symmetry, namely: On the contrary, in f / b notation the following relations hold:

EXPERIMENTAL DATA
The experimental results available for this work come from the ultra high pressure (UHP) test rig of GE Oil & Gas, Florence.The test seal was a smooth-rotor honeycomb-stator seal, with convergent clearance, diameter 220 mm and length 65 mm.The ranges of values for the test parameters are given in Table 1.

Testing Apparatus
Only the most relevant information will be given here.
A more detailed description of the testing apparatus can be found in [8].Two identical seals are tested at once.They are mounted in the test cell in a symmetrical back-to-back configuration, see Fig. 4. The relevant conditions for a compressor seal are reproduced, namely the high and low pressures, respectively upstream and downstream the seals, the rotational speed of the rotor, and the swirl of the gas before entering the seals (preswirl).
The gas used in the plant is nitrogen.The design pressure is 400 bar.
The rotor is supported by active magnetic bearings (AMB).The magnetic bearings can be controlled in such a way as to impose the desired displacements to the axis of the rotor.The relative position between the shaft and each of the bearings is measured, and so is the force actuated by each bearing.The orbit imposed to the rotor axis is given by the superposition of harmonic displacements at different concurrent frequencies (multi-frequency excitation).Forces and displacements are measured in two configurations: with a pressure difference across the seals, and without pressure.The reason for this is the need to separate the dynamic coefficients of the seals from those of the whole system.
The main contribution to the dynamic coefficients of the system, besides that of the seals, is due to the inertia of the rotor.This contribution can be computed through the displacements and the known mass distribution of the rotor.However, a direct measurement is preferred.
The preswirl ratio, i.e. the ratio of the circumferential velocity component of the gas and the peripheral speed of the rotor, is generated through the swirler ring which collects the gas flow from an inlet plenum and injects it toward the seals through aerodynamic nozzles evenly distributed along the circumference.The current swirler is designed to produce a high preswirl, in the range of 0.8 to 1 at 10 krpm.

Identification Methodology
The measurements of forces and displacements are sampled at a sufficiently high rate to obtain a reliable spectral decomposition in the range of the excitation frequencies.
For each spectral component, Eq. ( 1) holds.It consists of two scalar equations, that are not enough to determine the four impedances of the general formulation.However, assuming isotropy, there are only two unknowns The unknown impedances can be put in evidence by rearranging the terms: The shape of the orbit of each selected frequency should not be circular or near to circular, because the matrix of coefficients in Eq. ( 12) would become illconditioned.Flat straight-line orbits are preferred with this identification method.
To obtain the four impedances of the general case, two experiments are necessary, differing only for the excitation components and not for pressures, temperature and preswirl.Strict acceptability criteria must be met to couple two experiments.For the two coupled experiments, we can write: which is immediate to solve.

Bulk Flow Model
The bulk flow model with two control volumes, devised by Kleynhans and Childs [4] and Kleynhans [5] is an evolution of the bulk flow model used by Nelson [9] for smooth annular seals.A second control volume was added to take into account the gas trapped (actually recirculating) in the cavities of the honeycomb.The effective speed of sound is slowed down by the exchange of mass between bulk flow and cavities.
The bulk flow model reduces the dimensionality of the problem, from three to two.The balance equations are written in integral form along the radial direction, and the mean values of the properties of the flow are assumed to be close to the values in the bulk.The bulk flow variables are thus functions of only two geometrical coordinates, axial and circumferential.
In Fig. 5 the two control volumes are shown.Volume A is the clearance and volume B comprises the cavities of the honeycomb, or the holes of the hole pattern.
In the balance equations (Eq.( 14)) the symbols have the meaning shown in Section 6 and commas indicate partial derivatives with respect to the variables on their right.The subscript A has been omitted from the fluid properties of the corresponding control volume.The unknowns are ρ, T, u θ , u z , ρ B , T B .Any other thermodynamic quantity is derived from these using the equations of state.In this work, the fluid was modeled as a perfect gas with constant heat capacity, viscosity and thermal conductivity.
Radial fluxes φ m etc. between control volumes, friction stresses τ Rθ etc. and heat fluxes q R etc. must be determined through additional empirical equations.Different models result from different assumptions.
Two assumptions are shared by all the variants of the model.The pressure is assumed to be the same for the two control volumes at a given location.The temperature too is assumed to be the same for the two volumes.The assumption of thermal local equilibrium is not well supported, either theoretically or experimentally, but a more complicated model did not seem justified at this stage.Future experimental or theoretical feedback might suggest a better hypothesis.
The equations are first solved to find the steady state, independent of t, θ.Then, a first-order perturbation problem is set up.A precessional motion is imposed to the rotor, that entails a harmonic perturbation of the clearance, in time and along the circumferential coordinate.
To complete the problem, boundary conditions must be set for the fluid at inlet and outlet.Four conditions are needed: two on the pressure, at inlet and outlet, one on the inlet swirl velocity, and one on the inlet temperature.
The perturbation problem silently assumes that the flow can be described for every position z, θ and for every instant t by the same state variables as for the steady problem, namely ρ, T, u θ , u z , a.This assumpion, however, is much stronger than for the steady problem.It implies that, at the frequencies of interest, the perturbations are so slow that the vorticous flow in a honeycomb cavity differs little from the steady one identified by the same state variables.

Variants of the Bulk Flow Model
The variants of the bulk flow model implemented for this work are shown in Fig. 6.The options for the three submodels can be chosen independently.
submodel options heat exchange isothermal non-isothermal adiab.n q = 0 n q > 0 isoth.wall n q → ∞ entrance and exit loss incompressible compressible radial convective momentum flux

Fig. 6. Variants of bulk flow model
The model described in [4] and [5] can be reproduced using the options in column two (isothermal, incompressible, f p = 1).This model has been chosen as a reference state-of-the-art model.All the other options implement original submodels, to the authors' knowledge.
The radial fluxes φ pθ , φ pz are written as: where a non-dimensional f p , chosen by the user, has been introduced to estimate the convective part of the momentum fluxes.By imposing p = p B , T = T B , and using Eq. ( 15), we can dispose of the unknowns p B , T B and write the bulk flow equations in the form: where The isothermal model does not use the energy equation.

Friction Model
The friction model is the same for all the variants.Friction is assumed to be isotropic: the friction coefficient is independent of flow direction and the shear stresses are parallel to the relative velocity between the bulk flow and each wall.The magnitude of the relative velocity is indicated in the following formulas by u r .The same friction model is used for the rotor and the stator: where Eq. (18) defines the Fanning friction factor and Eq. ( 19) is Blasius friction formula.Different parameters m, n are used for the stator and for the rotor.Values for the rotor can be found in the manuals.For the stator, there is no simple predictive model.Experimental values were obtained on a flat plate tester [10] for some geometries and Reynolds numbers.

Heat Transfer Model
The most important phenomenon that the heat equation and the heat transfer model should be able to describe is the amount of heat produced by friction and how much of it is transferred through the walls.The effect of the rate of heat transfer on the speed of sound is also important.The temperatures of rotor and stator are assumed to be known and equal to the temperature of the fluid upstream.
For the heat transfer, a simple correlation has been posited between Nusselt and Reynolds numbers Eq. ( 20).The heat transfer model is summarized by the formulas: where m q = 0.8 and the extreme cases for n q are 0 (adiabatic) and ∞ (isothermic wall).The adiabatic wall temperature T aw , as defined by Shapiro [11], takes into account the fact that, for a compressible fluid, the temperature in the bulk differs from the temperature in the viscous layer near the wall, even when the wall is adiabatic.The adiabatic wall temperature differs a little from the stagnation temperature, and their ratio R w Eq. ( 20) depends mainly on Prandtl number.Since Prandtl number little for air, nitrogen, and other fluids of interest in compressors, a fixed value R w = 0.9 was used.The correlation in Eq. ( 20) can be viewed as a simplification of Sieder-Tate correlation: where the dependence of viscosity, thermal conductivity and heat capacity on pressure and temperature have been neglected.
The "isothermal wall" option is the limit case for n q → ∞, and has constant adiabatic wall temperature T aw .The bulk temperature is not constant, but decreases with velocity.

Boundary Conditions
The fluid, upstream, before being dragged by the movement of the rotor, can already have a circumferential velocity component, or swirl.It has been shown that the circumferential velocity affects greatly the dynamic coefficients and the stability of honeycomb seals.The effective damping is decreased by a positive swirl, i.e. in the same direction of the shaft's rotation, and is increased by a negative swirl.
Two of the boundary conditions reflect the assumption that circumferential momentum and total enthalpy are conserved at the entrance: The condition on enthalpy does not apply for the isothermal model.
The conditions on pressure are expressed in terms of loss coefficients, assigned by the user.The loss coefficients are the ratio of the amount of energy transformed irreversibly into heat and the kinetic energy in the conduct.For a perfect gas, assuming an adiabatic transformation, the following boundary conditions result [9]: where ξ ξ α γ γ For the isothermal case, conditions in Eq. (23) become: A formulation for incompressible fluids can be chosen by the user, identical for isothermal and nonisothermal models, This formulation was used in [4] and serves here as reference.It can be viewed as an approximation for low velocities, as in that case compressibility effects become less important.In fact, the axial Mach number at the entrance never exceeded 0.23, even in choked conditions, and no significant differences have been found using compressible and incompressible boundary conditions.
Care must be taken to recognize a choked flow.The axial velocity u z increases with the axial coordinate, as in a Fanno flow, even if the clearance is slightly divergent, because the expansion due to friction forces prevails.The flow cannot reach the critical velocity at any point along the seal, except at the exit, where it is: When the choked condition is reached, the second boundary condition of Eq. ( 26) must be replaced by Eq. ( 27).In choked condition, the flow is insensitive to the downstream pressure p D .

RESULTS
The data from nine experimental tests were compared with the predictions of the model with the options and parameters of Table 2.The loss coefficients were set to k I = 0, k E = 1, meaning a reversible transformation at the entrance, and p D = p at the exit, i.e. no pressure recovery.The friction factor parameters for the stator m s , n s were not available from dedicated tests.An estimate was made using the measured mass flow rates.
The estimate was made by running the simulation with "isot" options and minimizing the quadratic percentage error of the mass flow rates.Table 3 shows the estimated n s as a function of arbitrarily chosen m s .As can be seen from the table, the estimate of n s fits quite well the experimental data for a wide range of m s values.The choice of m s has little influence on the dynamic coefficients.Experimental tests with different pressures but equal pressure ratio show a similar behavior, as expected.Fig. 7 shows experimental data in non-dimensional form for different pressures but similar pressure ratios p D / P U = 0.79÷0.83.Angular frequencies and impedances are non-dimensionalized respectively with ω ref and K ref defined by: It should be noted that the non-dimensionalized frequency ω / ω ref , scales with the inverse of the velocity of sound c -1 , and hence with T -1/2 .Different tests were run with the same excitation frequencies, but slightly different temperatures, since the ambient temperature was not controlled.The temperature shift of the non-dimensionalized excitation frequencies is evident in Fig. 7.It can also be noted that this method of non-dimensionalization is indeed effective in comparing different tests by similarity.
The impedances H bf , H fb , that express the coupling between forward and backward precessions, are expected to be zero in case of isotropic dynamic behavior.The experimental data confirm the isotropy, within the uncertainties of the measurements.
Simulation results are shown for two cases, one with unchoked flow and one with choked flow, see Table 4.
The Mach number in the table was computed with the isothermal model.The simulation results are displayed in Figs. 8 and 9.It can be noted that the graphs of Im(H ff ) do not pass through the origin, but are negative at zero frequency.This is the effect of the positive swirl, and the reason why swirl brakes have a stabilizing effect.
In both cases, the adiabatic model predicts a temperature rise up to 14 K in the bulk, due to friction.The comparison between isothermal and adiabatic models is worth commenting.In the adiabatic model, the temperature should increase along the flow, because the high friction near the walls converts a sizable amount of mechanical energy into heat.The analysis has shown that this effect is partially compensated by a temperature decrease due to expansion.In fact the net result was a variation of the absolute temperature of about 5% in the worst case.
Another difference is in the speed of sound, that in the adiabatic model is about 20% greater than in the isothermal model.Yet its variation seems to have little impact.
The differences among the models are more evident for the choked case.This is probably due to the greater compressibility effects, which increase the sensitivity to the thermal hypotheses.
The graph of the choked case shows an almost constant offset between measured and simulated Re(H ff ).This is an indication that the actual clearance, during the experiment, was possibly narrower than at rest.The elastic deformations of the seal and its housing are a possible cause.Negative static stiffness can lead to instability.This kind of instability has been observed in long seals working in choked condition.
In analyzing the Im(H ff ) graph, it is useful to separate two kinds of discrepancies: a constant offset and difference in slope.A constant error can be ascribed to an error in swirl prediction, which, in turn, can be justified by anisotropic friction of the honeycomb, not difficult to implement.On the contrary, in order to justify a different slope in the graph a deeper revision of the model is needed, probably related to the observation reported at the end of section 3.1.

CONCLUSIONS
To predict the dynamic behavior of honeycomb gas seals, an isothermal bulk flow model is currently used.A new simulation tool has been developed which adds the possibility to make different assumptions on the transport of heat and momentum.
New experimental data at higher pressures than previously available have been compared with the bulk flow simulations.The results showed that the new models have no major impact on the dynamic coefficients, in the frequency range of interest, except for the effective stiffness of choked flows.
The fact that heat generation and transport have little effect on the temperature of the fluid and on the dynamic coefficients of the seal is a remarkable nontrivial result.
Since neither the isothermal model nor the new tentative bulk-flow approaches have been satisfactory in explaining the available experimental data, it is likely that a key factor has been overlooked.The assumptions underlying the perturbation equations are probably too restrictive and need a revision.This does not necessarily mean that bulkflow approaches have to be abandoned.The ease of calculation is still a strong point in their favor and obliges to investigate unexplored model tuning parameters.

Fig. 7 .
Fig. 7. Impedances in non-dimensional form from four comparable experiments; the legend indicates entrance/exit pressures (bar)

Table 1 .
Test parameters

Table 2 .
Simulation options for comparison with experimental data

Table 3 .
Estimate of friction coefficients

Table 4 .
Cases shown in the graphics