CFD Analysis of Pressure
Instabilities in Stator-Rotor
Disc Cavity Systems
Pedro Santiago Parras Bl´azquez
Link¨opings universitet Institutionen f¨or ekonomisk och industriell utveckling
Link¨opings universitet Institutionen f¨or ekonomisk och industriell utveckling ¨
Amnesomr˚adet Mekanisk v¨armeteori och str¨omningsl¨ara Examensarbete 2017|LIU-IEI-TEK-A–17/02936-SE
CFD Analysis of Pressure
Instabilities in Stator-Rotor
Disc Cavity Systems
Pedro Santiago Parras Bl´azquez
Academic supervisors: Magnus Andersson & J¨org Schminder Industrial supervisor: Niklas Edin
Examiner: Hossein Nadali Najafabadi
Link¨oping universitet SE-581 83 Link¨oping, Sverige
Abstract
The continuous demand to improve turbine performance has led manufactures to focus on aspects that have been previously considered of secondary importance such as the secondary air system. The purpose of this system is to cool the components and prevent ingestion of hot gas into the stator-rotor cavity that could lead to low frequency pressure fluctuations called Cavity Flow Instabilities. These instabilities could cause unpredictable rotor vibrations and damage several components.
A CFD method capable of detecting cavity flow instabilities in a rotor-stator disk cavity system is investigated, based on the 360◦ model of the cavity without any stator vanes and rotor blades. Boundary conditions are simplified by considering steady and uniform flow in the main gas path. Different turbulence models are tested such as Realizable k − ε, k − ω SST , DDES, and SAS. In order to test the performance of the method, different purge flow levels are simulated.
The most successful results, are predicted by the Realizable k − ε turbulence model. This model predicts two rotating low pressure structures in the cavity, for low purge flow levels. These pressure structures rotate at approximately 80% of the rotor speed. Furthermore, the spectra analysis of the pressure shows a reasonable agreement with the experimental results in terms of the frequency, showing a dis-tinct region of low frequencies pressure instabilities. Nonetheless, this method over predicts the amplitudes by a factor of 3-7 depending on the frequency. In addition, regions of one order of magnitude higher in frequency is also predicted. The DDES model shows similar findings but the amplitudes in the pressure spectra associated to the low frequencies are lower. Additionally, SAS also predicts the pressure in-stabilities but, in this case, the amplitudes are closer to RANS simulations, yet the high frequencies disappear. Unfortunately, k −ω SST , did not predict these pressure instabilities.
Further research is still needed in many of the aspects of this work, from the simplifications up to the turbulence models. However, it is concluded from this work that this method could be a useful tool for turbine design as it decreases the need for testing and prototype manufacturing.
Acknowledgements
This Master Thesis was carried out at GKN Aerospace in Trolh¨attan between Jan-uary and June 2017 under supervision of Link¨oping University and I would like to thank all the people involved for giving me the opportunity to work in such an interesting field.
I would like to express my deepest gratitude to my supervisor from GKN Aerospace Niklas Edin for his guidance and constant support providing me with useful and practical advise throughout the course of this thesis.
My special thanks to Hans Abrahamsson for his insightful comments, recom-mendations and ideas during the project meetings as well as his feedback.
I would also like to thank Pieter Groth and Staffan Brodin for their motivation and fruitful discussions that we had during the meetings.
Another person that I would like to thank is Lars Ljungkrona for his help and tips during the meshing process as well as his recommendations throughout the course of the thesis.
Last but not least, thanks to my academic supervisors at Link¨oping University, Magnus Andersson and J¨org Schminder as well as my examiner Hossein Nadali Najafabadi for their support, advise and feedback on the work.
Nomenclature
Abbreviations and Acronyms
Abbreviation Meaning
CFD Computational fluid dynamics
DDES Delayed Detached Eddy Simulation
DES Detached Eddy Simulation
DNS Direct Numerical Simulation
LES Large Eddy Simulation
RANS Reynolds-Averaged Navier-Stokes Equations
SAS Scale Adaptive Simulation
SBES Stress Blended Eddy Simulation
SRS Scale Resolving Simulation
SST Shear Stress Transport
URANS Unsteady Reynolds-Averaged Navier-Stokes Equations
CI Combined Ingress
EI Externally-induced Ingress
RI Rotationally-induced Ingress
ICAS-GT Internal Cooling Air Systems of Gas Turbines
Latin Symbols
Symbol Description Units
a1 k − ω SST model constant [−]
b Rim seal radius [m]
cr Radial component of the velocity [ms−1]
cx Axial component of the velocity [ms−1]
cθ Tangential component of the velocity [ms−1]
CDES Constant of the DES model [−]
Cp Pressure coefficient [−]
CW Non-dimensional purge flow [−]
Cµ k − ε model constant [−] C1ε k − ε model constant [−] C1ε k − ε model constant [−] Cω1 k − ω SST model constant [−] Cω2 k − ω SST model constant [−] DH Hydraulic diameter [m]
E Total energy per unit mass [m2s−2]
FDDES Shielding function for the DDES model [−]
f Body force per unit mass [ms−2]
f Frequency [Hz]
Symbol Description Units
f1 Blending function in the k − ω SST
model
[−]
I Turbulent intensity [−]
k Turbulence kinetic energy per unit
mass
[m2s−2]
kT Thermal conductivity [W m−1K−1]
Lt Turbulence length scale [m]
Lvk von K´arm´an length scale [m]
˙
mcool Purge flow mass flow [kgs−1]
P Pressure [P a]
Pref Reference pressure [P a]
ReDH Reynolds number based on the
hy-draulic diameter
[−]
Reφ Rotational Reynolds number [−]
Rex Axial Reynolds number [−]
S Strain tensor rate [s−2]
Sφ Scalar source term [kgm−3s−1]
t Time [s]
T Temperature [K]
u Velocity [ms−1]
u∗ Friction velocity [ms−1]
u+ Non dimensional velocity [−]
U0 First derivative of the velocity [s−1]
U00 Second derivative of the velocity [m−1s−1]
v Velocity [ms−1]
x Coordinate [ms−1]
y First cell height [m]
y+ Non-dimensional wall distance [−]
Yk Dissipation term [kgm−3s−1]
Greek Symbols
Symbol Description Units
βR Radial velocity ratio [−]
βT Radial velocity ratio [−]
β1 k − ω SST model constant [−]
β2 k − ω SST model constant [−]
β∗ k − ω SST model constant [−]
Γφ Scalar diffusion coefficient [m2s−1]
δ Kronecker delta [−]
∆ Cell edge length [m]
ε Dissipation rate [s]
Symbol Description Units
ζ2 SAS model constant [−]
ζ3 SAS model constant [−]
κ k − ω SST model constant [−]
µL Dynamic viscosity [kgm−1s−1]
µT Eddy viscosity [kgm−1s−1]
ν Kinematic viscosity [m−2s−1]
ρ Fluid density [kgm−3]
σk k − ε and SAS models constant [−]
σk1 k − ω SST model constant [−]
σk2 k − ω SST model constant [−]
σε k − ε model constant [−]
σφ SAS model constant [−]
σomega1 k − ω SST model constant [−]
σomega2 k − ω SST model constant [−]
τ Viscous stress tensor [s−1]
τF Favre-averaged Reynolds-stress tensor [s−1]
τw Wall shear stress [P a]
φ Variable in the SAS scale equation [s−1]
φ Scalar representing a flow parameter [−]
ω Specific dissipation rate [s−1]
Contents
1 Introduction 1 1.1 Background . . . 1 1.2 Aim . . . 3 1.3 Literature Survey . . . 3 1.4 Limitations . . . 6 2 Theory 7 2.1 Computational Fluid Dynamics . . . 72.2 Governing Equations of Fluid Flow . . . 8
2.3 Turbulence Modelling . . . 9
2.3.1 Reynolds-Averaged Navier-Stokes Equations and Classical Tur-bulence Modelling . . . 10
2.3.2 Scale Resolving Simulation . . . 13
2.3.3 Wall Treatment . . . 15
3 Method 17 3.1 Geometry and Mesh . . . 17
3.2 Boundary Conditions . . . 19 3.3 Solver Settings . . . 21 3.4 Convergence Criterion . . . 23 3.5 Evaluated Cases . . . 24 3.6 Post-Processing . . . 24 4 Results 27 4.1 Realizable k-ε model . . . 27 4.1.1 Unsteady Results . . . 27 4.1.2 Time-averaged Results . . . 31
4.1.3 Unsteady Pressure Analysis . . . 34
4.2 k-ω SST model . . . 36
4.2.1 Unsteady Results . . . 36
4.2.2 Time-averaged Results . . . 37
4.2.3 Unsteady Pressure Analysis . . . 37
4.3 Scale Resolving Simulation . . . 38
4.3.1 Unsteady Results . . . 38
4.3.2 Time-averaged Results . . . 39
4.3.3 Unsteady Pressure Analysis . . . 40
5 Discussion 43 5.1 Method . . . 43
5.2 Results . . . 45
5.2.1 URANS . . . 45
5.2.2 Scale Resolving Simulation . . . 47
7 Outlook 51
8 Perspectives 53
Appendices 57
A First appendix 57
A.1 Realizable k-ε Further Results . . . 57
A.1.1 Instantaneus Contours . . . 58
A.2 Experimental Unsteady Pressure Analysis . . . 61
B Second appendix 63 B.1 SRS Mesh Refinement . . . 63 B.1.1 Method . . . 63 B.1.2 Results . . . 63 B.1.3 Discussion . . . 66 C Third appendix 67 C.1 Stress Blended Eddy Simulation . . . 67
C.1.1 Theory . . . 67
C.1.2 Method . . . 67
1
Introduction
1.1
Background
The turbomachinery industry is faced with a continuous demand to improve engine performance, meet strict environmental and safety regulations as well as reducing development time and cost of new products. Considerable advances have been made during the last 50 years in order to achieve those objectives. However, improvements have become more difficult to achieve, for instance in component efficiencies. In ad-dition to this, development in other areas has become expensive. Thus, the scope has been shifted towards other aspects that had been previously considered of secondary importance.
One of the paths that have been followed to improve the performance of engines is the increase of turbine inlet temperature. The higher the temperature of the gas entering the turbine the greater the efficiency of the engine. However, this leads to gas temperatures that are higher than the critical yield limit of turbine materials. Also, the possible ingestion of hot gas into the turbine disc cavities or wheel-space (Fig. 1) could cause overheating of the discs. It is therefore necessary to have a secondary air system that is able to cool the blades and prevent ingestion to limit the deterioration of these components.
Fig. 1: Typical high-pressure turbine stage showing rim seal and wheel-space [1]
The analysis and study of the secondary air system is crucial in the life ex-pectancy and performance of turbines and engines. The air mass flow is usually purged from previous stations, such as the compressor, which causes a decrease in overall performance of the engine. Therefore it is essential to reduce it to a minimum while still preventing the ingestion of hot air into the cavity.
As mentioned earlier the function of the secondary air system is to prevent ingestion into the cavity which can cause not only material deterioration due to high temperatures, but also low frequency pressure fluctuations that could cause unpredictable rotor vibrations. These pressure fluctuations are called cavity flow instabilities (CFI) and they play a major role in turbine design specially in space
turbines due to its high loading. Hence it is important to understand the phenomena of hot gas ingestion.
Ingress from the gas channel into the cavity can occur due to several reasons, for instance as a consequence of the nonaxisymmetric pressure variation when the flow passes over the vanes and blades [2]. The interaction of the vanes and blades generates high-pressure regions where the static pressure of the hot gas exceeds the cavity pressure (Fig. 2). This type of ingestion is called Externally-induced ingress (EI ingress).
Fig. 2: Static pressure distribution in a turbine annulus [2]
Another type of ingress that still can occur, even if the external flow were ax-isymmetric, is called Rotationally-induced ingress (RI ingress). The fluid in the cavity while rotating creates a radial pressure gradient which can cause the pressure to be lower than in the annulus. This causes the egress of fluid near the the rotor where the centrifugal effects are strongest (disk pumping effect [2]).
When the effects of EI ingress and RI ingress are significant the term Combined ingress (CI ingress) is used. This type of ingestion plays an important role when the engine is operating ’off-design’ conditions [2].
In order to reduce the amount of hot-gas ingestion into the cavity rim seals are usually fitted at the periphery of the wheel-space. Typical rim seal geometries are shown in Fig. 3.
Fig. 3: Different rim seal configurations: (a) axial-clearance, (b) radial-clearance, (c) double radial-clearance, and (d) engine-representative double-clearance seal [2]
The growing trend in the industry is to use complex 3D unsteady computational fluid dynamics code to study the mechanisms of ingress, its consequences and how to
avoid it. However, as with any fluid simulation, experimental research is necessary to validate these codes. Even though there is some understanding of the basic parameters influencing the CFI phenomena, there is no method capable of correctly predicting them. Several attempts have been made to replicate test conditions where instabilities were found with partial success [3] [4].
1.2
Aim
The main objective of this thesis work is to develop a computational method that is able to predict the cavity flow instabilities as measured by experimental analysis. This computational method should be applicable for different operating conditions and stator-rotor cavity geometries in order for it to be used as an engineering and design tool for future turbine development.
1.3
Literature Survey
During the past years several theoretical, experimental and numerical studies have been conducted in order to better understand the low-frequency pressure fluctuations that occur in the stator-rotor disc cavity. These studies included different cavity ge-ometries, rim seals and operating conditions so that a complete understanding of the phenomena and the parameters affecting could be obtained.
ICAS-GT1
The European project Fluid Flow and Heat Transfer within the Rotating Internal Cooling Air Systems of Gas Turbines (ICAS-GT) was developed in order to improve the performance of gas turbines as well as reducing the design cycle costs and time. The approach adopted in this project included advanced experimental work as well as an evaluation of flow and heat transfer models.
The experimental investigation of the hot gas ingestion phenomena was carried out by Aachen University in a 1.5 stage turbine test rig (Fig. 4). Steady and un-steady pressure as well as carbon dioxide concentration measurements were done. These measurements were taken for several different configurations including a base-line and alternative rim seal geometries in the front and back cavity for five different operating points.
Smout et al. [3] found that the hub pressure distribution in the main gas path channel upstream of the rim was influenced by the amount of cooling air ejected into the hot gas path. Moreover, a periodic flow structure inside the front cavity was detected. The frequency of these pressure fluctuations was not coupled to the rotor frequency as its frequency was an order of magnitude lower than the blade passing frequency.
All measurement techniques implemented, detected the existence of a transition point, where the low-frequency pressure fluctuations in the cavity disappeared. This transition point was detected by the authors for all points of operation depending on the magnitude of sealing air mass flow.
Fig. 4: 1.5 Stage Turbine Rim Sealing Rig used by Aachen University for the experimental analysis during ICAS-GT1 and ICAS-GT2 [3]
As part of the numerical part of the project 2D unsteady simulations were con-ducted on the main air flow. Steady 3D simulations were also done for a simplified 22.5◦ sector of the cavity as well as for a stator-blade sector including all stages and cavities. For this last model, an unsteady simulation was also performed.
The authors [3] concluded that the rotation of the core in the cavity was not modelled correctly. They could confirm, however that a transition point existed as measured in the experiments. They also detected that their computations under-predicted the radial pressure gradient therefore over-predicting the sealing efficiency. It was believed that the reason was a poor computation of the static pressure field. Finally, due to the use of periodic conditions and only simulating one stator-blade sector they concluded that it was not possible to capture the unsteady effects with longer wave-lengths than the blade-to-blade distance. Thus, no low-frequency pres-sure fluctuations were found.
ICAS-GT2
When the ICAS-GT1 program was finished, many unanswered questions still re-mained with regard to the accuracy of the numerical work. During ICAS-GT1 there were considerable differences between the experimental results and CFD calcula-tions. These differences could not be explained at the time. Therefore, the first phase of ICAS-GT2 numerical work (continuation of ICAS-GT1) was focused on the resolution of these differences.
The industrial partners of the ICAS-GT2 program conducted a detailed numer-ical study in order to identify appropriate calculation methods for the ingestion of hot gas into rotor stator cavities. Jakoby et al. [4] used three different types of com-putational domains. The first one consisted of a 22.5◦ sector of the full geometry with the hot gas path with nozzle guide vanes, the rotor blades and the outlet guide vanes as well as both cavities. The second studied domain included only a 22.5◦ sector of the front cavity. Finally, the last grid consisted of the front cavity, but this time the full circumference, 360◦, was used in order to fully capture those unsteady pressure fluctuations.
Fig. 5: Static cavity pressure distribution at mid-plane for CW = 7000, showing three distinguishable low pressure regions [4]
The authors detected a large scale rotating structure, that revolved with 80% of the rotor speed in the front cavity (Fig. 5). Furthermore, these low-frequency fluc-tuations were not induced by the vane and blade interaction since their frequency was an order of magnitude lower than the blade-passing frequency, similar results to what Aachen University found in their measurements [3]. They concluded that this rotating structure significantly influenced the ingestion of hot gas into the turbine wheel space.
Other Work
Cao et al. [5] modelled the flow inside an unshrouded disc cavity without con-sidering the vanes and blades. They used the full 360◦model as well as a 90◦ sector. The geometry considered had a simple axial gap between rotor and stator disks with a uniform cavity. Their study showed unsteady pressure structures that rotated at speeds of 90-97% the rotor speed. The authors believed that the large-scale struc-tures were associated with the interaction between the main annulus and cavity flow. Moreover, they found that as the gap width was reduced, keeping the same net rim seal flow, the ingestion was suppressed and the unsteadiness weakened.
Julien et al. [6] simulated 74◦ sector including the disk cavity, vanes and blades in order to study the dynamic influence of vane-blade interaction on cavity flows for: no purge flow, low purge and high purge flow rates. They detected large scale flow structures for zero and low purge flow rates. For these two cases, the most energetic frequencies were lower than the blade frequency and concluded that they were associated with those large scale structures rotating at a lower angular speed than the rotor. For the high purge case, the authors found two competing frequencies one of them associated with the blade frequency.
Another numerical study performed by Wang et al. [7] included the complete annulus with all vanes and blades as well as a double seal forming an outer rim and an inner disc cavity. The results showed unsteady pressure structures, rotating in the same direction and at 86% of the rotor speed. These structures were detected for low and medium purge flow rate. For low purge flow rate however, the cells were found to be slightly disturbed by the ingress and egress of the flow.
O’Mahoney et al. [8] performed URANS and LES for a 13.3◦ model including guide vanes and rotor blades. The authors found that the pressure profiles in the cavity were not in good agreement with experiments. However, they detected the presence of low-frequency pressure variations in the URANS and LES results, with smaller amplitudes for the latter. They also performed URANS simulations on larger sector models, showing small changes in the analysis of the unsteady pressure varia-tions. The authors concluded that the sector size is not as important as turbulence modelling in the prediction of rim seal flows.
Finally Arzpeima [9] performed a preliminary study of the flow inside the disc cavity of a test turbine. The author modelled the main annulus stage separately in order to obtain the boundary conditions for the cavity. The annulus included one stator and rotor blade passage and the cavity was the full 360◦ model. No unsteady rotating pressure structures were observed for that configuration and op-erating points.
1.4
Limitations
The software used through this thesis are ANSYS R ICEM CFD v18.0 for meshing, ANSYS R Fluent v18.0 to set up and run the simulations and ANSYS R CFD-Post v17.1 for post processing the results.
In the beginning ANSYS R CFX v17.1 was the chosen software to set up and run all the cases. However, after running the first cases with unsuccessful results the switch to Fluent was done as it included the Realizable k − ε turbulence model that had been previously used with successful results [4], which CFX does not include.
In order to use these programs, 14 licenses are available to share among the engineering department. No license problems have been encountered, except in rare occasions.
Simulations are run in a Linux cluster with a total of 1024 nodes. A maximum of 40 nodes have been used for each simulation. With this amount of computer power, unsteady RANS simulations took approximately 6 hours to complete 1 rotor revolution.
Furthermore, additional post processing has been done with MATLAB R 2015a. There are 14 licenses also shared among the engineering department, even though no license problem has been encountered.
The experimental and numerical data from the mentioned projects ICAS-GT1 and ICAS-GT2 is not as complete as expected and some problems have been en-countered reading some of the experimental data files. Moreover, there is limited information regarding measurement techniques, thus no conclusions can be drawn regarding the uncertanties of experiments.
In addition to this, the numerical data from the industrial partners is also lim-ited. Even though detailed reports of the numerical analysis were available from Volvo Aero, the most successful analysis is not documented in detail. For instance solver settings and boundary conditions are not fully detailed and therefore no direct comparison could be made.
2
Theory
2.1
Computational Fluid Dynamics
Computational fluid dynamics (CFD) is the numerical analysis of systems involving fluid flow as well as heat transfer and chemical reactions. The equations describing the behaviour of the fluids are discretised and solved in a domain which is divided into smaller sub-domains, known as a mesh or grid.
CFD has numerous advantages, for instance, it presents the perfect opportunity to study specific terms in the governing equations in detail [10]. In addition to this, it provides considerable reduction of lead times and cost in new designs. CFD has the capability of simulating conditions that are not possible or difficult to achieve in experimental facilities such as natural disasters. Finally, CFD complements experi-mental an analytical approaches by providing an effective means of simulating real fluid flows where high-quality experimental facilities can be more expensive.
However, all this advantages do not suggest that CFD will replace experimental testing for design purposes in the near future. Numerical errors exist in computations and differences between simulations and reality exist. Furthermore, it is not a fully mature technique, there is still much to be investigated regarding many aspects, such as turbulence or discretisation techniques. Finally, deep knowledge and a careful analysis of the case to study is needed as there is not one method that is universally applicable.
CFD codes are structured around the numerical algorithms that can tackle fluid flow problems. All codes contain three main elements or modules: a pre-processor, a solver and a post-processor [11].
1. Pre-processor: The first step in any CFD analysis is the definition and creation of the geometry of the flow region. The second step, is the mesh generation, which is one of the most important steps in the pre-process stage. The domain is subdivided into smaller subdomains where the flow physics is solved. The accuracy of a CFD solution is strongly influenced by the number of cells and the type of the mesh. The following step is the selection of physics and fluid properties. The last step in this module is to specify the boundary conditions which should mimic the real physical representation of the fluid flow [10].
2. Solver: This module solves the flow equations in the domain. In order to solve the governing equations, they must be first discretised. There exist three main numerical discretisation techniques: finite element, finite difference and spectral methods. The main technique used in commercial CFD solvers is finite volume. The equations are solved in an iterative manner, thus the solution needs an initial guess which is crucial to the iterative procedure. If the initial guess is close to the final solution, the iterative process will converge faster. To asses whether the solution has been successful (converged) the imbalances in the conservation of flow properties can be tracked.
3. Post-processor: In this last step, the results are analysed and reviewed. Most modern commercial CFD software, include the capability of producing high quality colourful pictures in order to help with the analysis of the solution of the flow field. Many of these post-processing techniques include, vector plots, contour plots and animations.
2.2
Governing Equations of Fluid Flow
The equations that govern the motion of a fluid are derived from simple principles such as the conservation of mass, conservation of momentum and conservation of energy. They can be formulated in differential or integral form, here the former formulation is used as many CFD.
The continuity equation, 1, which is derived from the principle of conservation of mass, expresses that for single-phase-fluids, mass cannot be created in a system, nor it can disappear.
∂ρ ∂t +
∂ ∂xi
(ρvi) = 0 (1)
The momentum equation, 2, is derived from Newton’s second law which states the fact that variation of momentum is caused by the net force acting on a mass element. The forces that are applied to a fluid volume can be divided into external or body forces (gravitation, buoyancy, etc) and surface forces that result from the pressure distribution and shear and normal stresses from friction.
∂ ∂t(ρvi) + ∂ ∂xi (ρvivj) = − ∂p ∂xi + ∂ ∂xi τij + ρfi (2)
The last governing equation, 3 is derived from the first law of thermodynamics. It states that any change in time of the total energy inside the volume is caused by the rate of work of forces acting on the volume and by the net heat flux into it. The total energy is the sum of internal energy, e and kinetic energy per unit mass v2/2. The energy equation is then:
∂(ρE) ∂t + ∂ ∂xi (ρviE) = − ∂ ∂xi (pui) + ∂ ∂xi kT ∂T ∂xi + ∂ ∂xi (ujτij) + ρuifi (3)
where the viscous stress τij, assuming that the fluid is Newtonian, is given by [12]: τij = µ ∂ui ∂xj + ∂uj ∂xi −2 3 ∂uk ∂xk δij (4)
In order to solve the system of equations, two additional equations are needed that relate the state variables (T, p, ρ).The formulation depends on the type of fluid considered. Two of the most common used methods are the perfect gas and the real gas formulations [13].
In addition to this, sometimes it is useful to solve an additional scalar transport equation that represents some flow parameter that it is not solved with the Navier-Stokes equations presented above. An example of this, is the concentration of any species that needs to be tracked through the domain. The equation to be solved is given by: ∂ ∂t(ρφ) + ∂ ∂xi ρuiφ − Γφ ∂φ ∂t = Sφ (5)
where Γφ and Sφ are the diffusion coefficient and source term for the scalar to be solved.
2.3
Turbulence Modelling
The strong chaotic motion of the flow under certain conditions causes the various layers of the fluid to mix together. Although these chaotic fluctuations of the flow variables are of deterministic nature, the simulation of turbulent flows still continues to present a significant problem [13]. Therefore, the effects of turbulence need to be accounted for in an approximate manner.
A classification of the different approaches (Fig. 6) can be done depending on the amount of turbulence that is resolved or modelled, as well as by the steady or unsteady nature of the models.
The classical turbulence model approach is based on the Reynolds-averaged Navier-Stokes (RANS) equations which models all turbulence. Thus, attention is focused on the effects that turbulence has on mean flow properties [11]. Because of this, these models include simplifications that lead to approximate solutions in complex flows keeping the required computing resources low for reasonable accuracy.
RAN S
U RAN S
LES/U RAN S/RAN S (DES, SAS) LES DN S SRS All turbulence modelled All turbulence resolved Steady U nsteady
Fig. 6: Different methods of turbulence modelling. Adapted from [14]
On the other side of the spectrum lies Direct Numerical Simulation (DNS) which aims at resolving all turbulence by solving the unsteady Navier-Stokes equations in a sufficiently fine spatial grid and small time step. As a consequence, the compu-tational cost of DNS is very high and the compucompu-tational resources required would exceed the capacity of the most powerful computers currently available for typical engineering problems [15].
The gap between RANS and DNS approaches is bridged by Scale Resolving Simulation (SRS), where only some of the turbulent scales are resolved. In this cat-egory Large Eddy Simulation (LES) is included, where the Navier-Stokes equations are spatially filtered and only the larger turbulence scales are resolved, while the effect from the smallest scales are modelled. The demand on computing resources for LES is large [11] and as a consequence new hybrid LES/RANS models have been developed, such as Detached Eddy Simulation (DES). These new models are capa-ble of switching from a LES behaviour in the flow regions where the grid spacing is sufficiently fine to RANS mode in the computational domain regions where LES is not required or the grid is not suitable.
2.3.1
Reynolds-Averaged Navier-Stokes Equations and
Classical Turbulence Modelling
The first approach for the approximate treatment of the turbulent flows was pre-sented by Reynolds in 1895 [13]. The methodology is based on the decomposition of the flow variables into a mean ( ¯Φ) and fluctuating part (Φ0).
Φ = ¯Φ + Φ0 (6)
There are three different forms of the Reynolds averaging: time averaging, ap-propriate for stationary turbulence; Spatial averaging, apap-propriate for homogeneous turbulence and ensemble averaging, suitable for general turbulence.
In case the density is not constant, it is recommended to apply the density (mass) weighted or Favre decomposition to some quantities. The most convenient way is to use Reynolds averaging for density and pressure and Favre averaging for other variables [13].
The Favre average of the velocity, for instance, can be obtained from the relation: ˜ vi = 1 ¯ ρT →∞lim 1 T Z t+T t ρvidt (7)
where T is a chosen period of time, and ¯ρ is the Reynolds averaged density. The Fravre decomposition is then:
vi = ˜vi+ vi00 (8)
where ˜vi represents the mean value and vi00 the fluctuating part of the velocity. Applying the Reynolds averaging to density and pressure, and Favre averaging to the remaining flow variables in the governing equations, the averaged continuity and momentum equations are:
∂ ¯ρ ∂t + ∂ ∂xi ( ¯ρ ˜vi) = 0 (9) ∂ ∂t( ¯ρ ˜vi) + ∂ ∂xj ( ¯ρ ˜vjv˜i) = − ∂ ¯p ∂xi + ∂ ∂xj ˜ τij − ¯ρ gvi00vj00 (10) These set of equations are identical to the Navier-Stokes equations, 1 and 2, with the exception of an additional term τFij = ¯ρ gvi00vj00, which is called the Favre-averaged Reynolds-stress tensor.
Here is where the Boussinesq hypothesis is introduced, 11. This hypothesis as-sumes that the turbulent shear stress depends linearly on the mean rate of strain, as in laminar flow, with a proportionality factor which is the eddy viscosity, a modelling parameter. The task of the turbulence model is to compute the eddy viscosity µT.
τFij = 2µTS˜ij− 2µT 3 ∂˜vk ∂xk δij− 2 3ρ˜¯kδij (11) where Sij = 1 2 ∂ui ∂xj + ∂uj ∂xi (12) The k-ε Model
The k − ε model belongs to the two-equation model class, where two additional transport equations are solved for two turbulence variables, in this case the turbulent kinetic energy (k), given by 14 and turbulent dissipation (ε), defined by 13.
ε = 1 ¯ ρτij ∂u00i ∂xj (13) k = 1 2u 0 iu0i (14)
The model is the most widely used and validated turbulence model. The model performs well in confined flows where the Reynolds shear stresses are most impor-tant. This includes a wide range of flows with industrial engineering applications, which explains its popularity [11].
However, the model shows only moderate agreement in unconfined flows, it is reported not to perform well in weak shear layers and it will have problems in swirling flows, secondary flows in long non-circular ducts and adverse pressure gradients [11]. The additional transport equations for turbulent kinetic energy, 15, and turbu-lent dissipation, 16, are in their differential form [13]:
∂(ρk) ∂t + ∂ ∂xj (ρvjk) = ∂ ∂xj µL+ µT σk ∂k ∂xj + τFij · Sij − ρε (15) ∂(ρε) ∂t + ∂ ∂xj (ρvjε) = ∂ ∂xj µL+ µT σk ∂ε ∂xj + C1ετFij · Sij− C2ερ ε2 k (16)
The turbulent eddy viscosity is then computed as: µT = ρCµ
k2
ε (17)
This model has five constants which have been obtained by data fitting for a wide range of turbulent flows:
In order to improve this model and bypass some of its disadvantages, some changes have been done. One of the variations of the Standard k − ε model is the Realizable version which contains a new formulation for the turbulent viscosity and dissipation source term. The term Realizable means that the model satisfies cer-tain mathematical constraints on the Reynolds stresses, which are consistent with the physics of turbulent flows. It introduces a variable Cµ to take into account for example, the mean strain and and the turbulence parameters (k, ε). A benefit of the Realizable model is that in almost every comparison, it demonstrates a superior ability to capture the mean flow of the complex structures [16].
The k-ω SST Model
The k − ω SST model also belongs to the two-equation model class, where the two additional turbulence variables are the turbulent kinetic energy (k) and specific turbulent dissipation (ω) defined as:
ω = 1
ε (18)
The k-ω SST turbulence model merges the k −ω model with the k −ε formulation in order to combine the positive features of both models. The k − ω approach is employed in the inner parts of the boundary layer while the k − ε is used in the outer region of the boundary layer as well as the free stream, far from the wall.
One distinct feature of the SST turbulence model is the modified turbulent eddy-viscosity function. The purpose is to improve the accuracy of prediction of flows with strong adverse pressure gradients and of pressure-induced boundary layer separation. The transport equations for the turbulent kinetic energy and the specific dissi-pation of turbulence are in differential form [13]:
∂(ρk) ∂t + ∂ ∂xj (ρvjk) = ∂ ∂xj (µL+ µTσk) ∂k ∂xj + τFij· Sij− β∗ρωk (19) ∂(ρω) ∂t + ∂ ∂xj (ρvjω) = ∂ ∂xj (µL+ µTσω) ∂ω ∂xj +Cωρ µT τFij · Sij − βρω2+ 2(1 − f 1) ρσω2 ω ∂k ∂xj ∂ω ∂xj (20) The turbulent eddy viscosity is computed by:
µT = ρk ω 1 max(a1 1, Sf2 a1ω) (21) The coefficients of the turbulence model β, Cω, σk and σω are obtained by blending the coefficients of the kωmodel, denoted as φ1with those of the transformed k − ε model, φ2.
Finally, the model constants are:
σk1= 0.85 σω1= 0.5 β1= 0.075 Cω1= 0.533 σk2= 1.0 σω2= 0.856 β2= 0.0828 Cω2= 0.440 a1 = 0.31 β∗ = 0.09 κ = 0.41
2.3.2
Scale Resolving Simulation
Even though most of today’s CFD simulations are based on RANS turbulence mod-els, it is becoming clear that certain flows are better predicted by models in which all or part of the turbulence spectrum is resolved in some regions of the domain as Scale Resolving Simulation [15]. As a consequence, SRS models provide additional information and increased accuracy.
However SRS models are very challenging in their proper application to indus-trial flows. The models typically require special attention to various details such as: model selection, grid generation and numerical settings, among others [15]
Scale Adaptive Simulation
The Scale-Adaptive Simulation (SAS) is an improved RANS formulation capable of resolving turbulent structures if the flow is unstable. For instance, in highly separated regions or vortex interactions. However, in attached or mildly separated boundary layersand undisturbed pipe flows, the SAS model returns a steady solution. The difference between standard RANS and SAS models lies in the treatment of the scale-defining equation. In classic RANS models, the scale equation is modelled based on an analogy with the k-equation using dimensional analysis. On the other hand, the SAS model is based on an exact transport equation for the turbulence length scale, with the variable φ =√kLt for the scale equation [15]:
∂(ρk) ∂t + ∂ ∂xj (ρvjk) = Pk− Cµ3/4ρ k φ2 ∂ ∂xj µT σk ∂k ∂xj (23) ∂(ρφ) ∂t + ∂ ∂xj (ρvjφ) = φ kµTSijSij ζ1− ζ2 Lt Lvk 2 − ζ3ρk + ∂ ∂xj µT σφ ∂φ ∂xj (24) The turbulent eddy viscosity is computed by:
µT = Cµ3/4ρφ (25)
The auxiliary relations to complete the model are
Lt= √ k Cµ1/4ω , Lvk = κ U0 U00 , U0 = S = q 2 ¯SijS¯ij, U00= s ∂2Ui ∂x2 k ∂2Ui ∂x2 j (26)
The model constants are:
Large Eddy Simulation
Large Eddy Simulation (LES) is based on a spatial filtering operation, which de-composes any flow variable Φ into a filtered part ¯Φ and into an unresolved part Φ0, such that:
Φ = ¯Φ + Φ0 (27)
The reason of performing this spatial filtering is based on the concept of resolving the larger scales of turbulence, which are difficult to model and problem dependent. On the contrary, smaller scales are more universal and isotropic, thus they can be modelled more easily [15].
The filtering operation is defined as: ¯ Φ = Z ∞ −∞ Φ(~x0)G(~x − ~x0)d~x0, Z ∞ −∞ G(~x − ~x0)d~x0= 1 (28) Applying this filter to the incompressible form of the Navier-Stokes equations, the continuity equation becomes:
∂ ¯vi ∂xi
= 0 (29)
whereas the filtered momentum equation is given by: ∂ ∂t(ρ ¯vi) + ∂ ∂xi (ρ ¯viv¯j) = − ∂p ∂xi + ∂ ∂xi ¯ τij+ τLESij (30) The equations feature an additional stress term due to the filtering operation called the subgrid stresses:
τLESij = ρ ¯viv¯j− ρvivj (31) As mentioned above the basic idea in LES is to resolve (large) grid scales (GS) and to model the (small) subgrid scales (SGS). Therefore, a subgrid model is needed for the turbulent scales which cannot be resolved by the grid and the discretisation scheme. The simplest model is the Smagorinsky model [17]:
τij − 1
3δijτkk= −2(Cs∆) 2S ¯S
ij (32)
where ∆ is a measure of the grid spacing, and S = q
2 ¯SijS¯ij is the strain rate scalar and Cs is a constant.
Detached Eddy Simulation
The Detached Eddy Simulation (DES) approach is to automatically operate either in RANS or in LES mode depending on the local grid spacing and resolved flow characteristic. Typically, RANS models are used in the boundary layer while the LES treatment is applied to the separated regions.
This capability is introduced by a function that can switch between methods, based on the ratio between the turbulent integral length scale and the local grid spacing [15]: Lt CDES∆max ( < 1, RAN S mode ≥ 1, LES mode (33)
where CDES is a constant and ∆max is the maximum edge length of the cell. Therefore when the grid resolution is larger than the turbulent integral length scale the model operates in RANS mode. On the other hand, when the grid spacing is finer than the integral length scale, LES mode is activated.
The integral length scale Lt is calculated from the RANS turbulence model magnitudes, for example from the k − ε model the length scale becomes:
Lt= k3/2
ε (34)
The dissipation term is then modified in the turbulent kinetic energy transport equation to account for the switch capability.
Yk= ρk3/2 Lt max 1, Lt CDES∆max (35) It is important to note that the DES limiter can be activated by grid refinement inside the attached boundary layers. This can lead to Grid-Induced Separation, where the boundary layers can separate depending on the grid spacing which is an unphysical effect. In order to avoid this limitation, the DES concept has been extended to Delayed Detached Eddy Simulation (DDES) by protecting the boundary layer from switching to LES mode with a shield function FDDES [15].
Thus the turbulent length scale for the k − ε model is redefined as: Yk= ρk3/2 Lt max 1, Lt CDES∆max (1 − FDDES) (36)
2.3.3
Wall Treatment
Turbulent flows are significantly affected by the presence of walls where the velocity field tends to zero and large gradients of flow variables appear. Therefore it is important to have an accurate representation of the flow in that region for a correct prediction of near-wall flows. However, this is not always feasible as a fine and dense mesh is needed close to the wall. Sometimes, this region is not of interest as it has no significant effects on the flow field.
The boundary layer is typically divided into two main regions: the inner (near-wall region) and the outer regions. The former can be further divided into three layers (Fig. 7): the innermost layer, viscous layer, is characterised by an almost laminar flow and the molecular viscosity plays a dominant role in momentum and heat or mass transfer; the outer layer, fully-turbulent layer, is where turbulence plays a major role and the buffer layer a bridge region between them where the molecular viscosity and turbulence are equally important [18]. The outer layer is dominated by the inertia effects from the main flow and insensitive to viscous effects.
There are two approaches in modelling the near-wall region, that depend on the mesh resolution near the wall. When the mesh is not sufficiently fine, characterised by a large non dimensional wall distance (y+), semi-empirical formulas called wall functions are used to bridge the viscosity-affected region between the wall and the fully turbulent region (Fig. 7). Examples of this are standard, scalable and non-equilibrium wall functions. The main disadvantage of all wall functions is that the numerical results deteriorate under refinement of the grid in the wall normal direction [18].
For instance the velocity profile can be estimated with the law of the wall, Fig. 7, in the fully turbulent layer by:
u+ = 1 κln y
+ + B (37)
and close to the wall in the viscous layer:
u+= y+ (38)
where y+ = u∗y ν , u
∗ is the friction velocity defined as qτW
ρ and u+ = u ∗/u. Additionally κ is the von Karman’s constant and B is another constant.
y+ 10-1 100 101 102 103 u + 0 10 20 30 40 50
Viscous Layer Buffer Layer
Fully Turbulent Layer
u+ = y +
u+ = 2.5ln(y +) + 5.45
Fig. 7: Subdivisions of the Near-Wall region, y+= u∗y/ν, where u∗ is the friction velocity
defined asqτW
ρ , u
+= u∗/u. Adapted from [18]
In another approach, near-wall modelling, the turbulence models are modified to enable the viscosity-affected region to be resolved all the way to the wall.
In addition to this, there is another different approach available in ANSYS Fluent called Enhanced Wall Treatment. The model divides the domain into two sub-domains, a viscosity-affected region and a fully turbulent region, where the equations are solved differently. If the mesh is fine enough the equations all the way down to the wall are solved. On the other hand, if the mesh is not sufficiently fine, enhanced wall functions will be used, which are an improvement of the standard wall functions.
3
Method
In this chapter a description of the method followed in this thesis is described. The method followed here is based in the approach by Alstom Power under ICAS-GT2 project [4]. Some enhancements are developed in order to improve its performance.
3.1
Geometry and Mesh
The geometry of the fluid domain together with the mesh are created with commer-cial software ANSYS R ICEM CFD v18.0.
The geometry generated is only a 2-D section of the cavity in order for it to be easier to mesh. After the 2-D section of the mesh is created a rotation is done in order to obtain the full 3-D meshed model (Fig. 8).
Fig. 8: 3-D mesh of the full cavity model
The 2-D sector of the domain is meshed using a multiblock strategy, in order to obtain a structured and high quality mesh (Fig. 9). The objective is to create a wall function mesh, as it has been seen that the use of wall functions is sufficient in order to capture the pressure instabilities [3]. Thus, the first node is placed at y+> 30 so that wall functions do not provide an inaccurate solution.
Another important aspect is to have a sufficient grid resolution in the region near the rim seal where the process of ingress and egress occur. The gap is meshed with 20 elements in the radial and axial direction.
Once the 2-D mesh is obtained, it is rotated along its axis in order to obtain the 3-D model of the cavity. However, this is not a straightforward process. Due to a problem in ANSYS R ICEM CFD v18.0, a two step method is used. First only a 45◦ rotation is used, with 25 elements in the tangential direction. After this rotation is done, the mesh is copied 7 times to obtain the full model (Fig. 8).
Fig. 9: 2-D section of the meshed cavity and some regions in detail
Furthermore a mesh sensitivity study was carried out in order to asses how the results depended on the mesh resolution. Two successive refinements were done from the initial mesh, increasing the number of elements by 20% and 66% respectively in the cavity region. The case selected for mesh independence study is Cw = 7000 and Realizable k − ε turbulence model.
All three meshes were created taking into account their quality: a low aspect ratio for the cells, a minimum angle of 27◦ and the quality measure in ICEM CFD (determinant 3 × 3 × 3) being higher than 0.4. The limit values for the different meshes are gathered in Tab. 4. It is worth mentioning that the quality in the region of interest, the cavity, is close to excellent an the values presented in Tab. 4 are the worst values found in the domain.
Table 4: Characteristics of the meshes tested showing the number of cells, the minimum element quality and angle as well as maximum aspect ratio
Mesh properties
Number of cells Min. Quality Max. AR Min. Angle
Coarse 1,510,992 0.514 36.59 31.79◦
Medium 1,810,944 0.503 33.69 31.79◦
Fine 2,499,336 0.482 33.06 30.18◦
The results of the mesh independence study compared are the amplitude and frequency of the three most dominant pressure modes in the cavity. Results show that there are no significant differences between the three meshes (Fig. 10). Bigger differences are seen in the amplitude with errors between 2-3% than for the frequency (≈ 1%), Tab. 5.
Taking into consideration these results, the coarser mesh is selected. Even though, a coarser mesh could be tested, the resolution near the rim seal and the cavity is adequate with this mesh. Taking into consideration how the mesh is
con-structed, lowering the resolution would lead to coarsening the region in the rim seal, which is undesired.
(a)
Coarse Medium Fine
Amplitude [Pa] 0.7 0.9 1.1 1.3 High Freq Mid Freq Low Freq (b)
Coarse Medium Fine
f/f R [-] 1 2 3 4
Fig. 10: Amplitude (a) and frequency (b) of the three most dominant modes for three different meshes
Table 5: Error, as a percentage, of the frequency and amplitude with respect to the results of the previous mesh, of the three most dominant modes in the cavity
Frquency (f /fR) Amplitude (kPa)
Medium Fine Medium Fine
Mode 1 2.99 3.14 1.18 0.78
Mode 2 3.61 3.22 1.02 0.76
Mode 3 2.12 1.42 1.16 0.95
3.2
Boundary Conditions
The boundary conditions considered for the cavity studied are rather simplified. The reason for this is that it has been showed [4], that some degree of simplification produces a reasonable fluid flow structure that resembles what experiments capture. Despite of this simplification, imposing boundary conditions that are closer to a real stator-rotor interaction, leads to results closer to experiments as Jakoby et al. showed [4]. However, this information is not available at the time, hence simplified conditions are used.
Uniform and steady conditions are considered in all boundaries of the cavity (Fig. 11). The main gas path inlet is characterised by a constant axial and tangential velocity together with the static temperature. At the outlet constant pressure with radial equilibrium is imposed and the total temperature of the possible backflow.
The purge flow inlet is considered to have no swirl thus only the axial component of the velocity is prescribed together with the static temperature.
For this three boundaries, turbulence properties are also needed. The chosen pa-rameters to specify are the hydraulic diameter together with the turbulence intensity,
calculated with 39, obtained from [19]. Here the simplification of fully developed flow is made. I = u 0 ¯ u ≈ 0.16(ReDH) −1 8 (39)
The Reynolds number based on the hydraulic diameter is calculated by: ReDH =
ρuDH
µ (40)
Fig. 11: Boundaries for the cavity geometry considered in this problem. Only showing a sector.
The cavity walls, the stator and rotor, are considered as adiabatic and no-slip walls. The stator is stationary while the rotor has a fixed rotational speed, Ω = −942.5 rad/s. Furthermore, the shroud is consider as symmetry.
Finally, as a scalar equation is solved for the concentration, boundary conditions are also needed. For the stator and rotor, a zero flux condition is set, while for the main annulus a constant value of 1 is prescribed at the inlet and 0.98 at the outlet, as seen from experiments. The purge flow has a value of 0 for the scalar.
All the boundary values used for the annulus inlet, outlet and purge flow inlet are gathered in Tab. 6. The velocity of the purge flow air, as well as the turbulence intensity are not constant and depend on the purge flow level (CW). Five different cooling mass flows are studied and the different inlet values in the purge flow channel are gathered in Tab. 7. The dimensionless mass flow (CW) is calculated by:
CW = ˙ mcool
µb (41)
where ˙mcool is the purge mass flow, µ is the air viscosity and b is the radius at the rim seal.
Table 6: Boundary conditions properties for the annulus and purge flow inlet
Boundary conditions
Annulus Inlet Annulus Outlet Purge Flow Inlet
cx [ms−1] 64 - ucool[2] cθ [ms−1] -228 - 0 cr [ms−1] 0 - 0 T [K] 348 376[1] 298 P [P a] - 231728 -Concentration 1 0.98 0
Turbulence boundary conditions
I [%] 2.56 2.56 Icool[2]
DH [m] 0.06 0.06 0.01
[1] In this case it represents the total temperature if any backflow appears [2] The velocity and intensity depend on the purge mass flow specified (Tab. 7)
Table 7: Velocity and turbulence intensity boundary conditions for the inlet purge flow depending on the mass flow.
Inlet purge flow boundary values
CW 7000 10000 12000 15000 25000
u [ms−1] 3.61 5.16 6.19 7.74 12.89
I [%] 5.48 5.24 5.12 4.98 4.67
3.3
Solver Settings
All the calculations are done using the stationary reference frame. Even though previous studies have used a rotating formulation due to a better convergence ([4], [5]), it was found, in this case, that convergence was not only more difficult to achieve but the results were not close to the experiments and previous simulations. The pressure-based solver is chosen, as this is a low Mach number problem, with an absolute velocity formulation.
The energy equation is solved together with continuity, momentum and turbu-lence equations. Furthermore, a scalar variable is added and its transport equation is also solved. This scalar represents the quantity of hot gas present in the cavity and it helps to track the amount of ingestion.
The material used in this study is air with an ideal gas formulation and con-stant specific heat, thermal conductivity and viscosity. Regarding the operating conditions, the outlet pressure is set as the operating pressure.
Furthermore, two types of turbulence modelling are used, URANS and SRS. The settings and the discretisation schemes are different for each model and discussed below.
URANS turbulence models
The Realizable k − ε model is chosen due to its improved performance with respect to the Standard k − ε model [19]. Moreover, it has been seen that this model gives results within an acceptable range with the experimental results ([4], [5]). The model constants are set to their default values in FLUENT. Standard and enhanced wall functions were used during the initial simulations, yet the latter are chosen for the successive computations. The reason is that it is recommended to use an insensitive wall treatment in case there is a further refinement of the mesh.
In addition to this, the model k − ω SST is also compared. Despite that the model full potential is achieved by covering the boundary layer with 10-20 cells and the mesh used here does not fulfil this requirement, evidence of previous successful work has been found Wang et al. by [7]. As with the previous case, the model constants are set to their default values in FLUENT.
The QUICK scheme is used for all variables, as it might provide better accuracy for rotating or swirling flows than second order and it also gives better convergence rate [4]. Least squares cell based scheme is used for the gradient which is the default and recommended in FLUENT. It is known to be as accurate as the node-based gradient for irregular unstructured meshes, but less expensive to compute [19]
Body weighted scheme is used for the pressure interpolation as it is recommended for problems involving large body forces [19]. The PISO scheme is used for the pressure-velocity coupling. Even though it takes more CPU power, it can decrease the number of iterations needed for convergence, specially in transient problems [18]. Unsteady RANS simulation are started from a steady state converged case, while steady state simulations are started with standard initialisation with values from the inlet main flow.
Time advancement is done with second order implicit integration. The advantage of using an implicit scheme is that it is unconditionally stable with respect to time step size [18]. The time step was carefully study taking into account that the CFL number should be ≈ 1. Thus the time step chosen is 6.66 µs, where one rotor revolution corresponds to 1000 time steps.
Under relaxation factors are set as default as they give a good convergence rate. If they are increased closer to 1, the solution and residuals oscillate, making conver-gence more difficult.
Scale Resolving Simulation
SAS, based on the k − ω model and DDES, on the Realizable k − ε, are the next models chosen under the Scale Resolving Simulation category. They can provide more accuracy in the results without increasing the computational costs as much as a LES simulation would.
The standard numerical schemes developed for RANS models cannot be used in SRS due to their dissipative nature. A numerical scheme with low numerical dissipation needs to be used to avoid the damping of the smallest scales [15]. Thus, Bounded Central Difference (BCD) scheme is chosen. First order upwind scheme is used for the spatial discretisation of turbulence, as the discretisation is not critical in SRS [15]. The rest of the variables were discretised with the QUICK an pressure interpolation is done with the body weighted scheme.
The gradient method chosen is the Least Square Method. Even though the gradi-ent method is not of much relevance on SRS simulations on high quality hexahedral meshes this method is recommended for the SAS model [15] and pressure-velocity coupling is done with SIMPLEC scheme, which is the optimal choice in most cases [15].
Time integration is carried out with the Bounded Second Order Implicit Euler scheme which has proven of sufficient accuracy for a wide range of applications [15]. The time step is also chosen taking into account that the CFL number should be lower than 1. The time step is therefore similar to URANS simulations.
All cases are started from a converged solution obtained from a URANS simu-lation after 8 rotor revolutions.
Under relaxation parameters also depend on the mesh chosen. However, for SRS they can be increased to values closer to 1.
3.4
Convergence Criterion
The convergence criterion for the simulations varies according to the type of analysis performed, whether it is a steady state or transient simulation.
Two conditions are needed to be fulfilled in order to achieve convergence for the steady simulations. The first one is that the RMS of the scaled residuals need to reach a value lower than a set target, 10−4 in this case. The second and more im-portant condition is that the flow inside the domain reaches a steady state situation where the properties do not change significantly in each iteration. In order to asses this condition, pressure and temperature are monitored in the domain.
Iterations (a) 0 1000 2000 3000 Pressure [kPa] -10 0 10 20 30 40 Iterations (b) 0 500 1000 1500 2000 2500 3000 RMS 10-4 10-2 100 101 Continuity x-velocity energy k
Fig. 12: Monitored static pressure (a), at a point in the cavity, and (b) residuals values during the steady state simulation
Despite this, these two conditions could not be fulfilled simultaneously. The residual levels for some equations, such as continuity, do not achieve the convergence criterion and show some oscillations (Fig. 12b). However, the properties monitored, for instance the pressure, show a stable solution after a certain number of iterations (Fig. 12a). This difficulty in convergence is due to the unsteady nature of the problem.
The transient simulation convergence approach is slightly different. For converge inside each time step, a decrease in two orders of magnitude is achieved for the scaled residuals in less than 20 iterations per time step, more iterations do not decrease significantly the residuals and add computational cost. The most difficult equations to converge are continuity and the transport equation for turbulent kinetic energy (Fig. 13b). The overall solution convergence is assessed by monitoring the pressure at some points in the cavity, for instance at the stator wall (Fig. 13). Once a periodic pattern is seen, the solution is said to be converged. A converged solution is obtained after around 3 rotor revolutions for URANS simulations and slightly longer for SRS cases.
Time [Rev] (a) 0 2 4 Pressure [kPa] -15 -10 -5 0 5 Iterations (b) 0 50 100 150 200 RMS 10-6 10-4 10-2 100 101 Continuity x-velocity energy k
Fig. 13: Monitored static pressure in time (rotor revolutions) (a), at a point in the cavity, and (b) residuals values during the transient simulation
3.5
Evaluated Cases
One operating point is only considered in this work as well as one geometry config-uration for the rim seal.
Five different purge flow levels are simulated for this configuration with turbu-lence model Realizablek − ε: CW = 7000, CW = 10000, CW = 12000, CW = 15000 and CW = 25000.
In addition to this, a total of four turbulence models are tested: Realizable k −ε, k − ω SST , DDES and SAS. The simulations are carried out only for one purge flow level, CW = 7000
3.6
Post-Processing
The post-processing of the results is carried out with the objective of comparing it to previous experimental analysis [3].
In order to do the comparison, pressure and concentration is tracked at several points in the cavity, as in the experiments. The concentration in this case is already a
non-dimensional parameter (given by a scalar equation) and the pressure coefficient is given by:
Cp =
2(P − Pref)
ρΩ2b2 (42)
where Pref is the reference pressure, Ω is the rotational speed and b is the rim radius.
Additionally radial (βR) and swirl (βT) velocity ratio, given by 43 and 44 re-spectively, are also used to present the results.
βR= cr Ωr (43) βT = cθ Ωr (44)
where Ω is the rotational speed of the rotor and r is the local radius.
The unsteady pressure analysis is done with the Fast Fourier Transform function (FFT ) given in MATLAB R 2015a. The location of the points where the unsteady pressure is tracked during the simulation (Fig. 14) is the same as in the experimental analysis [3].
PGap PHub
4
Results
4.1
Realizable k-ε model
4.1.1
Unsteady Results
The evolution of static pressure, relative to the outlet pressure boundary condition (Tab. 6 and15), in a monitoring point located at the stator wall at r/b = 0.63 shows a stable and periodic pattern. The pattern established after three revolutions since the start of the unsteady calculations and it differs depending on the level of purge flow. Low purge flow levels (CW = 7000) are characterised by low frequency and high amplitude oscillations, while for higher purge mass flows, (CW = 25000) the amplitudes are higher in frequency but lower in amplitude.
Time [Rev] 0 1 2 3 4 P [kPa] -20 -15 -10 -5 0 5 10 15 C W = 7000 CW = 25000 Monitor Point
Fig. 15: Static pressure evolution, relative to the outlet pressure, for two different purge flow levels, CW = 7000 and CW = 25000. The monitoring point is located in the stator wall
at r/b = 0.644, the red point in the figure.
Instantaneous pressure, concentration and temperature contours are shown in order to easily asses the fluid flow characteristics in the cavity and main annulus (Fig.16, 17 and 18 respectively). These contours are on an axial plane mid-point between the rotor and stator and they are taken after 8 revolutions.
The properties in the cavity are heavily affected by the conditions at which the turbine is running. The pressure contours show that for high and low purge flow levels the properties at the main annulus are quite uniform ( 16). Moreover, at the rim gap, alternate regions of ingress and egress occur around the circumference.
On the other hand, the picture inside the cavity is different depending on the purge mass flow. For low purge flow levels, (CW = 7000), two rotating low pressure structures appear in the inside of the cavity (Fig. 16a). They rotate at 80% of the rotor speed)and have an impact on the ingress of hot gas towards the cavity. Hot gas from the main gas channel is driven into those rotating low pressure areas (Fig. 17a) and as a consequence the temperature is also higher (Fig. 18a).
Fig. 16: Static pressure contour, relative to the outlet pressure, at a mid-plane in the cavity and annulus for two cooling mass flows: (a) CW = 7000 and (b) CW = 25000.Taken after 8 revolutions
Fig. 17: Hot gas concentration contour at a mid-plane in the cavity and annulus for two cooling mass flows: (a) CW = 7000 and (b) CW = 25000.Taken after 8 revolutions
Fig. 18: Temperature contour at a mid-plane in the cavity and annulus for two cooling mass flows: (a) CW = 7000 and (b) CW = 25000.Taken after 8 revolutions
However, when the level of purge mass flow is increased beyond a certain limit, the pressure field inside the cavity is nearly uniform (Fig. 16b). This situation leads to no ingress of hot gas into the cavity (Fig. 17b) and therefore, the temperature in this case is also homogeneous in the cavity (Fig. 18b).
As mentioned previously, alternate regions of ingress and egress occur around the perimeter of the rim seal. A clear approach to present this information is by plotting the instantaneous radial velocity ratio (βR) and swirl ratio (βT), given by Eq. 43 and 44 respectively, at different radius (Fig. 19). The data presented is also taken after 8 revolutions.
Vane Phase [-] (a) 0 4 8 12 16 β R [-] -1 -0.5 0 0.5 1 r/b = .55 r/b = .68 r/b = .93 Vane Phase [-] (b) 0 4 8 12 16 β T [-] -0.5 0 0.5 1 1.5 2
Fig. 19: Instantaneous (a) radial and (b) tangential velocity ratio at different radial positions along the full circumference for CW = 7000. The vane phase is defined as θ/22.5◦
and βR and βT by 43 and 44, respectively. Taken after 8 revolutions
Radial velocity profiles for low purge flow levels (CW = 7000) show that near the rim seal, the radial velocity ratio is small and does not have large variations around the circumference (Fig. 19a). Despite this, some regions of ingress (βR< 0) can be seen around vane phases 4 and 10. In addition to this, egress (βR> 0) of gas from the cavity into the main annulus is seen around vane phases 8 and 14. Lower radius, show clear and distinct regions where the flow field changes its radial direction, from the main gas path into the channel (βR< 0) and vice versa (βR< 0). These lower radius show higher radial velocity than upper regions in the cavity.
Additionally, swirl ratio has a different behaviour depending on the radial dis-tance (Fig. 19b). The swirl ratio is considerably higher close to the rim seal than at lower parts of the cavity. Moreover, at lower parts of the cavity the flow decelerates and rotates in the opposite direction of the rotor (βT < 0), which means that this radial position is under the low pressure region. A detailed vector contour of the low pressure structure shows that the flow inside it recirculates and there are regions under it where the tangential velocity is opposite to the rotor speed (Fig. 20).
As a point of comparison, high purge flow levels (CW = 25000) do not show a similar behaviour, as expected. Because of the non-ingress situation into the cavity, no low pressure regions appear thus the radial velocity and swirl ratio profiles show stable patterns through the circumference (Fig. 21). Near the rim seal, the flow is characterised by high frequency oscillations but with low amplitude.
Fig. 20: Instantaneous pressure contour and detail of velocity vectors in the low pressure region for low purge flow level, CW = 7000. Vectors are coloured by the static pressure and they are normalised. Even though they do not represent the real velocity magnitude, their direction is kept. Taken after 8 revolutions
Vane Phase [-] (a) 0 4 8 12 16 β R [-] 0 0.02 0.04 0.06 r/b = .55 r/b = .68 r/b = .93 Vane Phase [-] (b) 0 4 8 12 16 β T [-] 0.06 0.07 0.08 0.09
Fig. 21: Instantaneous (a) radial and (b) tangential velocity ratio at different radial positions along the full circumference for CW = 25000. The vane phase is defined as θ/22.5◦
and βR and βT by 43 and 44, respectively.Taken after 8 revolutions. Note that different
scales have been used in comparison to Fig. 19
The concentration contours of hot gas in the rim seal vicinity (Fig. 22) show what has been previously mentioned, a different behaviour depending on the purge flow level. Low purge flow levels (CW = 7000) are characterised by a large ingestion of hot gas into the cavity(Fig. 22a), specially through the stator side.
High purge flow levels (CW = 25000) however, are characterised by egress of gas from the cooling flow into the main gas channel disturbing the flow, which is not desirable either (Fig. 22b).
The difference between stator and rotor hot gas concentration can be easily appreciated through the concentration contour for CW = 7000 (Fig. 23). The con-centration in the stator is higher due to the disk pumping effect which moves the flow upwards through the rotor wall and downwards through the stator causing a larger ingestion of hot gas into the latter.
