1. Introduction
The secondary air system (SAS) is essential for aero-engines moving toward a higher thermally efficient economy [
1] and operation safety [
2]. As an important part of the SAS, rotor-stator cavities provide cooling air from the compressor to different turbine components, either for cooling or sealing [
3]. Heat transfer in the rotor-stator cavity has a significant influence on the safety of gas turbine disks, which are engine-life-limited parts (ELLP) in compliance with the requirement of FAR33.70. The fracture of a gas turbine disk could directly induce non-containment of high-energy debris, resulting in catastrophic events like loss of aircraft and human life [
4]. Thus, a deep understanding of the flow structure and heat transfer characteristics in rotor-stator cavities is significant for disk cooling optimization and the reliable operation of aero-engines, which has been the subject of a large number of recent studies.
Extensive research has shown that swirl flow in the rotor-stator cavity is one of the major features that affect heat transfer and windage characteristics. Swirl flow in the rotor-stator cavity, determining the relative tangential velocity between the fluid and rotating disk surface, plays an important role in the velocity distribution of the boundary layer. The shear stress caused by the velocity gradient can generate frictional drag and heat the cooling air, which increases the power input for overcoming windage and the coolant air for necessary cooling requirements [
5]. For the swirl flow characteristics, valuable publications from different perspectives (swirl ratio, windage heating, and moment coefficient) have been conducted by experimental investigation and numerical simulation. Fernando et al. [
6] investigated the effect of surface roughness on swirl flow with the Large Eddy Simulation (LES) method, illustrating that tangential velocity profiles in the rotor-stator cavity are sensitive to even a little roughness level on the rotor surface alone. Quan et al. [
7] revealed the effect of ultrahigh rotational speed on the tangential flow and heat transfer in a millimeter-scale rotor-stator cavity. Ren et al. [
8] employed the particle image velocimetry (PIV) technique for velocity measurements in the rotor-stator cavity and showed swirl ratio distributions under different rotational Reynolds numbers and non-dimensional mass flow rates. Zhang et al. [
9] compared the axial distribution of tangential velocity under different cooling airflow rates. Zhao et al. [
10] investigated the influence of blade fracture on the distribution of the swirl ratio in the rotor-stator cavity. Tao et al. [
4] experimentally and numerically studied windage heating in a shrouded rotor-stator cavity with the superimposed flow. Moghaddam et al. [
11] conducted Computational fluid dynamics (CFD) predictions of the swirl flow in a rotor–stator cavity with a superimposed radial flow and investigated the effect on the moment coefficient and velocity distributions of different rotating bolt numbers. Kuntze et al. [
12] adopted two measuring procedures to obtain the friction torque and compared the moment coefficients for two procedures at different rotational Reynolds numbers and gap ratios. The importance of swirl flow has been widely recognized by scholars from different perspectives.
For the heat transfer characteristics, a large number of investigations focused on the heat transfer coefficient (HTC) and temperature distribution. Zhang et al. [
13] reported the effect of the mass flow rate and the rotation speed on the heat transfer of the pre-swirl rotor–stator cavity. Luo et al. [
14,
15,
16] conducted experimental research on heat transfer of rotor-stator cavities and pointed out that the dimensionless mass flow rate of cooling air is one of the key factors affecting the Nusselt number of the rotating disk. Karabay et al. [
17] investigated the effect of the swirl ratio of the cooling air on heat transfer in the rotating cavity and confirmed an optimal swirl ratio to minimize the Nusselt number. Liao et al. [
18,
19] found the pressure ratio and rotational Reynolds number exerted an influence on the magnitude of the HTC, and the local Nusselt number could be correlated with the local rotational Reynolds number by a power law. Liao et al. [
20] also numerically investigated heat transfer characteristics in the rotor-stator cavity using steady-state flow analysis, conjugate CFD, and the fully coupled finite element method (FEM)/CFD. The distributions of temperature and the HTC predicted by the three solvers were compared. Volkov [
21] conducted a conjugate heat transfer (CHT) analysis in a rotor-stator cavity, and temperature distributions of the rotor and stator surfaces were presented. Jia et al. [
22] also used the CHT method to solve for the heat transfer in a rotor-stator cavity, and the temperature rise in the turbine disk was studied. In the field of heat transfer in the rotor-stator cavity, the necessity and the applicability of the CHT method have been confirmed by many scholars.
Although a series of research has been conducted to investigate swirl flow and heat transfer in the rotor-stator cavity with significant achievements, most studies do not discuss the influence of the upstream components owing to the simplified modeling. In the typical structure of a turboshaft aero-engine, the mass flow of the cooling air in the rotor-stator cavity is controlled by the inlet labyrinth seal. Structural analysis indicated that the labyrinth seal clearance is the most susceptible component to external loads. Up to 15% of the average standard deviation is found to be due to the thermal and centrifugal load [
23]. Thus, the effect of inlet seal thermal deformation should not be neglected, especially for the high-temperature turbine disk cavity. Many scholars [
24,
25,
26,
27] have confirmed that the labyrinth clearance would experience dynamic change with combined centrifugal and thermal loads [
28], even producing non-uniform deformation along the axial direction [
29] during engine operation. Additionally, the change in the labyrinth clearance has a direct impact on the leakage flow [
30,
31], outlet swirl [
32,
33], and total temperature rise caused by windage heating [
26,
34], which all have a non-negligible influence on the heat transfer in downstream rotor-stator cavities. However, there is still a lack of in-depth and comprehensive evaluation of the effect of inlet seal thermal deformation on the heat transfer in rotor-stator cavities.
The contribution and novelty of the current research focus on two primary features: (1) the non-uniform hot-running profile of the inlet labyrinth was found, and its effect on swirl flow and heat transfer in the rotor-stator cavity is first discussed; (2) a clearance compensation method is proposed and validated to weaken the adverse effects caused by the non-uniformity thermal deformation of the seal clearance.
In this paper, a numerical framework was established by integrating CHT analysis and structural FEM analysis, considering two-way aero-thermo-elasto coupling interaction among elastic deformation (affecting seal clearance), leakage flow (affecting HTC of the disk), and heat transfer (affecting solid temperature). A comparison of the flow structure, swirl flow, and heat transfer under the cold-built clearance (CC) and the actual non-uniform hot-running clearance (ANHC) was conducted. To highlight the influence of non-uniformity of the labyrinth clearance, the nominal uniform hot-running clearance (NUHC) was defined as the arithmetic mean of tip clearances and also compared against the ANHC. Finally, a clearance compensation method was proposed and validated to weaken the adverse effects caused by the non-uniformity of the seal clearance. The present investigation establishes the foundation for subsequent engineering applications for the safety design of the aero-engine SAS.
2. Methodology and Modeling
2.1. System Configuration and the Simplified Model
This study aimed at the swirl flow and heat transfer characteristics of the rotor-stator cavity in consideration of the thermal deformation of the inlet labyrinth seal. In the typical structure of turboshaft engines, as shown in
Figure 1 [
35], a labyrinth is installed between the outlet of the centrifugal compressor disk cavity and the inlet of the turbine disk cavity. The labyrinth is situated below the circumfluence combustor, metering the mass flow rate of the cooling air from the compressor outlet to the turbine disk cavity, which is used to cool the turbine disk and prevent gas from invading the disk cavity. However, the labyrinth seal at the inlet of the rotor-stator cavity experiences harsh thermal loads, which would lead to thermal deformation of the labyrinth seal and result in heat transfer degradation in the rotor-stator cavity.
To facilitate this study, the structure of the turbine disk cavity and inlet labyrinth seal was simplified as follows:
- (a)
The influence of the compressor diffuser and the turbine guide vane was ignored, and only the thin-wall construction between the compressor disk and the turbine disk was analyzed, which is the main factor in the thermal deformation of the inlet labyrinth seal.
- (b)
Gas turbine blades were ignored, and an equivalent centrifugal load was adopted to the turbine disk rim.
- (c)
The supporting function of the coupling sleeve teeth and bearings was simplified to a simple displacement constraint.
The simplified computational domain of the rotor-stator cavity and inlet labyrinth seal is illustrated in
Figure 2. As can be seen in
Figure 2, the green part, the orange part, and the gray part are the stationary assembly, the rotating assembly, and the fluid domain, respectively. The main geometric parameters are specified in
Table 1. As a general design experience, it is recommended that the seal clearance be 0.15–0.3% of its radius [
36]. In the current study, the labyrinth was located at the inlet of the turbine disk cavity, which is a hot component, so the upper value of the recommended seal clearance (0.3%) was adopted.
2.2. Governing Equations
The governing equations included Navier-Stokes equations for the fluid domain and thermoelastic equations for the solid domain.
For the fluid domain, Navier-Stokes equations are shown as Equations (1)–(3).
In the above equations, V is the velocity vector, f is the body force, Q is the heating rate of internal heat source, and is the dissipation function.
For the solid domain, thermoelastic equations are shown as Equations (4)–(6).
In the above equations, J is the Jacobi determinant of coordinate transformation, u is the displacement vector, e = u is the volume strain, is Poisson’s ratio, is thermal expansion coefficient, and , are Lame constants.
For fluid-solid interface, force and energy balance equations are shown as Equations (7)–(9).
and are position vectors at the solid boundary and the fluid boundary respectively, is the external normal direction of the fluid-solid interface, and is the force vector of the fluid on the solid.
In this investigation, the body force of the air was ignored, and there was no internal heat source in the computing domain. The steady-state simulation (time-dependent analysis) was carried out numerically, so the time derivative terms in the above equations are all zero. In Equation (6), is the coupling term of the thermoelastic equation, reflecting the temperature change caused by elastic deformation. The thermoelastic coupling term of the turbine disk is a tiny amount compared to the heat transfer term, so the thermoelastic coupling effect was ignored in the solution process.
2.3. Numerical Method
For the rotor-stator cavity and inlet labyrinth seal, there is a complex multi-physics coupling relationship among the labyrinth leakage flow, heat transfer in the disk cavity, solid conduction, and structural deformation, as shown in
Figure 3. The labyrinth clearance is the flow boundary of leakage flow, which determines the mass flow rate, total temperature, and swirl ratio of cooling air in the disk cavity, thereby affecting the convective heat transfer of the turbine disk. The distribution of the solid temperature field is closely related to the heat transfer of the turbine disk, which will cause thermo-elastic deformation of the structure, and in turn, affect the labyrinth clearance. In order to accurately solve the two-way aero-thermo-elasto coupling relationship among leakage flow, convective heat transfer, solid heat conduction, and structural deformation, a computational strategy was developed by integrating conjugate heat transfer (CHT) analysis and structural finite element method (FEM) analysis under the platform of ANSYS Workbench (a commercial software), as shown in
Figure 4. The CFX module simultaneously solves fluid convective heat transfer and solid heat conduction by finite volume method (FVM), realizing CHT analysis. The fluid field in the rotor-stator cavity and the temperature field of the rotor and the stator were evaluated by considering fluid-solid interaction. Moreover, the Static Structural module solves the structural elastic deformation by FEM, realizing structural deformation analysis. The radial displacements of the rotor and stator and, consequently, the new seal clearance profile were evaluated by considering the temperature distribution and the centrifugal load. The simulation was a non-linear analysis because the physical properties of the fluid and the solid are temperature dependent. The fluid was regarded as ideal air, which is a compressible fluid. The dynamic viscosity and thermal conductivity of the fluid were subject to the Sutherland formula in
Appendix A. The effects of temperature on the thermal and mechanical properties of solid materials are considered in
Appendix B.
The solution process is shown in
Figure 5. Solid temperature and labyrinth clearance were transferred between the two modules, and the solution was iteratively solved until the solution converged in terms of seal clearance. The seal clearance convergence criterion was the clearance difference between the two adjacent iterations and shall not exceed ±5 μm. Firstly, the labyrinth clearance was estimated, and the CHT analysis was carried out under the initial clearance. Then the solid temperature field was transferred to the FEM solver for structural deformation analysis. Finally, the calculated seal clearance was compared with the initial clearance. If the convergence criterion was satisfied, the calculation stopped; otherwise, the geometry and mesh were updated according to the new seal clearance, and then the next iteration cycle began.
Figure 6 shows the computational grids of the solid domain and fluid domain. Hexahedral grids were generated for the fluid domain, while hexahedral and triangular prism hybrid grids were adopted for the solid domain. For the sake of completely taking into account the inlet seal thermal deformation effect, the mesh was refined specifically for the labyrinth seal to capture the temperature gradient and radial displacement variation of the labyrinth seal. Mesh refinement near the wall was carried out to ensure that y
+ was approximately 1, which was suitable for the Shear Stress Transport (SST) turbulence model. Taking into account axis symmetry, a 10° periodic sector model (1/36 sector) was selected in both the CHT and FEM analyses to reduce the time cost of the aero-thermo-elasto coupling calculation. In the CHT analysis session, a rotational periodicity interface model was performed for the periodic symmetry planes. The upwind scheme was used to discretize the turbulence terms, and the convection term was discretized by the high-resolution scheme. The air was regarded as a compressible fluid, and the total energy method was selected to simulate heat transfer. In the FEM analysis, the periodical symmetrical constraints were assigned for the circumferential symmetry planes. Other constraints were set in the specific location, as depicted in
Figure 2. A specific angular velocity,
Ω, was loaded to the rotating disk, which was defined by the rotational Reynolds number
. The materials of the rotor and stator were two kinds of high-temperature Ni-based superalloy.
Appendix B shows the thermal expansion coefficient, elastic modulus, and Poisson’s ratio of the solid materials in detail.
2.4. Dimensionless Parameters and Boundary Conditions
In order to improve the generality of the results, the boundary parameters and main operation parameters involved in this paper were subjected to dimensionless processing as follows. Radial structural deformation (
) and labyrinth clearance (
c) were normalized by the cold-built labyrinth clearance (
):
As for the fluid flow, the mass flow rate of through-flow cooling air (
) was expressed by the flow coefficient, which is defined as:
where
R is the gas constant, 287 J/(kg·K) for the air, and
is the area of the cold-built seal clearance.
Pressure distribution is related to flow resistance, which is an important flow characteristic of the rotor-stator cavity. The pressure coefficient was adopted to describe the pressure distribution of the rotating disk surface, which is defined as:
where
r is the radius.
Swirl ratio is the ratio of air swirl velocity (
) to local disk rotating speed, reflecting the relative circumferential velocity between cooling air and rotating disk.
The moment coefficient is the dimensionless friction moment of the rotating disk.
where
is the friction moment of the rotating disk.
For the aspect of heat transfer, the Nusselt number represents the thermal resistance ratio of heat conduction and convective heat transfer in the fluid layer near the disk surface, reflecting the intensity of convective heat transfer.
where
is the convective heat transfer coefficient,
is the thermal conductivity of the fluid,
is the heat flux of the rotating disk surface, and
is the temperature of the rotating disk surface.
The solid temperature of the rotating disk (
T) is expressed by dimensionless metal temperature.
In the present research, the cruise state, maximum continuous state, and take-off state were selected as the typical operating states of the turboshaft engine. The aero-thermo-elasto coupling analysis was conducted to explore the swirl flow and heat transfer characteristics in the rotor-stator cavity in consideration of the inlet seal thermal deformation effect. The centrifugal load of the turbine blade,
, was equivalent to the uniform load on the turbine disk rim. The boundaries ➀–➆ in
Figure 2 were all given the boundary conditions of convection heat transfer. The reference temperatures of convective heat transfer were the fluid temperatures near the wall, which could be gained from the SAS flow network solver [
2]. Convective heat transfer coefficients were derived from the correction formulas, as shown in
Appendix C. The boundary parameters of these three typical operating points are shown in
Table 2.
2.5. Model Verification
A grid sensitivity analysis was also performed to define the mesh size. Three sets of grids with different element numbers were established. The element numbers of coarse, medium, and refined meshes were 0.31 million, 0.62 million, and 1.20 million respectively. Computational convergence was accepted when the residuals of the continuity equation, momentum equation, and energy equation converge to the order of 2 × 10
−6.
Figure 7 shows the dimensionless metal temperature on the surface of the turbine disk under cruise state, in which the differences among coarse, medium, and fine grids are no more than 2.5%. The dimensionless metal temperature results can be considered independent from the mesh size. Considering the compromise between simulation accuracy and calculation cost, the medium grids were selected for subsequent calculation and analysis.
The analysis of convective heat transfer on the surface of the rotating disk is one of the most critical links in the aero-thermo-elasto coupling calculation. In order to verify the accuracy of the numerical simulation method in this paper, experimental data [
37] of heat transfer in an air-cooled rotor-stator cavity were used to implement the reliability test. The rotating-disk rig comprised two disks of 762 mm diameter. The axial clearance between the two shrouds was less than 4 mm, and the axial gap ratio for the disks was 0.12. A radial outflow of air could be supplied through the center of the stationary disk (100 mm diameter). The back face of the rotating disk, which was made from steel, was heated by 30 stationary, radiant-heater elements located at an axial distance of approximately 30 mm from the disk. Ten fluxmeters were located at different radial positions, and ten copper-constantan thermocouples were embedded at the same radial locations but displaced tangentially from the fluxmeters. The fluxmeters were calibrated using free-disk tests. The swirl ratio in the cavity and the Nusselt numbers on the surface of the rotating disk were measured in the experiment.
The numerical model was established, as shown in
Figure 8. The simulation results were compared with the experiment data in
Figure 9 and
Figure 10. Three turbulent models (k-ε, k-ω, and SST) were examined in three sets of experiment cases (
Cw = 2530,
Cw = 6100, and
Cw = 9680, in which
Cw =
). Swirl ratios from different turbulent models were validated by experimental data in
Figure 9. The SST model performed best, especially at the transition between the core region and near-wall region. Nusselt number distributions from all three turbulent models presented a similar tendency with the experiment data and are shown in
Figure 10. The calculation results indicate that the numerical method in the present work can simulate the heat transfer in the rotor-stator cavity with reasonable accuracy. Generally, the SST model had a better performance for swirl flow and heat transfer simulation in the rotor-stator cavity. The SST model modified the definition of the eddy viscosity coefficient and considered the transmission effect of turbulent shear stress. Therefore, the SST model has higher accuracy and credibility in a wide range of flow fields, especially in near-wall free flow. Considering the above advantages, the SST model was used in the subsequent numerical simulation.