# Aeroelastic stability of a labyrinth seal coupled to a flexible stator, with a one control-volume bulk-flow model including temperature fluctuations

M. Fleury<sup>1,2</sup>, F. Thouverez<sup>1</sup>, L. Blanc<sup>1</sup>, P. Girard<sup>2</sup>

<sup>1</sup> Ecole Centrale de Lyon, Laboratoire de Tribologie et Dynamique des Systèmes, 69134 Ecully, France

e-mail: marine.fleury@ec-lyon.fr

<sup>2</sup> Safran Aircraft Engines Villaroche, 77550 Moissy-Cramayel, France

## Abstract

In turbomachinery, labyrinth seals are key components to reach high performance of the engine. To that end, manufacturers of sealing systems tend to reduce the radial clearance between stationary and rotating parts of the turbine to increase efficiency. They also reduce the thickness of structural parts to get lighter designs. These new designs could lead to aeroelastic instability issues. To prevent such instability problems, an accurate prediction of the aeroelastic effective damping is necessary. Stability criteria used in the industry are mainly based on empirical observations and show limitations. A more representative prediction can be done using CFD calculations, but it remains computationally expensive. In this paper, a semi-analytical model is developed to analyze the fluid-structure interaction inside a multi-fins straight labyrinth seal. Fluid equations are written using a one control-volume bulk flow model, and assuming temperature fluctuations. This last feature is rarely addressed in the literature and is expected to provide better fidelity.

## 1 Introduction

In gas turbomachinery, labyrinth seals are used to control the leakage between high-pressure and low-pressure regions. They are located between the rotating and stationary parts of the structure. They are composed of several radial fins, forming annular cavities with gas recirculation inside. The axial leakage is allowed by a clearance of a few hundreds of microns, that leads to pressure drops all through the seal. New designs of seals present more deformable structural parts and even reduced clearances. Under certain operating conditions, this can lead to aeroelastic instabilities, such as flutter, which entails fatigue cracks. Such phenomenon has been historically observed on different labyrinth seals [1, 2] and becomes nowadays even more sensitive. Thus, it is necessary to achieve an accurate prediction of the seal behavior, to prevent damage.

Sealing systems stability has been studied for many years. In the 1960s, Alford [3] and Ehrich [4] are the first to highlight different parameters of influence. Radial clearance and fin-support position (high-pressure or low-pressure side) are shown to have both critical impact. Alford also proposes a stability criterion [1], considering radial deflection, and based on empirical results. This model neglects a lot of physical parameters, such as radial clearance, and cannot be used to predict flutter in seals. In the early 1980s, Lewis et al. [2] and Abbott [5] propose an improvement of Alford model. Abbott's criterion suggests that in the case of a high-pressure side support, the seal is stable as long as the mechanical frequency is lower than the acoustic frequency, at same nodal diameter. The phenomenon is reversed in the case of a low-pressure side support. This model, even if employed in the industry, assumes a weak coupling between fluid and structure, and remains poorly predictive in some cases.
Then, in the 1990s, a semi-analytical model, based on Iwatsubo work [6], was proposed by Childs and Scharrer [7]. The continuity and momentum equations are written for a control volume in each cavity, assuming isothermal flow along the seal. Resolution is carried out in two steps, assuming small perturbations around the steady state, leading to linearized equations for the resolution of the unsteady problem. The main flaw of this model is to consider a whole rigid structure, whereas new geometries of aircraft engines present thin parts, that might undergo deformations.

Then, developments in CFD calculations allowed better understanding of labyrinth seal behavior. In their different works, Rhode [8], Kwanka, Tsukuda and Hirano [9, 10] perform predictions of rotordynamic coefficients of the seal. Although reliable, CFD approach remains computationally expensive with respect to bulk-flow model, and cannot be used for pre-design studies. Recently, Cangioli et al. [11, 12] propose a one control-volume bulk-flow model involving thermodynamics, by introducing the energy equation at steady-state, with an enthalpic formulation. Corral and Vega [13, 14] also propose a stability criterion based on an explicit expression of the work per cycle: stability is assessed through the comparison between discharge time and acoustic propagation time along cavity.

The purpose of this paper is to propose a new semi-analytical model, rewritting the Navier-Stokes continuity and momentum equations, adding the energy equation to take the effect of temperature fluctuations into account, at both steady and unsteady states - this last point being the most remarkable novelty. One considers a simplified but representative design of phenomena occuring in labyrinth seals. For confidentiality issues, tested geometries presented are not real industrial seals. The model includes a rigid shaft, supporting the seal, and a flexible stator. It allows for different configurations of seal (multi-cavities, shaft rotation, inlet circumferential velocity). The structural motion induces pressure and velocity fluctuations inside the cavities of the labyrinth seal. The structure dynamics is described using a cylindrical shell model taking into account the support position (high-pressure side or low-pressure side, hereinafter 'HP' or 'LP'). A monolithic resolution of fluid and structural equations is carried out, leading to a strongly coupled model. While models historically used to be based on isothermic approaches [7, 15], or adiabatic ones [13], the approach in this paper is more general. With this model, these two limit cases can be identified and a better interpretation of seal behavior is obtained.

# 2 Fluid governing equations

## 2.1 Problem definition and hypotheses

This model is based on the one-control volume bulk-flow model of Childs and Scharrer [7]. It includes multi-cavities and shaft rotation. Each cavity $i$ is represented by one control-volume (Fig. 1). The fluid characteristics assumptions are as follows :

1. Air is supposed to be a perfect gas ; consequently, the ideal gas law will be used to describe the fluid state;

2. Air is a compressible gas, its density will vary in time and space.

3. A viscous flow is considered : the circumferential shear stresses will be calculated using a Colebrook model [16];

4. Pressure variations inside a cavity is negligible with respect to pressure variations between the cavities;

5. A choked flow condition is assumed on the last fin of the labyrinth seal.

Thermodynamically, rotor and stator are supposed to be adiabatic walls. The main contribution of this work is to introduce the energy equation to the model, allowing to consider the influence of temperature fluctuations through the seal. Indeed, this new approach enables to be more representative of industrial configurations, where shaft rotation might introduce heating due to wall friction. Moreover, a formulation in temperature has been preferred to an enthalpic one [11], because it enables a better physical interpretation of the results, as it is a measurable variable. The temperature fluctuation is taken into account at both steady and unsteady
states. The implementation of energy equation at unsteady state represents the main improvement of this model with respect to other one control-volume models such as Cangioli's one [11].

<figure><img src="image_1.png" /><figcaption>Figure 1: Control Volume of the i-th cavity with thermodynamic and kinematic variables</figcaption></figure>

Fluid behavior is described by the Navier-Stokes equations in each cavity, that leads to a highly non-linear system. Variables describing the system in a cavity i (Fig. 1) are: pressure $P_i$, temperature $T_i$, circumferential velocity $U_i$, fluid density $\rho_i$, radial displacement $H_i$ and mass flow rate $\dot{m}_i$. Nominal radial clearance $C_r$ denotes the distance between teeth tip and stator at rest.

## 2.2 Navier-Stokes equations

The model governing fluid behavior within a cavity is inspired by Childs work [17, 7] and based on Navier-Stokes equations. Unlikely to Childs hypotheses, one does not consider an isothermic flow. For each control volume of a cavity i, the fluid integrated governing equations are the following:

* Continuity equation:

$$ \frac{\partial}{\partial t}(\rho_i A_i) + \frac{1}{R_s} \left[ \frac{\partial}{\partial \theta}(U_i A_i \rho_i) \right] + \dot{m}_{i+1} - \dot{m}_i = 0 \quad (1) $$

* Circumferential momentum equation:

$$ \frac{\partial}{\partial t}(\rho_i U_i A_i) + \frac{1}{R_s} \frac{\partial}{\partial \theta}(\rho_i A_i U_i^2) + U_i \dot{m}_{i+1} - U_{i-1} \dot{m}_i = - \frac{A_i}{R_s} \frac{\partial P_i}{\partial \theta} + \tau_{r,i} a_{r,i} L - \tau_{s,i} a_{s,i} L \quad (2) $$

The fluid shear stress $\tau_{r,i}$ and $\tau_{s,i}$ (respectively for the rotor and the stator) are expressed using a Colebrook formulation [16]. This formulation as been preferred to a Blasius one as it shown better results, especially for high Reynolds number [15].

* Energy equation:

$$ \begin{aligned} & \frac{\partial}{\partial t} \left( \left(C_p T_i + \frac{U_i^2}{2}\right) \rho_i A_i \right) + \frac{\partial}{\partial \theta} \left( \left(C_p T_i + \frac{U_i^2}{2}\right) \frac{\rho_i A_i U_i}{R_s} \right) + \dot{m}_i \left(C_p T_i + \frac{U_i^2}{2}\right) - \dot{m}_{i-1} \left(C_p T_{i-1} + \frac{U_{i-1}^2}{2}\right) \\ & = \tau_{r,i} a_{r,i} L_i R_s \omega + \frac{\partial}{\partial t} (P_i A_i) \end{aligned} \quad (3) $$

This last equation is based on the White's enthalpic formulation [18], assuming we have:

$u_i = h_i - \frac{P_i}{\rho_i}$, and $h_i = C_p T_i + h_0$, where $h_0$ is an arbitrary constant.

Here, a temperature formulation is chosen as it is a physical and measurable variable.

Finally, to represent fluid leakage behavior in axial direction, a mass-flow equation is necessary. This last
one induces pressure coupling between cavity i and cavity i − 1. Here, mass flow-rate is modeled using a modified form of the Neumann’s empirical formula [19]:

$$
\dot{m}_i = \mu_0 \mu_i H_i \sqrt{\frac{P_{i-1}^2 - P_i^2}{RT_i}}, \qquad (4)
$$

with $\mu_i$ being the kinetic carry-over coefficient, and $\mu_0$ the discharge coefficient.

The system of equations (1), (2), (3), (4) provides a fully-coupled and highly non-linear problem to describe
fluid behavior within the seal. Small perturbation method is carried out below.

2.3 Linearized problem

To predict system’s dynamics, equations (1), (2), (3), (4) are developed using a perturbation method. It assumes small variations around a mean value. In cavity i, thermodynamic and kinematic variables are written as sums of a mean term (zeroth order) and a fluctuation term (first order), where ϵ denotes the perturbation parameter:

$$
\begin{align*}
& P_i = P_{0i} + \epsilon P_i', \\
& T_i = T_{0i} + \epsilon T_i', \\
& U_i = U_{0i} + \epsilon U_i', \\
& \rho_i = \rho_{0i} + \epsilon \rho_i', \\
& H_i = C_r + \epsilon H_i', \\
& \dot{m}_i = \dot{m}_{0i} + \epsilon \dot{m}_i'.
\end{align*}
\tag{5}
$$

Equations (5) can be solved separately, cascading linear algorithms to solve zeroth and first epsilon orders.

2.3.1 Steady-state problem

Equations at steady-state are established here to provide mean distributions of pressure, circumferential
velocity and temperature within the seal. At steady-state, Navier-Stokes equations are reduced to :

• Continuity equation:

$$
\dot{m}_{0i+1} = \dot{m}_{0i} = \dot{m}_0 \quad \forall i \tag{6}
$$

• Circumferential momentum equation:

$$
\dot{m}_0 (U_{0i} - U_{0i-1}) = \tau_{r0,i} a_{r0} L - \tau_{s0,i} a_{s0} L \quad (7)
$$

with

$$
\tau_{r0,i} = \frac{1}{2} \frac{P_{0i}}{RT_{0i}} (R_s \omega - U_{0i})^2 n_r \left( \frac{|R_s \omega - U_{0i}| D_h}{\nu} \right)^{m_r} \operatorname{sgn}(R_s \omega - U_{0i}) \quad (8)
$$

and

$$
\tau_{s0,i} = \frac{1}{2} \frac{P_{0i}}{RT_{0i}} U_{0i}^2 n_r \left( \frac{U_{0i} D_h}{\nu} \right)^{m_r} \operatorname{sgn}(U_{0i}) \quad (9)
$$

• Energy equation:

$$
\dot{m}_0 \left(C_p T_{0i} + \frac{U_i^2}{2}\right) - \dot{m}_0 \left(C_p T_{0i-1} + \frac{U_{i-1}^2}{2}\right) = \tau_{r0,i} R_s \omega L \quad (10)
$$

To solve this system, we need to provide fluid input parameters and evaluate the mass flow rate assuming a choked flow hypothesis [20].

Problem inputs are inlet and outlet pressures, inlet temperature, inlet velocity, and shaft rotation velocity. To obtain solution at steady-state, a first approximation of the flow-rate has to be achieved. To that end, choked
flow conditions are assumed on the last fin [20]. This enables to define the ratio between the outlet pressure $P_S$ and the last cavity pressure $P_{0N-1}$ as 0.528 [15]. then, the steady mass flow-rate in the last cavity is given by Fliegner's formula [21]:

$$ \dot{m}_0 = \frac{0.510\mu_0}{\sqrt{RT_0}} P_{0N-1}C_r. \qquad (11) $$

By means of this equation and Neumann's one (eq. (4)), we get a closed system of fully-coupled equations.

The addition of the energy equation induces a coupling between the three equations at zeroth order. Moreover, the cavities are linked to each other : the state of the cavity $i$ clearly depends on the cavities $i-1$ and $i+1$. This implies a monolithic resolution of the steady state problem, involving a Newton-Raphson algorithm. To perform the resolution, a well chosen initial point is necessary to prevent convergence issues. The solution result of the Navier-Stokes equations without the energy equation is chosen as initial point.

It is to be noticed that the choked flow condition is checked after each iteration of the whole procedure, and a correction on the flow-rate is made until satisfactory convergence on the pressure distribution. This first step provides pressure, circumferential velocity and temperature distributions at steady-state.

### 2.3.2 First-Order problem

As stated before, the governing equations (1), (2), (3) are expanded using a perturbation analysis method (eq. (5)). The linearized equations at order $\epsilon$ are presented in Appendix A of this paper.

Due to the geometry, the solution is expanded as Fourier series (according to the angle $\theta$); each fluctuation variable $y'_i$ is written as follows:

$$ y'_{i} = \sum_{n=1}^{J} [y_{n,i}^{c} \cos(n\theta) + y_{n,i}^{s} \sin(n\theta)], \qquad (12) $$

where $n$ is the spatial harmonic and $J$ is the truncature order.

Pressure, circumferential velocity, radial displacement and temperature variables are expressed by means of the relation (12) and substituted into the linearized system (cf. appendix A) leading to a trigonometric equations system. The resolution of the new system is carried out using a Galerkin approach in order to obtain a linear algebraic problem. We may note that such projection lead to a separation of the harmonics between them.

## 3 Structural equations

The structural model assumes a rigid shaft and a deformable stator. This last point constitutes the main difference with the Childs' structure modelling [17] and is expected to provide a more representative description. Fins can be located either on rotor or stator. We assume there is no eccentricity between rotor and stator axis; the displacement fluctuation is due to the flexibility of stator part. A flexural cylindrical shell theory [22] is used to describe the stator dynamic behavior.

The structure undergoes flexural vibration along the axis, as illustrated in figure 2. This kinematics is introduced adding a quadratic factor depending on axial position and clamping condition to the radial displacement.

The radial displacement is expanded as a Fourier series. Due to the fact that the spatial harmonics are uncoupled, we will consider the harmonics one by one. With such formulation, each harmonic $n$ can be directly identified to the nodal diameter of the structure.

The radial displacement for a given harmonic $n$ takes the form (with function $z$ defined in Fig. 2):

$$ w(x, \theta, t) = z(x)[a_n(t)\cos(n\theta) + b_n(t)\sin(n\theta)]. $$
<figure><img src="image_2.png" /><figcaption>Figure 2: First bending mode for high-pressure support side (left) and low-pressure support side (right)</figcaption></figure>

The stator dynamics equations are obtained calculating the Lagrangian L of the problem:

$$L = E_c - E_d + W_{ext}.$$

with, $E_c$ being the kinetic energy, $E_d$ the strain energy and $W_{ext}$ the work of external forces.

The strain energy is calculated from the Frey model [22]:

$$E_d = \frac{1}{2} \iint D\left[\left(\frac{1}{R^2}\frac{\partial^2 w}{\partial \theta^2}\right)^2 + 2\frac{\nu}{R^2}\frac{\partial^2 w}{\partial x^2}\frac{\partial^2 w}{\partial \theta^2} + \left(\frac{\partial^2 w}{\partial x^2}\right)^2 + \frac{2(1-\nu)}{R^2}\left(\frac{\partial^2 w}{\partial x \partial \theta}\right)^2\right]Rd\theta dx. \quad (13)$$

The kinetic energy is expressed as follows:

$$E_c = \frac{1}{2} \iiint \rho \dot{w}^2(x, \theta, t) r dr d\theta dx. \quad (14)$$

The fluid work on the structure is given by:

$$W_{ext} = - \int_{0}^{2\pi} \int_{0}^{L_{tot}} P(x, \theta, t) \vec{u} \cdot \vec{n} R_s d\theta dx, \quad (15)$$

where $\vec{u} = w(x, \theta, t)\vec{e}_r$ is the radial displacement, and $\vec{n}$ is the normalized surface normal vector.

As the pressure is variable within the seal, the work calculation is performed cavity per cavity, the fluid work then takes the form of a sum of contributions:

$$W_{ext} = - \sum_{i=1}^{N} \int_{0}^{2\pi} \int_{L(x_{i-1})}^{L(x_{i})} P_{i}(\theta, t) \vec{u}(x_{i}) \cdot \vec{n} R_{s} d\theta dx. \quad (16)$$

Coefficients resulting from the calculation of Lagrangian equations are given in appendix B. We define $S_1$, $S_2$, $S_{3,i}$ and $S_{4,i}$ the coefficients respectively associated with kinetic energy, strain energy, work induced by pressure fluctuation and work induced by mean pressure. This enables to write the matrix system associated to the structural equation of one cavity i:

$$\begin{pmatrix} S_1 & 0 & 0 & 0 \\ 0 & S_1 & 0 & 0 \end{pmatrix} \begin{pmatrix} \ddot{a}_n \\ \ddot{b}_n \\ \ddot{p}_{n,i}^c \\ \ddot{p}_{n,i}^s \end{pmatrix} + \begin{pmatrix} S_2 + S_{4,i} & 0 & S_{3,i} & 0 \\ 0 & S_2 + S_{4,i} & 0 & S_{3,i} \end{pmatrix} \begin{pmatrix} a_n \\ b_n \\ p_{n,i}^c \\ p_{n,i}^s \end{pmatrix} = \mathbf{0}. \quad (17)$$

Natural frequencies of a given geometry of cylindrical shell have been calculated, and compared to results of experimental modal analysis. For a large frequency range (the first 6 nodal diameters have been considered), results between numerical model and experimental data were consistent, allowing to validate the structural
model. FE simulations also confirm the correct prediction of mechanical frequencies. Moreover, based on experimental results, a correction can be added to the flexural rigidity to obtain more accurate values of mechanical frequencies for a given nodal diameter.

# 4 Strong fluid-structure coupling

To couple fluid and structural equations, one have to link the kinematics variables describing the stator motion. For each fluid cavity $i$, the radial displacement $H'_i$ (corresponding to the clearance variation) can be expressed as a function of the generalized coordinates of the cylindrical shell $a_n$ and $b_n$, depending on the clamping condition, then:

$$H'_i(\theta, t) = w(x_i, \theta, t). \qquad (18)$$

This rearrangement allows to describe the kinematics of whole problem with the variables $a_n$ and $b_n$.

Considering (N-1) cavities, the fluid system provides a set of $6(N-1)$ equations. Two additional equations are added coming from the structural problem. The global system can be written as follows:

$$[M]\ddot{\mathbf{X}} + [C]\dot{\mathbf{X}} + [K]\mathbf{X} = \mathbf{0}, \qquad (19)$$

where $[M]$, $[C]$ and $[K]$ are $[6(N-1)+2] \times [6(N-1)+2]$ banded matrix. $\mathbf{X}$ is containing kinematic and fluid variables, with the following arrangement:

$$[\mathbf{X}]^{\mathrm{T}} = [\mathbf{a}_{\mathrm{n}} \ \mathbf{b}_{\mathrm{n}} \ \mathbf{X}_{1} \ \dots \ \mathbf{X}_{i} \ \dots \ \mathbf{X}_{\mathrm{N}}], \qquad (20)$$

with

$$[\mathbf{X}_i]^{\mathrm{T}} = [\mathbf{p}_{n,i}^{\mathrm{c}} \ \mathbf{p}_{n,i}^{\mathrm{s}} \ \mathbf{u}_{n,i}^{\mathrm{c}} \ \mathbf{u}_{n,i}^{\mathrm{s}} \ \mathbf{T}_{n,i}^{\mathrm{c}} \ \mathbf{T}_{n,i}^{\mathrm{s}}],$$

where $i$ is the cavity number and $n$ the harmonic number. We will seek the complex eigenvalues ($\lambda_k$) of (19) to analyze the evolution of the pulsation (associated to the imaginary part) and the stability (associated to the real part). In order to connect the results to the classical mechanical analysis [23], we will prefer display the following quantities: $\zeta_k = \frac{Re(\lambda_k)}{|\lambda_k|}$, $f_k = \frac{1}{2\pi} \frac{Im(\lambda_k)}{\sqrt{1-\zeta_k^2}}$. By convention, a positive value of $\zeta_k$ is associated to an unstable system.

# 5 Results

## 5.1 Global view of operating strategy

Different modal analyses are carried out to highlight critical parameters on seal stability. Results of these analyses are presented in this section. Tested geometries and operational fluid parameters are representative of industrial seals. Configuration data inputs are presented in table 1. The considered solicitation mode is the 2 nodal diameter one, as it is commonly among the most subject to aeroelastic issues. Sensitivity analyses presented here give some qualitative results, to be compared and correlated with empirical knowledge of seal stability. The analysis is performed on a range of fluid inlet circumferential velocities (preswirl). Figures show evolution of fully-coupled acoustic and structural frequencies, and the corresponding evolution of reduced damping ratios.

Also, to evaluate the coupling effect, acoustic frequency and natural mechanical frequency of the cylinder in vacuum are plotted with dashed lines. The theoretical acoustic frequency is given by the formula:

$$f_{ac} = \frac{nc}{2\pi R_s}, \qquad (21)$$

where $n$ is the nodal diameter and $c$ the sound velocity.

Firstly, a one-cavity seal is studied to highlight the influence of seal-cavity length on the structural and acoustic frequencies in section 5.2. The model enables to evaluate sensitivity to different notable parameters,
<table><caption>Table 1: Data inputs of tested configurations - invariant parameters: inlet pressure $P_{in} = 8$ bar, outlet pressure $P_{out} = 1$ bar, inlet temperature $T_{in} = 293$ K, shaft radius $R_s = 76.4$ mm, stator thickness $e = 4.5$ mm, stator length $L_t = 130$ mm, nodal diameter $n = 2$, varying parameters are given below.</caption><thead><tr><th>Configuration reference</th><th>C1</th><th>C2</th><th>C3</th><th>C4</th><th>C5</th><th>C6</th><th>C7</th></tr></thead><tbody><tr><td>Support side</td><td>HP</td><td>HP</td><td>HP</td><td>HP</td><td>HP</td><td>HP</td><td>LP</td></tr><tr><td>Nominal clearance $C_r$[µm]</td><td>150</td><td>150</td><td>150</td><td>150</td><td>500</td><td>150</td><td>150</td></tr><tr><td>Number of cavities $N_{cav}$</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>3</td><td>3</td></tr><tr><td>Length $L_1$[mm] (cf. Fig. 3)</td><td>5</td><td>5</td><td>5</td><td>5</td><td>5</td><td>5</td><td>5</td></tr><tr><td>Length $L_2$[mm] (cf. Fig. 3)</td><td>105</td><td>5</td><td>50</td><td>90</td><td>90</td><td>50</td><td>20</td></tr><tr><td>Cavity length $L_{cav}$[mm]</td><td>20</td><td>120</td><td>75</td><td>35</td><td>35</td><td>75/3</td><td>105/3</td></tr><tr><td>Shaft rotation $\omega$[rpm]</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>5000</td></tr></tbody></table>

such as inlet pressure, radial clearance, number of cavities, nodal diameter, clamping condition. In section 5.3, we present the stability sensitivity to the radial clearance and number of cavities. Finally, the influence of energy equation is illustrated for a 3 cavities seal, showing our results changes when the energy equation is removed.

## 5.2 Fluid-Structure coupling sensitivity to seal cavity length

In this section, attention is paid to the influence of cavity length on the fluid-structure coupling. For a given configuration (fixed fluid parameters and cylinder geometry), the varying parameter is the cavity length (Fig. 3). The evolution of frequencies appears to be very different from one case to another (Fig. 4). For

<figure><img src="image_3.png" /><figcaption>Figure 3: Cavity length description</figcaption></figure>

configuration C1 (Fig. 4.a), the smallest cavity, there is no visible coupling between the acoustic frequency of the cavity and the structural frequency of the cylinder. For this case, the structural frequency does not appear to vary with preswirl, and its value is equal to the uncoupled natural frequency (1086 Hz). A difference between the fixed acoustic frequency (absence of preswirl) and theoretical acoustic frequency (represented with dashed lines at 1430 Hz) is observed. This point is due to the presence of leakage flow in axial direction (in opposition to a closed acoustic cavity). Conversely, in configuration C2 (Fig. 4.b), the one-cavity seal presents a very large cavity length, almost equal to the cylinder length. In this configuration, the structural frequency changes according to the preswirl, and differs from the uncoupled natural frequency. A split between these two modes is well observed. In this case, interaction between the acoustic and structural waves is notable. In configuration C3 (Fig. 4.c), the cavity length is approximately equal to half the cylinder length. The coupling between the acoustic backward wave and structural wave is notable. Indeed, from a preswirl of 70 m/s, a mode shape exchange is visible between acoustic and structural modes. This result highlights the influence of cavity geometry on fluid-structure coupling. This point is consistent with other literature studies [24].
<figure><img src="image_4.png" /><figcaption>Figure 4: Acoustic and structural frequencies evolutions for 3 different cavity lengths
(—) Forward acoustics dominated mode, (—) Backward acoustics dominated mode, (—) Forward structure dominated mode, (c) Backward structure dominated mode, (- -) Uncoupled acoustic frequency, (- -) Natural mechanical frequency.</figcaption></figure>

**5.3 Stability analysis: results**

**5.3.1 Radial clearance influence**

In this section, we will present the influence of radial clearance on configurations C4 and C5. First, focusing on frequencies in figures 5.a and 5.b, we note a coupling between structural and acoustic waves: an exchange of mode shapes is visible at around 50 m/s. Frequencies evolution between the two configurations are similar. Thus, it can be noticed that the forward acoustic wave tends to evolve more with respect to the preswirl, in the case of a radial clearance of 500 µm. This can be explained by a greater mass flow rate through the seal, allowed by a bigger clearance. It is to be noticed that the choked flow condition is verified in both cases. Damping ratios evolutions (Fig. 5.c and 5.d) present notable differences. In the case of a smaller clearance, an instability appears for a preswirl of 50 m/s whereas the other configuration remains stable until a preswirl of 140 m/s. This instability is assigned to the acoustic wave. The high sensitivity to radial clearance is in agreement and confirms the known trends of seal stability [2, 13, 24].
<figure><img src="image_5.png" /><figcaption>(a) Configuration C4 - frequencies</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>(b) Configuration C5 - frequencies</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>(c) Configuration C4 - reduced damping ratios</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>(d) Configuration C5 - reduced damping ratios</figcaption></figure>

Figure 5: Frequencies and reduced damping ratios evolutions for a 150 and a 500 µm clearance seals (—) Forward acoustics dominated mode, (—) Backward acoustics dominated mode, (—) Forward structure dominated mode, (—) Backward structure dominated mode, (--) Uncoupled acoustic frequency, (--) Natural mechanical frequency.

### 5.3.2 Number of cavities influence

To provide a more comprehensive analysis, the model is extended to the study of a 3 cavities labyrinth seal. This section presents behavior differences between a single cavity and a 3 cavities seals (configurations C3 and C6 respectively). Focusing on the coupled frequencies of the 3 cavities seal (Fig. 6.b), one identifies this time 2 additional couples of forward and backward (hereinafter respectively fwd and bwd) waves, associated with the acoustic frequencies of second and third cavities. Also, at zero preswirl velocity, a closer look reveals slightly different values of frequencies from one cavity to another, which is due to the pressure distribution along the seal. This last point also explains the difference in the evolution of acoustic frequencies with preswirl. Considering respective damping ratios evolutions (Fig. 6.c and 6.d), a stabilizing effect is clearly visible between configurations C3 (single cavity) and C6 (3 cavities). Whereas the one cavity seal appears to be unstable for values of preswirl over 30 m/s, the 3 cavities seals remains stable at any preswirl value. We may note that minimal damping ratio values for C6 are one order of magnitude larger than for configuration C3.
<figure><img src="image_9.png" /><figcaption>Figure 6: Evolutions of acoustic and structural frequencies for a single cavity and a 3 cavities labyrinth seals (—) Fwd acoustics dominated mode 1, (__) Bwd acoustics dominated mode 1, (—) Fwd structure dominated mode, (__) Bwd structure dominated mode (—) Fwd acoustics dominated mode 2, (__) Bwd acoustics dominated mode 2, (__) Fwd acoustics dominated mode 3, (__) Bwd acoustics dominated mode 3, (--) Uncoupled acoustic frequency, (--) Natural mechanical frequency.</figcaption></figure>

### 5.3.3 Energy equation influence

This section focuses on the effect of energy equation. We will compare results obtained from the whole model (including energy equation) with the one without the energy equation, for configuration (C7). This configuration is a 3-cavities labyrinth seal, considering a shaft rotation of $ω = 5000$ rpm, and LP clamping condition. Results are given in figure 7. Figures 7.a and 7.b show acoustic and structural frequencies evolutions. At zero preswirl velocity, gaps between forward and backward waves are observed due to the shaft rotation $ω$. In this configuration, energy equation has a weak influence on frequencies evolution. But, we may note a strong dependency on the stability. Focusing on the damping ratios evolutions (Fig. 7.c and 7.d), the backward structural mode appears to be unstable on a larger range of preswirl for the whole model. A positive damping ratio is observed for preswirl values of 25 to 170 m/s compared to 30 to 100 m/s for the model without energy equation.
<figure><img src="image_10.png" /><figcaption>(a) Frequencies without energy equation</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>(b) Frequencies with energy equation</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>(c) Reduced damping ratios without energy equation</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>(d) Reduced damping ratios with energy equation</figcaption></figure>

Figure 7: Frequencies and reduced damping ratios evolutions of config. C7 with and w/o energy equation (—) Fwd acoustics dominated mode 1, (—) Bwd acoustics dominated mode 1, (—) Fwd structure dominated mode, (—) Bwd structure dominated mode (—) Fwd acoustics dominated mode 2, (—) Bwd acoustics dominated mode 2, (—) Fwd acoustics dominated mode 3, (—) Bwd acoustics dominated mode 3, (—) Uncoupled acoustic frequency, (—) Natural mechanical frequency.

# 6 Conclusions

This work proposes a new analytical model, considering a full coupling between fluid behavior and stator dynamics. It enables to take temperature fluctuations into account, including energy equation at both steady and unsteady states. The model allows to conduct quick stability analyses and identify notable influence parameters. Results show that coupling phenomenon and system stability are both very sensitive to parameters like radial clearance, cavity length, number of cavities, and inlet circumferential velocity. This novel approach agrees with literature trends. The model is expected to provide better fidelity results, especially for operating conditions where temperature fluctuations have significant impact. Temperature do not have notable effect on coupled frequencies, but we observed a strong dependency of energy equation on stability results. This last point constitutes the main novelty of this new model. At the moment, a lack of experimental data cannot allow to conduct a correlation between numerical and experimental results. This point will be addressed in further investigation.
## Acknowledgments

The authors are grateful to Safran Aircraft Engines for providing the financial support of the project, and for giving permission to publish this work.

## References

[1] J. S. Alford, "Protection of labyrinth seals from flexural vibration," *Journal of Engineering for Power*, vol. 86, no. 2, pp. 141-147, 1964.

[2] D. A. Lewis, C. E. Platt, and E. B. Smith, "Aeroelastic instability in FLOO labyrinth air seals," *Journal of Aircraft*, 2012.

[3] J. S. Alford, "Protecting turbomachinery from unstable and oscillatory flows," *Journal of Engineering for Power*, vol. 89, no. 4, pp. 513-527, 1967.

[4] F. Ehrich, "Aeroelastic instability in labyrinth seals," *Journal of Engineering for Power*, vol. 90, no. 4, pp. 369-374, 1968.

[5] D. R. Abbott, "Advances in labyrinth seal aeroelastic instability prediction and prevention," *Journal of Engineering for Power*, 1980.

[6] T. Iwatsubo, "Evaluation of instability forces of labyrinth seals in turbines or compressors," *NASA CP 2250 Proceeding of a workshop at Texas A M University entitled Rotordynamic Instability Problems in High Performance Turbomachinery*, pp. 205-222, 1980.

[7] D. W. Childs and J. K. Scharrer, "An iwatsubo-based solution for labyrinth seals: Comparison to experimental results," *Journal of Engineering for Gas Turbines and Power*, vol. 108, no. 2, pp. 325-331, 1986.

[8] D. L. Rhode, S. J. Hensel, and M. J. Guidry, "Three-dimensional computations of rotordynamic force distributions in a labyrinth seal," *Tribology Transactions*, vol. 36, no. 3, pp. 461-469, 1993.

[9] T. Hirano, Z. Guo, and R. G. Kirk, "Application of computational fluid dynamics analysis for rotating machinery—part II: Labyrinth seal analysis," *Journal of Engineering for Gas Turbines and Power*, vol. 127, no. 4, pp. 820-826, 2005.

[10] T. Tsukuda, T. Hirano, C. Watson, N. R. Morgan, B. K. Weaver, and H. G. Wood, "A numerical investigation of the effect of inlet preswirl ratio on rotordynamic characteristics of labyrinth seal," *Journal of Engineering for Gas Turbines and Power*, 2018.

[11] F. Cangioli, "Effect of energy equation in one control-volume bulk-flow model for the prediction of labyrinth seal dynamic coefficients," *Mechanical Systems and Signal Processing*, 2018.

[12] F. Cangioli, "On the thermodynamic process in the bulk-flow model for the estimation of the dynamic coefficients of labyrinth seals," *Journal of Engineering for Gas Turbines and Power*, vol. 140, 2018.

[13] R. Corral and A. Vega, "Conceptual flutter analysis of labyrinth seals using analytical models. part i: Theoretical support," *Journal of Turbomachinery*, vol. 140, 2018.

[14] A. Vega and R. Corral, "Conceptual flutter analysis of labyrinth seals using analytical models. part ii: Physical interpretation," *Journal of Turbomachinery*, 2018.

[15] A. Dairien, "Stabilité des écoulements et interaction fluide-structure dans les joints labyrinthes," phd thesis, 2019.
[16] C. F. Colebrook, "Turbulent flow in pipes, with particular reference to the transition region between the smooth and rough pipe laws." *Journal of the Institution of Civil Engineers*, vol. 11, no. 4, pp. 133-156, 1939.

[17] D. Childs, "Turbomachinery Dynamics, Phenomena Modeling and Analysis.", pp.290-354, 1993.

[18] F. M. White, "Integral relations for a control volume," in *Fluid Mechanics*. McGraw-Hill, 2011, pp. 184-186.

[19] Y. Dereli and D. Eser, "Flow calculations in straight-through labyrinth seals by using moody's friction-factor model," *Mathematical & Computational Applications*, vol. 9, 2004.

[20] R. Malvano, F. Vatta, and A. Vigliani, "Rotordynamic coefficients for labyrinth gas seals: Single control volume model," *Meccanica*, vol. 36, no. 6, pp. 731-744, 2001.

[21] E.A.J. John, "Gas Dynamics". Wiley, 1979.

[22] F. Frey, "Analyse des structures et milieux continus", volume 5. Presses polytechniques et universitaires romandes, CH-1015 Lausanne, 2003.

[23] M. Géradin, D.J. Rixen, "Mechanical Vibrations, Theory and Application to Structural Dynamics", Third Edition, Wiley, 2015.

[24] T. Miura and N. Sakai, "Numerical and Experimental Studies of Labyrinth Seal Aeroelastic Instability", *Journal of Engineering for Gas Turbines and Power*, vol. 141, 2019.

<figure><img src="image_14.png" /></figure>
# Appendices

## A Nomenclature

A

$C_p$

$C_r$

c

$D$

$D_h$

$E_d, E_c$

$f_{ac}$

$H$

$h$

$L$

$\dot{m}$

$N$

$n$

$P$

$R_s$

$R$

$T$

$U$

$W_{ext}$

$u$

$\omega$

$\rho$

$\tau_r, \tau_s$

$\nu$

$\gamma$

$$ \mu_0 = \left( \frac{N}{(1-j)N+j} \right)^{1/2} $$

$$ j = 1 - \left( 1 + 16.6 \frac{C_r}{L} \right)^{-2} $$

$$ \mu_i = \frac{\pi}{\pi+2-5s_i+2s_i^2} $$

$$ s_i = -1 + \left( \frac{P_{i-1}}{P_i} \right)^{\frac{\gamma-1}{\gamma}} $$

### Sub-scripts and Super-scripts

$i$

0

′

Cross-sectional area

Heat capacity at constant pressure

Nominal radial clearance

Sound velocity

Flexural rigidity

Hydraulic diameter

Strain energy, Kinetic energy

Uncoupled acoustic frequency

Radial displacement, also radial clearance

Enthalpy per mass unit

Seal cavity length

Mass-flow rate per circumferential length unit

Number of teeth

Nodal diameter, also spatial harmonic

Pressure

Shaft radius

Specific gas constant

Temperature

Fluid circumferential velocity

Work of external forces

Internal energy per mass unit

Shaft rotational speed

Fluid density

Rotor shear stress, Stator shear stress

Poisson ratio

Heat capacity ratio

Kinetic energy carry-over coefficient

Discharge coefficient for the i-th cavity

Of the i-th cavity

Steady state

Fluctuation term
## B Linearized fluid equations

The developments of the linearization of fluid equations, using the small perturbation method.

* Continuity equation:

$$ \rho_{0i}L\frac{\partial H_i'}{\partial t} + A_{0i}\frac{\partial \rho_i'}{\partial t} + \frac{1}{R_s}[U_{0i}A_{0i}\frac{\partial \rho_i'}{\partial \theta} + A_{0i}\rho_{0i}\frac{\partial U_i'}{\partial \theta} + U_{0i}\rho_{0i}L\frac{\partial H_i'}{\partial \theta}] + \dot{m}_{i+1}' - \dot{m}_i' = 0 \quad (22) $$

* Circumferential momentum equation:

$$ \rho_{0i}U_{0i}L\frac{\partial H_i'}{\partial t} + A_{0i}U_{0i}\frac{\partial \rho_i'}{\partial t} + A_{0i}\rho_{0i}\frac{\partial U_i'}{\partial t} + \frac{1}{R_s}[2\rho_{0i}A_{0i}U_{0i}\frac{\partial U_i'}{\partial \theta} + A_{0i}U_{0i}^2\frac{\partial \rho_i'}{\partial \theta} + \rho_{0i}U_{0i}^2L\frac{\partial H_i'}{\partial \theta}] + U_i'\dot{m}_{0i+1} - U_{i-1}'\dot{m}_{0i} + U_{0i}\dot{m}_{i+1}' - U_{0i-1}\dot{m}_{i-1}' = -\frac{A_{0i}}{R_s}\frac{\partial P_i'}{\partial \theta} + \tau_{ri}'a_{ri}L - \tau_{si}'a_{si}L \quad (23) $$

* Energy equation:

$$ \rho_{0i}A_{0i}[C_p\frac{\partial T_i'}{\partial t} + U_{0i}\frac{\partial U_i'}{\partial t}] + \frac{\rho_{0i}A_{0i}U_{0i}}{R_s}[C_p\frac{\partial T_i'}{\partial \theta} + U_{0i}\frac{\partial U_i'}{\partial \theta}] + \dot{m}_{i+1}'(C_pT_{0i} + \frac{U_{0i}^2}{2}) - \dot{m}_i'(C_pT_{0i-1} + \frac{U_{0i-1}^2}{2}) + \dot{m}_{0i+1}(C_pT_i' + U_{0i}U_i') - \dot{m}_{0i}(C_pT_{i-1}' + U_{0i-1}U_{i-1}') = \tau_{r,i}'a_{ri}L + R_s\omega_{rotor} + A_{0i}\frac{\partial P_i'}{\partial t} \quad (24) $$

## C Structural equation coefficients

The calculation of kinetic energy, strain energy and fluid work respectively give:

$$ E_c = \sum_{n=2}^{N} \frac{1}{2} \rho \frac{(e^2 + 2Re)}{2} \frac{\pi L}{5} (\dot{a}_n^2 + \dot{b}_n^2) \quad (25) $$

$$ E_d = \sum_{n=2}^{N} \frac{RD\pi}{2}(a_n^2 + b_n^2)[\frac{Ln^4}{5R^4} + \frac{4\nu n^2}{3LR^2} + \frac{4}{L^3} + \frac{8(1-\nu)n^2}{3LR^2}] \quad (26) $$

$$ W = \sum_{i=1}^{N} \left( -\pi R_s \frac{L(x_i)^3 - L(x_{i-1})^3}{3L_{tot}^2} (p_{n,i}^c a_n + p_{n,i}^s b_n) - \pi \frac{L(x_i)^5 - L(x_{i-1})^5}{5L_{tot}^4} P_{0i}(a_n^2 + b_n^2) \right) \quad (27) $$