# Design and Analysis of CFD Experiments for the Development of Bulk-Flow Model for Staggered Labyrinth Seal

Authors: Filippo Cangioli; Paolo Pennacchi; Leonardo Nettis; Lorenzo Ciuchicchi
Journal: International Journal of Rotating Machinery
 (2018)

---

## Abstract

Nowadays, bulk-flow models are the most time-efficient approaches to estimate the rotor dynamic coefficients of labyrinth seals. Dealing with the one-control volume bulk-flow model developed by Iwatsubo and improved by Childs, the “leakage correlation” allows the leakage mass-flow rate to be estimated, which directly affects the calculation of the rotor dynamic coefficients. This paper aims at filling the lack of the numerical modelling for staggered labyrinth seals: a one-control volume bulk-flow model has been developed and, furthermore, a new leakage correlation has been defined using CFD analysis. Design and analysis of computer experiments have been performed to investigate the leakage mass-flow rate, static pressure, circumferential velocity, and temperature distribution along the seal cavities. Four design factors have been chosen, which are the geometry, pressure drop, inlet preswirl, and rotor peripheral speed. Finally, dynamic forces, estimated by the bulk-flow model, are compared with experimental measurements available in the literature.


---

## Body

## 1. Introduction

Labyrinth seals and their sealing principles are commonplace in turbomachinery and they can arise in various configurations [1]. The most popular are the straight-through, the staggered, the slanted, and the stepped labyrinth seals. By their nature, labyrinth seals are noncontact seals that can harmfully rub against stationary parts [2]. Their working principle consists in reducing the leakage mass-flow rate by dissipating the flow kinetic energy via sequential cavities arranged to impose a tortuous path to the fluid (see the scheme of a staggered labyrinth seal in Figure 1). The speed and pressure at which they operate are bounded by their structural design.Figure 1
Schematic view of a staggered labyrinth seal.The clearance is defined by aerothermomechanical conditions that stave off the contact against the shroud under radial and axial excursions. In straight-through labyrinth seal, the angle at which the flow approaches to the teeth is usually 90°; the coefficients of discharge (Cf), peculiar of labyrinth seals, emphasize the effectiveness of the sealing element. The discharge coefficient is defined as the ratio between the actual mass-flow rate and that in ideal flow conditions. The relation between the sharpness of the tooth and the ability to restrict the flow has been given by Mahler [3] and more recently explored by means of computational fluid dynamics (CFD) in [4].Labyrinth seals are effective in reducing the flow but have a strong effect on the dynamic behaviour and often lead to dynamic instabilities. These problems have been addressed by several researchers starting with Thomas [5] and Alford [6]. They recognized that dynamic forces can lead to instabilities. Childs [7] addressed the root cause to the swirl velocity and introduced the swirl brake at the seal inlet to mitigate the circumferential velocity of the fluid.The current trends in the oil and gas market lead the manufacturers to maximize the turbomachinery efficiency. This target pushes the rotor dynamic design to be challenging. Therefore, the accuracy of rotor dynamic analyses has become of primary importance to avoid undesired instability phenomenon. Because labyrinth seals are one of the major sources of destabilizing effects, it is straightforward that the accurate prediction of labyrinth seal dynamic coefficients is crucial.The bulk-flow model is the most used approach for estimating the seal dynamic coefficients thanks to the low computational effort compared to CFD calculations. It was developed by Iwatsubo [8] and improved by Childs and Scharrer [9] in the 80s. The bulk-flow model is a control volume approach: the continuity and circumferential momentum equations are solved for each cavity, allowing the circumferential fluid dynamics within the seal to be investigated. For the axial fluid dynamics, empirical correlations are required, aiming at estimating the leakage and the pressure distribution. No turbulence model is employed in the bulk-flow model; therefore also the effect of the turbulence is demanded to correlations that estimate the wall shear stresses. The leakage strictly influences the swirl velocity profile along the seal and consequently the resulting dynamic coefficients.Dealing with staggered labyrinth seals, few models allow the dynamic coefficients to be estimated, despite the fact that this seal configuration is widely used in turbomachinery, especially in steam turbines. The advantage of this configuration, with respect to straight-through labyrinth seal, is the reduction of the leakage mass-flow rate given equal clearance.Wang et al. [10] use the straight-through bulk-flow model for the staggered labyrinth seal; Childs [7] considered the kinetic energy carry-over equal to the unity for all the cavities. Kwanka and Ortinger [11] compared the numerical results of the bulk-flow theory with experimental results performed on a staggered labyrinth seal. They revealed the strong connection between the circumferential flow and the axial one. The numerical results were not accurate compared to the experimental ones. The leakage correlation used by Kwanka et al. did not reproduce well the “zig-zag” behaviour of the discharge coefficients, because they considered an average value. Li et al. [12] investigated numerically and experimentally the effects of the pressure ratio and rotational speed on the leakage flow in a staggered labyrinth seal. In their paper, the authors do not report any correlation for the prediction of the leakage flow.This paper aims at filling the lack of analytical models for staggered labyrinth seals. Firstly, the authors performed a design and analysis of computer experiments (DACE) using steady-state CFD analysis to develop a new empirical correlation for the prediction of staggered labyrinth seal leakage. The CFD analysis investigates several seal geometries and various operating conditions, which are pressure drop, inlet preswirl, and rotor peripheral speed. To capture the nonlinear dependency of these factors on the flow conditions, three levels each are chosen, except two for the preswirl. The correlation obtained is finally verified numerically, comparing the results of the bulk-flow model with those of a CFD analysis in terms of leakage mass-flow rate, pressure distributions, and Mach number.In conclusion, the rotor dynamic coefficients prediction has been compared with the experimental measurements performed by Kwanka and Ortinger [11], highlighting the more accurate estimation of the cross-force with respect to other bulk-flow models.


## 2. Bulk-Flow Steady-State Solution

The estimation of the steady leakage, pressure, and circumferential velocity distribution represents the so-called zeroth-order solution in the bulk-flow model, as described in [13]. For a fixed pressure drop and assuming, at a first stage, an isenthalpic process, the leakage and the pressure distribution in the seal cavities are determined by solving iteratively the continuity equation, until the mass-flow rate is equal for each control volume (see Figure 2). Then, the enthalpy and the circumferential velocity in each cavity are determined by solving the coupled circumferential momentum and energy equation. Multivariate Newton-Raphson algorithm is employed to find the solution. At the following step, the leakage is calculated one more time, considering the enthalpy variation previously estimated. The iteration continues until the leakage mass flow-rate converges.Figure 2
Bulk-flow control volumes.The continuity, circumferential momentum, and energy zeroth-order equations for each cavity are stated as follows:Continuity Equation (1) m ˙ 01 ≡ ⋯ ≡ m ˙ 0 i ≡ ⋯ ≡ m ˙ 0 N T ≡ m ˙ 0 .Circumferential Momentum Equation (2) m ˙ 0 V 0 i - V 0 i - 1 = τ 0 r i a r i - τ 0 s i a s i .Energy Equation (3) m ˙ 0 h 0 i + V 0 i 2 2 - m ˙ 0 h 0 i - 1 + V 0 i - 1 2 2 = τ 0 r i a r i R s Ω ,where m˙0 is the steady leakage mass-flow rate, V0i and h0i are the steady circumferential velocity and enthalpy for the ith cavity. τ0ri and τ0si are the shear stresses on the rotor and stator parts, whereas ari and asi are the lengths of the rotor and stator parts, respectively.The leakage correlation (m˙i) used in the model is given by the Neumann formula, as described in [14]:(4)m˙i=RsiRsCfiμiHiρi-1Pi-1-Pi,where Cfi and μi are, respectively, the discharge coefficient and kinetic energy carry-over coefficient, Hi is the clearance, and Pi and ρi are the pressure and density in the ith cavity.The shear stresses both on the rotor and on the stator (τ0ri and τ0si) have been evaluated by implementing the Swamee-Jain [15] explicit formula for the estimation of the Darcy friction factor. Compared with the Blasius correlation, which is used in several models [9, 16–19], the Swamee-Jain correlation is more accurate in the calculation of the friction factor, because it takes into account the roughness of the surface and it is validated also for high Reynolds numbers. A comprehensive investigation of the correlation has been performed by Kiijarvi in [20]. The rotor and stator friction factors are defined as(5)fri=116log10⁡er3.7·Dhi+5.74Reri0.9-2fsi=116log10⁡es3.7·Dhi+5.74Resi0.9-2,where Dhi is the circumferential hydraulic diameter, μi is the kinematic viscosity, e is the absolute roughness of the rotor/stator surface, and Ω is the rotational speed of the rotor. The Reynolds numbers are defined for the rotor and stator parts, respectively, as follows:(6)Reri=DhiVRiμiResi=DhiVSiμi.The shear stresses generated by the fluid/structure interaction are evaluated considering both the axial and circumferential velocity, as described in [19] by Wyssmann et al., as follows:(7)τri=ρi2friRsΩ-ViVRiVRi=Ui2+RsΩ-Vi2τsi=ρi2friViVSiVSi=Ui2+Vi2,where Ui is the axial velocity of the fluid in the ith cavity, which is estimated as a function of the calculated mass-flow rate as(8)Ui=m˙iCfiHiρi.


## 3. Preliminary Results

To investigate the accuracy of the leakage correlation used for straight-through labyrinth seals in the case of staggered seals, the zeroth-order solutions are compared with the CFD results. Steady-state CFD analysis is very accurate in the prediction of the leakage, static pressure, and circumferential velocity within the seal; therefore CFD has been considered as the reference for the zeroth-order solution.A detailed description of the CFD analysis is shown in the next chapter of the paper. The seal geometry is shown in Figure3, and the geometrical parameters are reported in Table 1, using the nomenclature shown in Figure 4. The working fluid is the steam.Table 1
Seal geometry parameter.
Seal geometry Shaft radiusRs 189 Radial clearance (mm)s 0.6 Steps-to-casing radial distance (mm)G 2.00 Steps height (mm)H 2.00 Steps width (mm)F 3.00 J-strips width at tip (mm) W 0.25 J-strips pitch (mm) P 4.00 J-strip-to-step axial distance (mm) T 2.75 J-strips number NJ 47 Rotor steps numberNS 24Figure 3
Staggered seal geometry.Figure 4
Parametric view of staggered labyrinth seal.The leakage can be calculated in two manners: considering the kinetic energy carry-over coefficient (Wang’s model), as suggested by Wang et al. [10], or assuming that the kinetic energy carry-over coefficient is equal to the unity for all the constrictions (Childs’ model), as suggested by Childs in [7]. Both models use the empirical Neumann’s correlation (see (4)) for the leakage [14] and Chaplygin’s formula (see (9)) for the discharge coefficient. The kinetic energy carry-over coefficient is calculated using the semiempirical correlation developed by Hodkinson in [21]. (9)Cfi=ππ+2-5·βi+2·βi2where,βi=Pi-1Piγ-1/γ-1(10)μi=11-NJ-1/NJ·Hi/Li/Hi/Li+0.02.The results in terms of leakage mass-flow rate are shown in Figure5, while the swirl velocity and the pressure distribution are shown in Figure 6. The mass-flow rate is overestimated for both bulk-flow models: the relative error is, respectively, equal to +75% and +272% for Childs’ model and Wang’s model. The swirl ratios along the seal cavities are strongly influenced by the mass flow. Even if the asymptotic value of the swirl ratio for Childs’ model is equal to that of the CFD, the values in the first cavities are far from the CFD ones. The predicted static pressure distribution is very similar to Childs’ model.Figure 5
Mass-flow rate comparison.Figure 6
Comparison of the swirl ratio and of the static pressure in the seal cavities evaluated with the CFD, Childs’ model, and Wang’s model.Eventually, Childs’ model is the more accurate one-control volume bulk-flow model in the literature for staggered labyrinth seals; however, the results in terms of leakage, pressure, and circumferential velocity distribution are not accurate enough compared to CFD results. The zeroth-order solution strongly influences the prediction of the seal dynamic coefficients; hence, the authors decide to develop a new empirical correlation for staggered labyrinth seal.


## 4. Design and Analysis of Computer Experiments

Design and analysis of CFD experiments allow several staggered seal geometries and operating conditions to be investigated for the development of a new leakage correlation. A priori, it is not possible to know which are the geometric parameters that influence the leakage flow rate, as well as the operating conditions.Three seal geometries are chosen according to the OEMs design criteria, while additional two are selected in order to extend the domain of the investigation in the case of possible axial shifting of the rotor. Axial excursion of the rotor commonly happens during the start-up condition or even during the variation of the process parameters (i.e., the output power density). The validation of the correlation also for off-design seal operations allows extending the seal dynamic coefficients predictability. The geometric parameters that define the seal geometry are those reported in Table1. Three different rotor diameters have been investigated (220 mm, 500 mm, and 800 mm); the teeth are 20, while the rotor steps are 10. The additional two seal geometries are equal to those with the rotor diameter equal to 500 mm (geometry #4) and 800 mm (geometry #5) but with a different length of the rotor step. The three main geometries are shown in Figure 7; the steam flows from right to left.Figure 7
Three seal geometries investigated.Because CFD analysis is very time-consuming, a factorial DACE is performed, even if not all the combinations of the parameters have been investigated. The list of the CFD analysis performed is shown in Table2.Table 2
Design of computer experiments.
Test Pressure drop Geometry Peripheral speed Preswirl (bar) (m/s) — 1 50 Geom 1 69.1 1.9 2 20 Geom 1 69.1 1.9 3 10 Geom 1 69.1 1.9 4 50 Geom 2 69.1 1.9 5 50 Geom 3 69.1 1.9 6 20 Geom 2 69.1 1.9 7 20 Geom 3 69.1 1.9 8 50 Geom 4 69.1 1.9 9 50 Geom 5 69.1 1.9 10 50 Geom 1 69.1 0.5 11 20 Geom 1 69.1 0.5 12 50 Geom 2 69.1 0.5 13 50 Geom 3 69.1 0.5 14 50 Geom 1 115.2 1.9 15 20 Geom 1 115.2 1.9 16 50 Geom 2 115.2 1.9 17 50 Geom 1 161.3 1.9 18 50 Geom 1 115.2 0.5 19 20 Geom 2 115.2 1.9 20 50 Geom 2 115.2 0.5 21 50 Geom 3 115.2 0.5 22 20 Geom 2 69,1 0,5 23 20 Geom 3 69.1 0.5 24 20 Geom 1 115.2 0.5 25 50 Geom 4 161.3 1.9 26 50 Geom 5 161.3 1.9 27 10 Geom 4 69.1 1.9 28 10 Geom 5 69.1 1.9 29 50 Geom 4 69.1 0.5 30 50 Geom 5 69.1 0.5 31 50 Geom 1 161.3 0.5 32 10 Geom 1 69.1 0.5 33 10 Geom 1 161.3 1.9


## 5. CFD Analysis and Results

The domain discretization has been carried out using Ansys™ [22]. The meshing procedure consists of creating, at a first stage, the 2D computational domain by means of quadrilateral element. After that, the 2D domain is rotated circumferentially to obtain the 3D structured mesh, composed by hexahedra elements. The 3D domain corresponds to a five-degree sector along the circumferential direction. For instance, in the case of the seal geometry with the diameter equal to 220 mm, one element per degree has been used for the mesh rotation, while for the other diameters, two elements per degree are necessary to avoid excessive stretched elements in the circumferential direction.It is worth mentioning that the discretization has been refined close to the wall to achieve a low Reynolds approach in the boundary layer (see Figure8). Details concerning the discretization have been summarized in Table 3.Table 3
Mesh parameters.
Mesh parameters Height of first layer 1 ⋅ 10 - 3 mm Cell expansion ratio 1.2 Number of elements 1.25 million (∅ 220 mm) 2.50 million (∅ 500 and 800 mm)Figure 8
Details of the labyrinth seal discretization.Among the several approaches to solve numerically the Navier-Stokes equations, the most employed in the turbomachinery field is the Reynolds Averaged Navier-Stokes (RANS) due to the good compromise between accuracy and computational time required. The IAWPS steam tables have been employed to properly model the fluid properties [23].Turbulence has been modelled using theSSTk-ω model, firstly introduced by Menter [24] and widely tested in turbomachinery application in the last decade. This turbulence model is associated with the automatic near wall treatment, where a smooth transition between the low Reynolds number approach and the standard wall functions is adopted. However, due to the low value of y+ achieved, using the grids employed in this paper, this treatment corresponds to the low Reynolds model; hence, the viscous sublayer is not modelled but directly solved. This approach is surely preferable due to the final purpose of the present analysis, where the exact amount of flow rate for a given pressure drop must be assessed through the clearances.The boundary condition scheme is hereafter introduced in Figure9. The same set of conditions are imposed to all the test cases. The value of inlet total pressure is changed at each new iteration to target the desired pressure drop; moreover, the circumferential velocity component is adjusted to obtain the desired preswirl ratio.Figure 9
Boundary conditions scheme.The results of some simulations are shown in Figures10 and 11. Static pressure and temperature contours are reported in Figure 10. Dimensionless values are obtained using the downstream value equaling the unity. Furthermore, streamlines in the meridional plane are illustrated in Figure 11, focusing on the inlet and outlet locations. The flow structure before the first clearance is different from the one occurring upstream the consecutive teeth.Figure 10
Nondimensional static pressure (a) and temperature (b) contours.
(a) (b)Figure 11
Meridional velocity coloured streamlines at inlet (a) and outlet (b) location.
(a) (b)


## 6. Discharge Coefficients Correlation

The objective of the paper is the development of a new leakage correlation to accurately predict the leakage mass-flow rate and the pressure distribution in the seal cavities in steady-state condition (centred rotor). In a nozzle or other constrictions, the discharge coefficient is defined as the ratio of the actual mass-flow rate at the discharge end of the constriction to that in which an ideal working fluid expands with the same initial and exit conditions. Starting from the definition of the discharge coefficient, it is possible to define the leakage formula, which is equal to the product between the discharge coefficient and the ideal mass-flow rate. The ideal mass-flow rate in a constriction is derived from the Bernoulli equation [25], assuming that the outlet area is much smaller than the inlet one. The coefficient of discharge is defined as(11)Cfi=m˙i+1Ai2ρiP0i-1-P0i.Since the mass-flow rate, pressure, and density can be evaluated by the CFD, the discharge coefficient for each cavity can be evaluated. Using the CFD data of all the tests investigated by the DACE, the idea is to develop a new correlation based on the geometric parameters and on the operating conditions.The authors have observed that the discharge coefficients have different values depending on the constriction type (long teeth or short teeth coupled with the rotor step); however, the values along the seal orifices almost remain constant, for both the constriction types, expect for the first tooth (see Figure12). Since the discharge coefficient depends on the angle at which the fluid approaches to the tooth, the different discharge coefficient is likely due to the fluid, which approaches to the long teeth and to the short ones in two different manners. For instance, as can be observed in Figure 12, for test 1 shown in Table 2, the discharge coefficient in correspondence of the short teeth is about 0.3, whereas for long teeth it is practically constant and equal to 0.48. For the first tooth, the discharge coefficient is about 0.6.Figure 12
Discharge coefficients in the seal orifice for the test number 4 (seal diameter = 500 mm, pressure drop = 50 bar, inlet preswirl = 1.9, and rotor peripheral velocity = 66.19 m/s).After several investigations, the authors have observed that two geometrical parameters mainly influence the discharge coefficient. These parameters are listed in Table4. The parameter s/H+G is related to the ratio between the height of the constriction and the height of the seal. The other parameters are related to the axial position and length of the rotor step. In addition, one parameter regarding the operating conditions has been considered. From a primary screening of the CFD results, the authors observed that the only operating condition parameter that significantly affects the leakage is the pressure drop. The inlet preswirl and the rotor peripherical speed have practically no influence on the mass flow and therefore on the pressure distribution.Table 4
Nomenclature of the geometric parameters defined in Table1.
Geometric parameters for thej-strip Geometric parameters for the rotor step x j = s H + G x s = s H + G  y j = 1 - F + T - P P - W y s = 1 - T P - WThe authors assume that the cross-dependence terms between the geometrical parameters and operating conditions are negligible. Thus, at a first stage, the discharge coefficients have been fitted with a polynomial as a function of the geometrical parameters for both thej-strips teeth and the short j-strips teeth combined with the rotor steps. At a second stage, the residual between the discharge coefficients evaluated using the CFD data and the discharge coefficients calculated with the polynomial in (12)-(13) have been fitted using the operating parameter as shown in (14)-(15). The superscripts and subscripts “j”indicate the long j-strips (long teeth), whereas the superscripts and subscripts “s”indicate the shortj-strips combined with the rotor steps.(12)Cfj=q0+q1xj+q2yj+q3xj2+q4xjyj(13)Cfs=p0+p1xs+p2ys+p3xs2+p4xsys(14)Residualj=q5+q6ρinPin-Pout(15)Residuals=p5+p6ρinPin-Pout.For confidentiality reasons, the values of the coefficientsq0,…,q6 and p0,…,p6, shown in (12)–(15), are not reported in the paper.Figures13 and 14 show the polynomial surface that fits the discharge coefficients evaluated with the CFD analysis for both short teeth (Figure 13) and long teeth (Figure 14). The blue points represent the average value of the discharge coefficients, in correspondence of the rotor step, for each CFD analysis. The CFD points are well fitted, since the adjusted R2 [26] is about 0.97 for the short teeth and 0.98 for the long ones, whereas the root mean square error (RMSE) [26] is equal to 0.01406 and 0.01399, respectively. These statistical parameters highlight the goodness of fit of the data.Figure 13
Surface fitting of the discharge coefficients for the short teeth with the geometrical parameters.Figure 14
Surface fitting of the discharge coefficients for the long teeth with the geometrical parameters.Furthermore, Figures15 and 16 show the fitting of the residual by the operating parameter (ρinPin-Pout). The blue points represent the average value of the residual, in correspondence of the rotor step, for each CFD analysis. The adjusted R2 is about 0.9689 and the RMSE is 0.01406.Figure 15
Curve fitting of the residual for the short teeth with the operating parameter.Figure 16
Curve fitting of the residual for the long teeth with the operating parameter.The discharge coefficient in correspondence of the first tooth has been considered constant for all the tests and it seems to be uncorrelated from the geometric parameters and from the pressure drop.


### 6.1. Prediction of the Leakage and Pressure Distribution

The new discharge coefficient correlation has been implemented in the bulk-flow model. The bulk-flow leakage mass-flow estimations compared to the CFD leakage are shown in Figure17. The bulk-flow leakage is in the range of ±6% of relative error (the blue region defines the ±10% of relative error).Figure 17
Nomenclature of the geometric parameters.Moreover, the pressure distribution predicted with the bulk-flow model is very similar to that evaluated with the CFD. In Figure18(a), the pressure shows a “zig-zag” behaviour along the seal, which is due to the nonconstant discharge coefficients (see Figure 18(b)).Figure 18
Static pressure distribution in the seal cavities and discharge coefficients in the seal orifices for the test number 1 (seal diameter = 220 mm, pressure drop = 50 bar, inlet preswirl = 1.9, and rotor peripheral velocity = 66.19 m/s).
(a) (b)


### 6.2. Prediction of the Circumferential Velocity

In this section, the authors have compared the fluid-wall shear stresses estimated by the bulk-flow model with those computed using the CFD. For the sake of brevity, only the results related to test case 1 are reported (see Figure19) considering both cases with and without the axial velocity component in the bulk-flow model. The accuracy of the results is still long valid for the other tests.Figure 19
Rotor and stator shear stresses as a function of the axial coordinate for case 1 shown in Table2.The CFD shear stresses are averaged along the length of each cavity to be compared with the bulk-flow results. In Figure20, the circumferential velocity profile is shown. Typically, in the 1CV bulk-flow model (see bulk-flow models in [9, 10, 18]), the shear stresses are evaluated considering only the circumferential velocity, whereas the authors have considered also the contribution of the axial velocity (see (7)) as suggested in [19].Figure 20
Circumferential velocity as a function of the axial coordinate for case 1 shown in Table2.


### 6.3. Numerical Verification of the New Correlation

The new correlation has been numerically verified using the results of the staggered seal shown in Figure3, which has about twice the teeth of the DACE seal geometries and moreover the operating conditions are far from those investigated in the DACE. The relative error between the leakage evaluated with the bulk-flow model and the one with the CFD is equal to 10% (see Figure 21). Moreover, the pressure distribution is accurately reproduced, and the sonic flow condition is reached as in the CFD analysis (see Figure 23). The results of Childs’ model and Wang’s model are reported in Figures 21 and 22.Figure 21
Mass-flow rate comparison.Figure 22
Comparison of the swirl ratio and of the static pressure in the seal cavities evaluated with the CFD, Childs’ model, Wang’s model, and authors’ model.Figure 23
Mach number evaluated using the CFD and the authors’ model.


## 6.1. Prediction of the Leakage and Pressure Distribution

The new discharge coefficient correlation has been implemented in the bulk-flow model. The bulk-flow leakage mass-flow estimations compared to the CFD leakage are shown in Figure17. The bulk-flow leakage is in the range of ±6% of relative error (the blue region defines the ±10% of relative error).Figure 17
Nomenclature of the geometric parameters.Moreover, the pressure distribution predicted with the bulk-flow model is very similar to that evaluated with the CFD. In Figure18(a), the pressure shows a “zig-zag” behaviour along the seal, which is due to the nonconstant discharge coefficients (see Figure 18(b)).Figure 18
Static pressure distribution in the seal cavities and discharge coefficients in the seal orifices for the test number 1 (seal diameter = 220 mm, pressure drop = 50 bar, inlet preswirl = 1.9, and rotor peripheral velocity = 66.19 m/s).
(a) (b)


## 6.2. Prediction of the Circumferential Velocity

In this section, the authors have compared the fluid-wall shear stresses estimated by the bulk-flow model with those computed using the CFD. For the sake of brevity, only the results related to test case 1 are reported (see Figure19) considering both cases with and without the axial velocity component in the bulk-flow model. The accuracy of the results is still long valid for the other tests.Figure 19
Rotor and stator shear stresses as a function of the axial coordinate for case 1 shown in Table2.The CFD shear stresses are averaged along the length of each cavity to be compared with the bulk-flow results. In Figure20, the circumferential velocity profile is shown. Typically, in the 1CV bulk-flow model (see bulk-flow models in [9, 10, 18]), the shear stresses are evaluated considering only the circumferential velocity, whereas the authors have considered also the contribution of the axial velocity (see (7)) as suggested in [19].Figure 20
Circumferential velocity as a function of the axial coordinate for case 1 shown in Table2.


## 6.3. Numerical Verification of the New Correlation

The new correlation has been numerically verified using the results of the staggered seal shown in Figure3, which has about twice the teeth of the DACE seal geometries and moreover the operating conditions are far from those investigated in the DACE. The relative error between the leakage evaluated with the bulk-flow model and the one with the CFD is equal to 10% (see Figure 21). Moreover, the pressure distribution is accurately reproduced, and the sonic flow condition is reached as in the CFD analysis (see Figure 23). The results of Childs’ model and Wang’s model are reported in Figures 21 and 22.Figure 21
Mass-flow rate comparison.Figure 22
Comparison of the swirl ratio and of the static pressure in the seal cavities evaluated with the CFD, Childs’ model, Wang’s model, and authors’ model.Figure 23
Mach number evaluated using the CFD and the authors’ model.


## 7. Prediction of Staggered Labyrinth Seal Rotor Dynamic Coefficients

### 7.1. Bulk-Flow First-Order Mathematical Treatment

The bulk-flow model, developed by the authors, is based on the one-control volume bulk-flow model described in [13, 27–29] and here customized for staggered labyrinth seal. Each cavity of the seal is modelled by one-control volume. The model includes the steam properties evaluated with the IAPWS tables [23]. The model is based on the perturbation analysis, as described in [7].The fluid dynamics within the seal are governed by the continuity equation, circumferential momentum equation, and energy equation (only in the zeroth-order problem). The first two governing equations are stated as follows:(i)
Continuity equation [9]:(16)∂∂tρiAi+∂∂ϑρiAiViRs+m˙i+1-m˙i=0.(ii)
Circumferential momentum equation [9]:(17)∂∂tρiAiVi+∂∂ϑρiAiVi2Rs+m˙i+1Vi-m˙iVi-1=-AiRs∂Pi∂ϑ+τriari-τsiasi.


### 7.2. Perturbation Analysis

The rotor whirls within the seal in thex and y directions based on the equations of motion are as follows:(18)xt=ϵr0cos⁡ωt,yt=ϵr0sin⁡ωt,where ϵr0 is the radius of the orbit.By using the complex notation, the perturbation of the clearanceH1t,ϑ can be written as follows:(19)H1t,ϑ=-Rer0ejϑ+ωt+r0ejϑ-ωt.The solution of the first-order equation represented by the perturbed terms of pressure and velocity can be given by the same mathematical expression of the perturbed clearance due to the linearity of the system as follows:(20)P1it,ϑ=ReP1i+ejϑ+ωt+P1i-ejϑ-ωtV1it,ϑ=ReV1i+ejϑ+ωt+V1i-ejϑ-ωt,where Pi+,Pi-,Vi+, and Vi- represent the amplitude of the harmonic solution. Four linearized algebraic equations can be obtained for each cavity as follows:(21)Ai-1Yi-1+Ai0Yi+Ai+1Yi+1=Bi.The matrices in (21) are defined in the appendix. Lastly, the matrices of each cavity can be assembled in a 4N-1×4N-1 band matrix, where N is the number of cavities. The linear system has 4N-1 unknowns, that is, Pi+,Pi-,Vi+, and Vi-, for each ith cavity. Then, the complete seal linear system can be obtained as follows:(22)⋱⋮⋮⋮⋮⋮⋰⋯Ai-1-1Ai-10Ai-1+100⋯⋯0Ai-1Ai0Ai+10⋯⋯00Ai+1-1Ai+10Ai+1+1⋯⋰⋮⋮⋮⋮⋮⋱⋮Yi-2Yi-1YiYi+1Yi+2⋮=⋮Bi-1BiBi+1⋮.


### 7.3. Dynamic Coefficients

By considering the contribution of the shear stresses on the rotor surface, the dynamic force acting on the rotor [18] is given by(23)F=Fx+jFy=-ϵRs∑i=1N-1Li∫02πP1iejϑ-jariτr1iejϑ.The dynamic coefficients (K,k,C,c) of the seal are calculated by means of the following equations:(24)K=12ReZ++Z-k=-12ImZ++Z-C=-12ΩImZ+-Z-c=-12ΩReZ+-Z-,where Z+ and Z- are complex numbers that correspond to the dynamic force acting on the rotor and are produced by the forward and the backward orbit, respectively, which can be given as follows:(25)Z+=πRs∑i=1N-1Li∫02πP1i+1-jari∂τri∂P1i-jari∂τri∂V1iV1i++jari∂τri∂H1Z-=πRs∑i=1N-1Li∫02πP1i-1-jari∂τri∂P1i-jari∂τri∂V1iV1i-+jari∂τri∂H1.


### 7.4. Numerical Results

A preliminary investigation on the dynamic behaviour of staggered labyrinth seals has been performed in the paper. The numerical results are compared with the experimental measurements performed by Kwanka and Ortinger in [11]. The radial force generated by a fourteen-cavity staggered labyrinth seal was obtained by integrating the pressure measured along the circumferential direction. The maximum circumferential velocity at the seal inlet is equal to 22 m/s. The seal geometry and the operating conditions are reported in Tables 5 and 6.Table 5
Geometry parameters of the seal used in the experimental tests.
Shaft radius R s 111 [mm] Radial clearance s 0.50 [mm] Steps-to-casing radial distance G 4.00 [mm] Steps height H 5.00 [mm] Steps width F 10.50 [mm] J-strips width at tip W 0.30 [mm] J-strips pitch P 10.50 [mm] J-strip-to-step axial distance T 5.00 [mm] J-strips number N J 15 Rotor steps number N S 8Table 6
Operating conditions of the experimental tests.
Pressure drop Δ P 0.8 [bar] 1.80 [bar] 3.00 [bar] 4.60 [bar] Outlet pressure P o u t 1.01 [bar] Inlet temperature T i n 693 [K] Rotational speed Ω 3000 [rpm]Figure24 shows the measured cross-force as a function of the pressure drop, considering the rotor eccentricity equal to 0.4 mm. The geometry of the seal and the operating condition are far from those used for the development of the leakage correlation, but the prediction obtained by the authors’ model is more accurate than the numerical results reported in the paper by Kwanka and Ortinger.Figure 24
Comparison between the cross-forces evaluated experimentally and those predicted by three different bulk-flow models.


## 7.1. Bulk-Flow First-Order Mathematical Treatment

The bulk-flow model, developed by the authors, is based on the one-control volume bulk-flow model described in [13, 27–29] and here customized for staggered labyrinth seal. Each cavity of the seal is modelled by one-control volume. The model includes the steam properties evaluated with the IAPWS tables [23]. The model is based on the perturbation analysis, as described in [7].The fluid dynamics within the seal are governed by the continuity equation, circumferential momentum equation, and energy equation (only in the zeroth-order problem). The first two governing equations are stated as follows:(i)
Continuity equation [9]:(16)∂∂tρiAi+∂∂ϑρiAiViRs+m˙i+1-m˙i=0.(ii)
Circumferential momentum equation [9]:(17)∂∂tρiAiVi+∂∂ϑρiAiVi2Rs+m˙i+1Vi-m˙iVi-1=-AiRs∂Pi∂ϑ+τriari-τsiasi.


## 7.2. Perturbation Analysis

The rotor whirls within the seal in thex and y directions based on the equations of motion are as follows:(18)xt=ϵr0cos⁡ωt,yt=ϵr0sin⁡ωt,where ϵr0 is the radius of the orbit.By using the complex notation, the perturbation of the clearanceH1t,ϑ can be written as follows:(19)H1t,ϑ=-Rer0ejϑ+ωt+r0ejϑ-ωt.The solution of the first-order equation represented by the perturbed terms of pressure and velocity can be given by the same mathematical expression of the perturbed clearance due to the linearity of the system as follows:(20)P1it,ϑ=ReP1i+ejϑ+ωt+P1i-ejϑ-ωtV1it,ϑ=ReV1i+ejϑ+ωt+V1i-ejϑ-ωt,where Pi+,Pi-,Vi+, and Vi- represent the amplitude of the harmonic solution. Four linearized algebraic equations can be obtained for each cavity as follows:(21)Ai-1Yi-1+Ai0Yi+Ai+1Yi+1=Bi.The matrices in (21) are defined in the appendix. Lastly, the matrices of each cavity can be assembled in a 4N-1×4N-1 band matrix, where N is the number of cavities. The linear system has 4N-1 unknowns, that is, Pi+,Pi-,Vi+, and Vi-, for each ith cavity. Then, the complete seal linear system can be obtained as follows:(22)⋱⋮⋮⋮⋮⋮⋰⋯Ai-1-1Ai-10Ai-1+100⋯⋯0Ai-1Ai0Ai+10⋯⋯00Ai+1-1Ai+10Ai+1+1⋯⋰⋮⋮⋮⋮⋮⋱⋮Yi-2Yi-1YiYi+1Yi+2⋮=⋮Bi-1BiBi+1⋮.


## 7.3. Dynamic Coefficients

By considering the contribution of the shear stresses on the rotor surface, the dynamic force acting on the rotor [18] is given by(23)F=Fx+jFy=-ϵRs∑i=1N-1Li∫02πP1iejϑ-jariτr1iejϑ.The dynamic coefficients (K,k,C,c) of the seal are calculated by means of the following equations:(24)K=12ReZ++Z-k=-12ImZ++Z-C=-12ΩImZ+-Z-c=-12ΩReZ+-Z-,where Z+ and Z- are complex numbers that correspond to the dynamic force acting on the rotor and are produced by the forward and the backward orbit, respectively, which can be given as follows:(25)Z+=πRs∑i=1N-1Li∫02πP1i+1-jari∂τri∂P1i-jari∂τri∂V1iV1i++jari∂τri∂H1Z-=πRs∑i=1N-1Li∫02πP1i-1-jari∂τri∂P1i-jari∂τri∂V1iV1i-+jari∂τri∂H1.


## 7.4. Numerical Results

A preliminary investigation on the dynamic behaviour of staggered labyrinth seals has been performed in the paper. The numerical results are compared with the experimental measurements performed by Kwanka and Ortinger in [11]. The radial force generated by a fourteen-cavity staggered labyrinth seal was obtained by integrating the pressure measured along the circumferential direction. The maximum circumferential velocity at the seal inlet is equal to 22 m/s. The seal geometry and the operating conditions are reported in Tables 5 and 6.Table 5
Geometry parameters of the seal used in the experimental tests.
Shaft radius R s 111 [mm] Radial clearance s 0.50 [mm] Steps-to-casing radial distance G 4.00 [mm] Steps height H 5.00 [mm] Steps width F 10.50 [mm] J-strips width at tip W 0.30 [mm] J-strips pitch P 10.50 [mm] J-strip-to-step axial distance T 5.00 [mm] J-strips number N J 15 Rotor steps number N S 8Table 6
Operating conditions of the experimental tests.
Pressure drop Δ P 0.8 [bar] 1.80 [bar] 3.00 [bar] 4.60 [bar] Outlet pressure P o u t 1.01 [bar] Inlet temperature T i n 693 [K] Rotational speed Ω 3000 [rpm]Figure24 shows the measured cross-force as a function of the pressure drop, considering the rotor eccentricity equal to 0.4 mm. The geometry of the seal and the operating condition are far from those used for the development of the leakage correlation, but the prediction obtained by the authors’ model is more accurate than the numerical results reported in the paper by Kwanka and Ortinger.Figure 24
Comparison between the cross-forces evaluated experimentally and those predicted by three different bulk-flow models.


## 8. Conclusion

The paper focuses on the numerical modelling of staggered labyrinth seals to predict the rotor dynamic coefficients. By considering the one-control volume bulk-flow model, a new leakage correlation has been developed to estimate accurately the zeroth-order solution: leakage mass-flow rate, static pressure distribution, circumferential velocity distribution, and temperatures.Design and analysis of CFD experiments have been used to develop the leakage correlation using four design factors: the geometry, the pressure drop, the inlet preswirl, and the rotor peripherical speed. Only the geometry and the pressure drop have been confirmed to be effective on the estimation of the leakage and static pressure.The procedure used by the authors for the development of the correlation relies on the definition of the discharge coefficient. The discharge coefficients have different values in correspondence of the long teeth with respect to the short ones. Polynomial functions have been used to calculate the discharge coefficients as a function of the seal geometry parameters and of the pressure drop. Once these coefficients have been evaluated, the static pressure in the cavities is estimated till the leakage converges.Finally, a detailed description of the one-control volume bulk-flow model for staggered labyrinth seals has been reported in the paper. The model allows the rotor dynamic coefficients to be estimated. Numerical results are compared to experimental measurements by Kwanka and Ortinger. The results of the model developed by the authors are more accurate than the other numerical models.



---

*Source: 9357249-2018-05-15.xml*