Turbulence Models for Simulation of the Flow in a Rotor-Stator Cavity

Modelling of the flow in the cavities between rotor and stator in turbomachines (e.g. pumps or turbines) is a task of great interest. Correctly evaluated pressure and velocity fields enable calculation of the disk losses and therefore assessment of efficiency. It is also crucial for determination of axial thrust and thus design of the bearings. The study demonstrates abilities of various turbulence models to describe the flow in a narrow gap between rotating and stationary disks. Numerical simulations were performed in order to find out the ability of particular models to capture unstable structures appearing during specific operating conditions as well as to calculate the velocity profiles precisely. Large Eddy Simulation (LES), Scale Adaptive Simulation (SAS), Detached Eddy Simulation (DES), Reynolds stress model (RSM) and SST k − ω model were used. Obtained results were also compared with experimental measurement published by Viazzo et al. [1].


Introduction
The flow of a fluid between rotating and stationary disk received attention due to its applicability to many industrial and scientific problems. The investigations started at the beginning of the 20th century by study of Ekman [2] describing wind-driven ocean currents. Ekman was followed by Von Kármán [3], who obtained a similarity solution of Navier-Stokes equations for an infinite disk rotating in quiescent fluid. Bödewadt [4] followed his predecessor with study of flow near a stationary disk in rotating fluid. These findings became the background for comprehensive study of flow and boundary layers established between stationary and rotating disks.
Smith [5] noticed disturbances in form of waves or vortices dependant on Reynolds number. He carried out experiments with rotating disk boundary layer and observed fluctuations in a narrow range of Reynolds numbers below the transition to turbulence. The instability is referred to in literature as Type 1, Type B or crossflow instability. The instability was experimentally visualised and it was found out, that it appears in form of spiral vortices rotating anticlockwise. Faller [6] described the second instability occurring in lower Reynolds numbers, which is known as Type 2 or Type A instability. It also forms as spiral vortices, however it has opposite direction to Type 1 instability.
Instabilities in the vicinity of the rotating disk were first studied experimentally, i.g. [7][8][9][10], later on, numerical studies [11][12][13] occurred. With development of computing technology, computational fluid dynamics (CFD) became a widespread tool for flow simulation. New challenges connected with modelling of turbulence arose. Different approaches exist, nevertheless, the choice is always * e-mail: lucie.zemanova1@vut.cz a compromise between the accuracy and computational resources.
From the accuracy point of view, the best results can be achieved by direct numerical simulation (DNS). The turbulence is explicitly resolved at all relevant length and timescales. From that fact results also a disadvantage, since the computational domain has to be large enough to include the largest naturally-occurring eddies while a computational grid has to be fine enough to fully resolve all dissipation scales. This leads to limited usability for relatively low Reynolds numbers, which are orders of magnitude smaller than in common industrial applications. DNS is not implemented into commercial software due to the requirement of special numerical schemes [14].
For computation of engineering problems, simpler description of the turbulence based on averaging in time is used. Instead of finding the instantaneous flow field, Reynolds-averaged Navier-Stokes (RANS) equations are solved. The set of the equation is a result of Reynolds decomposition, which divides the variables into timeaveraged and fluctuating component. The problem of concept of averaging is that it produces additional variables, therefore the system of equations is not possible to solve. This is in literature referred to as the closure problem. By supplementing additional equations (turbulence models), RANS can be solved. Many different types of turbulence models exist. Depending on the number of additional equations, we can distinguish (0-4)-equation models. The two-equation models based on Boussinesq eddy-viscosity hypothesis, such as k − or k − ω are the most common in practical applications [15]. Common k − model offers good trade-off between the accuracy and computational costs. It is based on the two additional transport equations for turbulent kinetic energy k and turbulent dissipation . In its basic form, it belongs to so-called high-Reynolds models group. As the name suggests, it is applicable to high turbulent Reynolds numbers. Turbulent Reynolds number is given by the ratio of eddy viscosity to molecular viscosity, therefore it is a measure of the level of turbulence. In the areas with low Re, where the viscous effects dominate (e.g. near wall regions), high-Reynolds models are not appropriate and other approach needs to be employed. The common practise it to use wall-functions, where the integration is not performed up to the wall, but in these areas, it is replaced by empirically obtained equations. Later, modifications of k − which enables to solve the flow field even in the near wall regions without the wall-functions were developed (low-Reynolds version of k − model). However, the wall functions as well as the integration up to the wall are based on the empirical relationships derived for boundary layers developed during the steady flows over a flat plane. Therefore, validity in case of the complex flow involving boundary layer separation, pressure gradient, suction, blowing, roughness etc. is questionable [16].
Better results for solution of the flow field in a boundary layer gives low-Reynolds two-equation k − ω model. It adds two transport partial differential equations for k and ω (rate of the energy dissipation) to deal with the closure problem. It is suitable for flows with adverse pressure gradients and when the integration through the viscous sublayer is preferred (e.g. for solving boundary-layer transition) [17]. Modification of the model derived by Menter [18] SST k − ω is often implemented in commercial software due to the effective combination of high-and low-Reynolds approach. It couples the k−ω and k− turbulence model in a way that the k − ω is used in the region of the boundary layers and switches to the k − in the core flow. The next step further to describe precisely the turbulence is Reynolds Stress Model (RSM). Unlike the eddy viscosity models, it takes into consideration anisotropy of the turbulence. The individual components of the Reynolds stress tensor are directly computed, therefore it is able to resolve more complex turbulent interactions such as strongly swirling flows. Wall-function as well as low-Reynolds approach can be applied in connection with RSM [15].
Different approach to simplification of the Navier-Stokes equations than Reynolds averaging in time offers large eddy simulation (LES). Eddies of large scale are solved directly using Navier-Stokes equations, while the small scale (sub-grid scale) eddies are modelled. The outlined procedure leads to significant reduction of the computational costs compared to DNS. It is still more demanding than RANS, nevertheless this approach gives results of much better accuracy. The large eddies, which carry the most of the turbulent energy and are responsible for the most of the momentum transfer are captured directly. On the other hand, sub-grid scale eddies are from their nature more isotropic and homogeneous, therefore the modelling is sufficient enough to obtain high accuracy solution [19]. Moreover turbulence models for isotropic turbulence can be relatively simple. The disadvantage is a need to use time resolution orders smaller in comparison with RANS. Similar situation is also for the spatial resolution of the grid. Wall functions are not possible to use and strict re-quirements on element quality in boundary layer regions has to be met [20].
In order to overcome the computational demands of LES for higher Reynolds numbers, so-called hybrid models combining LES with RANS were developed. Detached eddy simulation (DES) utilizes LES for resolving the detached eddies (separated regions) far from boundaries and RANS model for the near wall flow. It exhibits good results in flows with large separation regions, cavities with simple geometry or flows connected with acoustic problems. On the other hand, the weakness is description of curvature streamlines and strong dependency of the results on the mesh [21]. The regions of separation solved by LES should be known prior to the simulation and meshed in an appropriate way, as describes [22]. Slightly different approach offers scale-adaptive simulation (SAS) method, which brings dynamic behaviour to the model. By introducing Von Kármán length scale into the scaledetermining equation of RANS turbulence models, automatic balancing of the contributions of modelled and resolved parts of the turbulent stresses is enabled. As a consequence, for unstable flows the model changes smoothly from LES model through various stages of eddy-resolution to steady RANS model [23].
For CFD simulations of the flow inside rotating cavities, different authors use various approaches to modelling the turbulence. Chew [24] solved the velocity profile of the flow inside the rotating cavity with high-Reynolds k − model, while authors in [25] used the low-Reynolds approach, Elena and Schiestel [26] involved RSM model. All these studies were focused on the solution of velocity profiles and no instabilities were mentioned. By comparison with measurements it was found out, that these models underestimate the boundary layer thickness. Poncet et al. [27] reported relatively good agreement with experiments by RSM, while in [28] there is described successful deployment of SST k − ω to rotating cavity flow problem. Lo et al. [29] compared measurements with simulations using LES. Apart from other turbulence models, they were able to detect even the instabilities emerging on rotor and stator disks. From its nature, promising results with reduced computational resources should be achieved using hybrid LES-RANS approaches.
Taking into consideration the need to capture instabilities, the structure of the flow is highly complex involving laminar, transitional and turbulent regions. Moreover, the turbulence is strongly inhomogeneous and anisotropic due to rotation effects. It leads to a very challenging task for turbulence models [30]. The aim of this study is to explore and compare the capabilities of particular models to capture instabilities in rotor-stator cavities and their capability to describe the velocity profile precisely. Resolving the flow in boundary layers is crucial for determination of friction torque and efficiency, since it directly influences the dissipation on the shroud and the hub, so-called disk friction. The findings will be later on applied on the flow in real shape rotor-stator cavity of a centrifugal pump and simulations leading to the design optimization.

Numerical model
The present study is based on the same geometry of rotating cavity as was analysed by Viazzo et al. [1]. Different approaches to modelling the turbulence were applied on this case, namely LES, SAS, DES, RSM and k − ω. Commercial software ANSYS Fluent 19.1 was used to perform the calculations.

Geometry and mesh
The analysed geometry with corresponding dimensions is shown in Fig. 1. It consists of two parallel disks of outer radius b = 140 mm. The disks are delimited by an inner cylinder (the hub) attached to the rotating disk and by an outer cylinder (the shroud) connected to the stator. The height of the gap is h = 20 mm. Due to the symmetry of boundary conditions and the resulting flow, it is possible to reduce the domain to the fraction and save computational resources, as was shown in [1]. The flow can be characterized by the aspect ratio of the cavity G = 5, the curvature parameter R m = 1.8 and the rotational Reynolds number Re = 4 × 10 5 , defined as: For the purposes of evaluation, dimensionless axial and radial location is defined: Where z * = 1 gives the position of the stator and z * = 0 is a location of the rotor. Dimensional radial location r * = 0 describes the hub and r * = 1 is the shroud. The mesh which would meet requirements even for LES was created as shown in Fig. 2. The spatial resolution was 180 × 120 × 80 elements in radial, tangential and axial direction with refinements near walls, enabling low-Reynolds approach. The mesh fulfills demands described in [31], which means x + 100 (stream-wise direction), y + 1 (wall-normal direction) and z + 30 (span-wise direction).

Boundary conditions
No-slip boundary condition was applied to all of the walls of the cavity. The shroud and the upper disk was stationary. The lower disk and the hub were set as rotating walls with constant angular velocity Ω = 20.4082 rad.s −1 , corresponding to the rotational Reynolds number Re = 4 ×10 5 , for which the instabilities occur. Periodic boundary condition was applied to the side walls of the domain.

Computational details
The fluid inside the rotating cavity was considered incompressible. The default model for water from ANSYS library was used. The calculation was first run as a steady state with k − turbulence model. SIMPLE scheme was used for pressure-velocity coupling. First orders of accuracy were set for advection terms in all transport equations. Enhanced wall-treatment (integration up to the wall without wall functions using one-equation model in the nearwall region) was used. After achieving the convergence, pressure, turbulent kinetic energy and turbulent dissipation rate was switched to the second orders of accuracy and momentum to QUICK scheme (3 rd order of accuracy). The results were then used as an initial condition for transient analysis. Pressure-velocity coupling algorithm was changed to PISO, which is recommended for unsteady flows. An appropriate time step was chosen with respect to the known nature of instabilities and particular model of turbulence.

Large eddy simulation (LES)
Based on the results previously reported in the literature, e.g. [1], [29], the time step for LES enabling to capture coherent vortical structures was set to 1×10 −5 s. Smagorinsky-Lilly subgrid-scale model was chosen. Second orders of accuracy with bounded central differencing method for momentum was considered.

Scale Adaptive Simulation (SAS)
The time step was initially set to 1×10 −5 s as well as in LES case. It was expected, that coarser time resolution would be sufficient, however the setting was used in order to ensure better convergence and to obtain an initial guess for Courant number.
where U is the magnitude of velocity, ∆t is time step and ∆x is the minimum cell size. According to Courant-Friedrichs-Lewy (CFL) condition, Courant number should not exceed 1 for explicit numerical method. Based on that, maximum time step for obtaining correct results is 4×10 −4 s. Second orders of accuracy were applied as in the LES case. Scale adaptive simulation was used in combination with the SST turbulence model.

Detached Eddy Simulation (DES)
According to CFL condition, even slightly larger time step can be applied in DES, however the difference is not significant, therefore the same time step as in previous case ∆t = 4×10 −4 s was applied. As a RANS model for resolving the parts of the domain where no separation is detected, SST k − ω was chosen. The case was solved with the second order of accuracy. DES used SST k − ω RANS model.

Reynolds stress Model (RSM)
For RSM, linear pressure-strain model was set, the second order of accuracy schemes for advection terms were applied. Time step, which ensure C < 1 is again ∆t = 4×10 −4 s.

SST k-ω
As in previous cases, appropriate ∆t = 4×10 −4 s and the second order of accuracy was set.

Instabilities
Coherent structures emerging in the boundary layers can be detected by Q-criterion (second invariant of velocity gradient tensor). For given rotational Reynolds number, Type II instability is expected, which can be observed as spirals in 15 • angle in tangential direction, rotating in opposite direction than the disk. Fig. 3. Experimental visualisation of the instabilities [33] In rotor-stator cavity, Ekman boundary layer is formed on the rotor, while Bödewadt type of the boundary layer appears on the stator. Bödewadt boundary layer is less stable than Ekman and therefore, instabilities occur first on the stator side and are more pronounced. The ability of particular turbulence models to capture the instabilities is demonstrated by figures in which Q-criterion on the stator is shown.
As can be seen in Figs. 4 and 5, capability of LES to describe the instabilities was confirmed. The spiral vortices in the middle part as well as stronger unstructured turbulence near the hub described in experiment [1] is captured. The other turbulence models failed to detect the coherent structures in rotor-stator cavity, as document Fig. 6-9. In case of SAS, indications of vortices can be seen in near hub region, where the strongest intensity of the turbulence takes place. DES, RSM and SST k − ω were not able to model spiral structures entirely. Even when the time step was lowered to the 1×10 −5 s as in LES case, no satisfactory results were obtained.

Velocity profiles
The flow inside the rotor-stator cavity can be characterized as a core with two thin boundary layers. The thickness of Ekman boundary layer created on the rotor is significantly lower than Bödewadt boundary layer on stator. It corresponds to the results published in [1] and [32]. The fluid is driven centrifugally outwards along the rotor, turns into axial direction and due to conservation of mass, it is forced to flow radially inwards along the stator.
The velocity of the core fluid can be divided into axial, radial and tangential directions. The axial velocity component approaches zero, therefore it is not visualised. Fig.  10 shows instantaneous radial (left) and tangential (right) velocity component.
The velocity of the core fluid can be divided into axial, radial and tangential directions. The axial velocity component approaches zero, therefore it is not visualised. Fig.  10 shows instantaneous radial (left) and tangential (right) velocity component.
The profiles were evaluated in the middle of the computational domain (mean radius of the disk cavity). Both, the axial coordinate and the velocity are dimensionless. The components of velocity were made non-dimensional according to following relationships: Reasonable agreement between numerical predictions and experimental data was reached. All mentioned models were able to calculate velocity in the core region precisely. The closer the walls, the more the discrepancies are visible. Unfortunately, in these regions, only few measured data were obtained and the measurement was averaged over time, thus unsteady phenomena were suppressed. As mentioned earlier, LES was the only model, which was able to capture instabilities. It propagates into the velocity profile, as can be observed especially on the rotor side (z * = 0). This is the reason for the difference in velocity profile between experimental data and LES model towards the rotor. In comparison with measured data, numerical models slightly underestimate the thickness of Ekman boundary layer and overestimate the thickness of Bödewadt layer. Velocity profiles in near wall region can be seen in detail in following figures. Fig. 11 shows radial and tangential velocity close to stator (Bödewadt boundary layer). In Fig. 12, similar results for rotor (Ekman boundary layer) are presented. Correct resolution of velocity profile in inner boundary layer is crucial in determination of the shear stress and consequently friction losses.
In Figs. 11 and 12, only numerical models are compared, since the closest measured point has y + out of the viscous sublayer -the main area of interest. It can be seen, that the velocity profiles computed by SAS, DES, RSM and SST k − ω models do not differ too much between each other in general. They are able to describe the profile with good accuracy in case of tangential velocity (see Figs. 11 and 12 right).
Even though the instabilities were not captured, the slope of the tangential velocity in viscous sublayer (y + <5) is very similar to the LES model, which included the unstable nature of the boundary layer.
In buffer layer region (5 < y + < 60), the differences are most significant. Out of the buffer layer, with increasing y+ the velocity profiles obtained by different models approach each other. Larger discrepancies are visible in case of radial velocity (Figs. 11 and 12 left). The shape of the velocity profile has the opposite direction on the stator and the rotor side, since the fluid flows radially inwards along the stator, whilst outwards orientation can be observed on the rotor. The slope of the curve differs in viscous sublayer as well as in buffer layer. SST k − ω, DES and SAS produce similar result on the stator as well as on rotor. RSM slightly differs on the more unstable stator side. The velocity profiles computed by these models are shifted from LES profile. On the stator, they predict larger gradient of radial velocity in viscous sublayer, on the other hand, on the rotor side the gradient of radial velocity is smaller in comparison with LES. The radial velocity magnitude from SAS, DES, RSM and SST k − ω calculation is smaller in viscous sublayer and buffer layer of the stator compared to LES. On the contrary, on the rotor side, the predicted radial velocity by SAS, DES, RSM and SST k − ω is larger than LES modelled. In this paper, capabilities of LES, SAS, DES, RSM and SST k − ω turbulence models to describe the flow inside rotor-stator cavity were demonstrated. The ability to capture the vortices emerging in unstable regime was assessed. For further practical application meaning calculation of the disk friction losses and efficiency, the numerical model should be able to solve velocity profiles precisely, especially in near wall regions. Therefore, the second part of the study dealt with comparison of velocity profiles obtained by particular turbulence models. It was found out, that less computationally demanding models, i.e. SAS, DES, RSM and SST k − ω, are not able to detect the instabilities at all. The vortical structure obtained by LES agrees with experimentally reported data in [1]. The model is able to capture the spiral vortices on rotor and stator as well as stronger turbulence with unstructured swirling near hub.
The radial and tangential velocity profiles computed by SAS, DES, RSM and SST k − ω are very similar. They agree with LES in the core of the fluid, however, in the near wall regions larger differences are obvious. The largest discrepancies were observed in radial component of velocity. Unfortunately, these regions are of the great interest for determination of friction losses. Also, these regions are important for the onset and further spreading of the instability.
In conclusion, the comparison between the models showed, that LES in inevitable level of modelling for such flows. Reduction of computational demands by using different turbulence model leads to significant reduction of accuracy and to the loss of the important information about instabilities.
Turbulence Models for Simulation of the Flow in a Rotor-Stator Cavity was supported by project FSI-S-17-4615.