Assessment of advanced RANS turbulence models for the stability analysis of low specific speed pump-turbines

Reversible pump-turbines, which can switch between pump and turbine operation, provide important support for the stabilization of the electrical grid. In order to fulfill their task, fast transition between the operation modes is required. This includes stable synchronization with the electrical grid during the process of turbine start at small guide vane openings (GVOs). However, the special geometry of pump-turbine impellers tends to produce unstable flow behavior, which can affect the synchronization process in a negative way. The prediction of the stability behavior with the help of CFD during the design process is not yet fully developed. The simulation of small specific speed pump-turbines operating at small GVOs sets high requirements on the underlying numerical model. To capture all relevant flow structures accordingly, 3-dimensional unsteady computation of the entire pump-turbine domain is required. Further, it has been shown that turbulence modeling is a key factor in order to reproduce the 4-quadrant characteristic in the unstable S-shape region. Vortex structures in the vaneless space and runner channels are considered the root-cause for the occurring instabilities. Standard two-equation models are not able to predict the stability behavior under these difficult circumstances (highly separated flows, small specific speed, small GVO angle). Models, which consider the anisotropic behavior of Reynolds-stresses have shown to produce better results. Considering the accurate prediction of mixing processes at off-design conditions, the focus was oriented to the anisotropy present in the turbulent structures. Standard models are often not sufficient to accurately predict vortices, which can have a huge impact on the performance, since based on the assumption of isotropic turbulence. Accordingly, they tend to dissipate and diffuse vortices too quickly. Improved models, which take the anisotropic nature of the Reynolds-stresses into account, can help in this context. The models can introduce the anisotropy via an explicit algebraic expression, or model directly the transport equations for the Reynolds-stresses. For this work the performance of a full Reynolds-stress model with improved coupling was investigated. In order to obtain the predictions using advanced turbulence models a particularly robust framework based on a pressure-based fully coupled approach was used. The goal of this work is the development and testing of an improved full Reynolds-stress model for the application in the stability analysis of pump-turbines. The focus lies thereby on separation behavior and mixing/vortex dissipation.

Abstract. Reversible pump-turbines, which can switch between pump and turbine operation, provide important support for the stabilization of the electrical grid. In order to fulfill their task, fast transition between the operation modes is required. This includes stable synchronization with the electrical grid during the process of turbine start at small guide vane openings (GVOs). However, the special geometry of pump-turbine impellers tends to produce unstable flow behavior, which can affect the synchronization process in a negative way.
The prediction of the stability behavior with the help of CFD during the design process is not yet fully developed. The simulation of small specific speed pump-turbines operating at small GVOs sets high requirements on the underlying numerical model. To capture all relevant flow structures accordingly, 3-dimensional unsteady computation of the entire pump-turbine domain is required. Further, it has been shown that turbulence modeling is a key factor in order to reproduce the 4-quadrant characteristic in the unstable S-shape region. Vortex structures in the vaneless space and runner channels are considered the root-cause for the occurring instabilities. Standard two-equation models are not able to predict the stability behavior under these difficult circumstances (highly separated flows, small specific speed, small GVO angle). Models, which consider the anisotropic behavior of Reynolds-stresses have shown to produce better results.
Considering the accurate prediction of mixing processes at off-design conditions, the focus was oriented to the anisotropy present in the turbulent structures. Standard models are often not sufficient to accurately predict vortices, which can have a huge impact on the performance, since based on the assumption of isotropic turbulence. Accordingly, they tend to dissipate and diffuse vortices too quickly. Improved models, which take the anisotropic nature of the Reynoldsstresses into account, can help in this context. The models can introduce the anisotropy via an explicit algebraic expression, or model directly the transport equations for the Reynoldsstresses. For this work the performance of a full Reynolds-stress model with improved coupling was investigated.
In order to obtain the predictions using advanced turbulence models a particularly robust framework based on a pressure-based fully coupled approach was used. The goal of this work is the development and testing of an improved full Reynolds-stress model for the application in the stability analysis of pump-turbines. The focus lies thereby on separation behavior and mixing / vortex dissipation.

Introduction
The changes in the electricity market in recent years and the outlook to future developments force operators of hydraulic power plants to adapt their business model. The increased percentage of subsidized renewable energy, mainly solar and wind, in combination with a politically driven reduction of conventional power plants, leads to (1) lower energy prices in the marked and (2) the necessity to stabilize the electrical grid with alternative methods. Hydraulic storage power plants are predestined to overtake the role of balancing the electrical grid with tasks like storage, peak production and other ancillary services. All require an increased dynamic operation of the machines and correspondingly short reaction times, which can be met by pump-turbines. The drawback are additional mechanical stresses since a considerable part of the time these are operated at off-design conditions. It is therefore necessary to further understand the behavior of the pump-turbine under these circumstances.
Pump-turbine runners can present unstable behavior in turbine mode in particular close to the speed-no-load region (turbine and turbine brake) and the reverse pump mode. This can produce for example oscillations in the hydraulic system [1] and slow down the synchronization process. It is therefore important to assess already during the design phase the stability, i.e. before model tests.
To capture the relevant flow features with CFD it is necessary to simulate the full machine in 3-D and fully unsteady. This can be achieved by computing individual operating points (OP's) along the 4Q-curve characteristic, as already been done in [2,3]. The results largely depend on the turbulence model, the GVO angle and the specific speed of the pump-turbine [4,5].
The commonly used turbulence models are SST, SAS-SST, standard and RNG k-, v2-f. In the last few years also anisotropic approaches have been used, such as the EARSM turbulence model. In general, it could be observed that almost all turbulence models deliver good results for large GVO angles, i.e. they can capture the S-shape characteristics sufficiently well independent on the specific speed of the investigated machine [6,7,8]. On the other hand there are serious difficulties and differences between the models when computing the S-region at small GVO angles, for example in final phase of load rejection or close to synchronization conditions. Here especially the SST model seems to be inappropriate to accurately compute the S-shaped [2,9].
The aim of this work is to analyse the unstable behavior of a pump turbine with low specific speed and at small GVO using an improved Reynolds-stress model, which adopts coupling of the individual stress components. This model is integrated in the HSLU in-house pressure-based and fully coupled solver. The use of a numerically robust framework is crucial for the simulation of pump-turbines in the S-shape region, especially when using advanced turbulence models. The results are compared with measurements and afterwards analyzed in more detail to find new insights into the underlying phenomena of the S-shape characteristic.

Solver
The operating points simulated in this work are located in the region of the turbine brake, an area which deviates strongly from design conditions. These operating points show a strong incidence at the guide vanes and runner blades and are dominated by complex and highly transient flow processes in the vaneless space. Simulation under these conditions impose high requirements on the solver in terms of stability and robustness. Particularly since a Reynolds-stress model is used for the simulations, which can capture the occurring vortex structures more sharply than standard two-equation models and thereby causes less turbulent diffusion/dissipation, which could otherwise have a stabilizing effect on the simulation. Simulation speed is the second key criterion for the choice of the solver besides stability. Since a relatively large simulation model (see next section) has to be calculated transiently, the computing effort can be substantial.
The solver used here is the block coupled solver by Mangani et al. [11]. The governing 30th IAHR Symposium on Hydraulic Machinery and Systems IOP Conf. Series: Earth and Environmental Science 774 (2021) 012020 IOP Publishing doi:10.1088/1755-1315/774/1/012020 3 equations are coupled implicitly and therefore the robustness and code speed criteria can be met. The system of equations shown in Equation (1) is solved simultaneously. Cell-values are denoted by the subscript "C" and neighbour contributions by the subscript "NB".

Computational Domain, Mesh and Simulation Setup
The simulation is based on a model-sized pump turbine, for which the test rig is set up in the HSLU hydro laboratory. Therefore the simulation can be directly compared with corresponding measured values. The simulation model consists of several components and is divided into four sections (see figure 1). The first section (A) contains the spiral case (SC) and the stay vanes (SV). This is followed by (B) containing the guide vanes (GV) and the (C) the runner (RU). The simulation model does not consider the runner side spaces of the real installation. This means that friction losses and leakages are neglected. The draft tube (DT) in (D) forms the last section. In order to improve the convergence behavior of the simulation, the outlet was extended with a pinch. All mesh components are structured grids with a resolution typical for industrial application (i.e. medium mesh resolution for efficient simulation times). The overall cell count is 5 million, and a boundary layer resolution with an average y + of approximately 30 is applied. Figure 1 shows the mesh topology between runner and volute for a GVO angle of 6 • , which corresponds to synchronization condition.
The discretization of the flow equations is second order in time and space, with a backward Euler temporal scheme and a high resolution spatial scheme.

Turbulence modeling
The turbulence model used in this work based on [20] and is an improved coupled Reynoldsstress model. The benefits of the model are its favorable stability and convergence behaviors. A brief outline of the model is given next. The underlying equations of the model are (2)(3)(4)(5)(6). As base serves the Stress-Omega model of Wilcox [21] with modifications of the ω equation suggested by Menter [22]. It is a so called SSG Reynolds-stress model (Speziale, Sarkar & Gatski).
The exact production term is defined as and the pressure strain term as The pressure strain expression is usually divided into two parts, the slow pressure-strain (term 1.) and the rapid pressure-strain (terms 2-3). P k , the production term of the turbulent kin. energy and D ij , the turbulence transport term, are defined as follows For details on the model coefficients and blending functions the reader is referred to [22] and [21]. The implementation uses several coupling terms derived from the dissipation and slow pressure-strain expressions to improve the stability properties of the model [20].

Simulation control
For the boundary conditions the following settings were used: At the spiral case inlet the mass flow rate, turbulence intensity of I T = 0.05 and a viscosity ratio (ν t /ν) of 10 are specified. An area-averaged static pressure of p = 0 is prescribed at the draft tube outlet. All walls are defined as no-slip walls, either rotating or stationary. Considering the runner domain, a constant rotational speed of n = 1500 is imposed for all computations. The individual hydraulic components are numerically coupled by Arbitrary Mesh Interfaces (AMIs). Maximum simulation duration corresponds to ca. 20 runner revolutions. All evaluated quantities are averaged over a period of ca. 10 runner rotations. The time resolution during the unsteady simulations was set to 1/20 of a runner-channel rotation, which corresponds to ca 2.6 • per time step. The time-step was set to 0.29 ms. The maximum number of iterations per time-step was set to 10 with the convergence criteria of residual at each time-step was in the order of 1e-5. During the numerical simulations, the results of steady RANS simulations were used as the initial flow field.
3. Discussion of the 4Q-characteristic Figure 3 and 2 show the comparison between simulation and measurement for the 4Qcharacteristic curve in the S-shape region at a GVO angle of 6 • and 20 • . The definition of the coefficients K cm1 and K u1 are given in (7).  The comparison between the simulation results and the measurements reveals several findings. On the one hand, the hypothesis can be confirmed that anisotropic turbulence models are better in predicting the physical processes in the S-shape than standard two equation models [4,5] for small GVO angles. The slope of the simulated characteristic curve corresponds relatively well to the behavior observed for the measurements.
Also the starting point of the S-shape characteristic is better predicted than with other models such as k − ε. However, the simulation still shows a significant downwards shift with regards to the measurements. The precise detection of the beginning of the S-shape region is essential for the decision whether a geometry is stable or unstable with regard to synchronization. In the case of figure 3 the speed-no-load line of the measurement is included. If the S-shape occurs sufficiently far below this line, than the pump-turbine can be considered as stable. In the case of the simulation results however, it must always be taken into account that the speed-no-load line changes as well. An important factor for this shift has to do with ignoring the friction losses of the runner side spaces. In turbine operation the runner torque is therefore constantly overestimated by the simulation. This gives some explanation why the simulated speed-no-load line is lower in the 4Q diagram compared to measurements. For the simulation results shown in Figure 3, the line passes just below OP B. Thus in the presented case CFD correctly predicts that the runner is unstable, especially with regards to synchronisation.

Flow features 4.1. Vortex structures in the vaneless space and separation behavior at the GV
The figures 4-7 show 3-D streamlines starting at the GV leading edge for the operating points A-D. Before the S-shape is formed in OP A, a strong shift of the streamlines to the side of the shroud can be seen ( figure 4). This shift is caused by vortex structures in the vaneless space which force this behavior on the flow in the GV channel. OP B differs only slightly from OP A. Also here, shifting of streamlines towards the shroud occurs, however less intense ( figure 5). For OP C, the first point within the S-shape region, no shift is observable anymore ( figure 6). The streamlines are symmetrical in span-wise direction. OP D shows the flow structures deep inside the unstable S-shape region. Remarkable is the distinct formation of horseshoe-like vortex structures at the leading edge of the GV (figure 7). These are accompanied by a stabilization of the vortex structures in the vaneless space and runner channels, which are shown below.     For OP D, the horseshoe-like vortices occur simultaneously in all GV channels, regardless of the runner position. Figure 8 visualizes this fact by showing an iso-volume of the helicity in the vaneless space and GV channels. The GV blades are colored with the static pressure and the iso-volumes with the helicity. The uniform behavior in circumferential direction also explains why the pressure fluctuations decrease along the S-shape curve.  All operating points analyzed in this work are strongly in the off-design range and show therefore strong incidence at the RU and GV leading edges. The figures (14)(15) show the vortex structures on the mid-plane for the operation points A to D using the LIC (Line Integral Convolution) method. A comparison of the different figures reveals that from OP A towards OP D the vortices are getting stronger, which leads to an increased blockage of the runner channels. Furthermore, it can be recognized that the deeper the operating point is located in the S-shape region, the more uniform appear the flow pattern in circumferential direction. While for OP A  figure 14) individual channels still appear partially free, OP D ( figure 15) shows more or less identical blockage for all channels. The vortex structures have settled and are relatively stable. The increased throttling effect caused by the vortices leads to the characteristic S shape curve. A decrease in mass flow is thereby associated with an increased pressure drop.    Mid-plane view with LIC streamlines colored by the vorticity (axial comp.): OP D.
In order to get a better impression of the flow conditions, 3-D streamlines through the large clockwise-turning runner vortex are shown in figure 14 for OP A and in figure 15 for OP D. It can be recognized that blocking of the runner channels does not only occur at the height of the mid-plane, but over the entire span of the runner. The comparison between figure 14 and 15 shows again the growth of the clockwise-turning vortex along the S-shape. Further on the pictures another feature is visible. Parts of the toroidal vortex, known under the term water-ring in the literature [23], can be seen in the vaneless space. It seems that a considerable part of the water-ring is used to feed the clockwise-turning runner vortex.

Rotating-stall phenomenon at low discharge
During the simulation of operating points in the S-shape region it was possible to detect a special flow behavior, which is also visible in the experiment, the phenomenon of rotating stall. As already shown in [24], the investigated pump-turbine model shows at small GVO angles and low discharge (therefore in the S-shape region) a rotating stall behavior, which appears at the guide vane channels and travels along the circumference of the entire guide vane apparatus.   Figure 16 shows the normalized pressure signal recorded at the guide vane blade (GV) in the volute casing (SC) and in the draft tube (DT) during the simulation of OP C. The placed pressure probes in the simulation correspond to the set-up of the real test-rig (see also [24]). The definition of p norm is defined in eq. (8). In the GV and in a weaker form in SC, a low-frequency oscillation can be identified, which is superimposed on the pressure signal. This oscillation cannot be detected for DT. An FFT analysis confirms this ( figure 17). The peaks at 7 and 14 times the runner frequency correspond to the blade passing and a multiple of it. The low-frequency oscillation is at about 0.7 times the runner frequency, which shows good agreement with the measurements [24] and typical values found in the literature for rotating stall. This means that the simulation based on the Reynoldsstress turbulence model not only provides (1) good results for the 4Q-characteristic and (2) an insight in relevant flow features, but can also reflect the local effect of rotating stall according to the measurements, which reinforces the confidence in the validity of the simulation results. Other anisotropic models like EARSM are also able to capture the characteristic and local flow features, but could not detect the rotating stall phenomenon properly.
In order to visualize the rotating stall cell figure 9 shows the contour plot of the radial velocity at the GV for OP D. The inward flow is marked blue and the backflow region red. The results are only shown on a plane near the hub, but the behavior on this plane is almost identical to that near the shroud. Figure 9 clearly shows a non-uniformity of the flow along the circumference. In the left half of the picture the GV channels are mostly blocked. On the right side, almost all GV channels are free of blockage and show flow in the correct (inwards) direction. The flow pattern is time dependent. As previously noted, the obstruction travels along the circumference at 0.7 times the impeller frequency.

Conclusive Remarks
In this work the performance and predictions of a full RS model for the simulation in the S-shape region of a small specific speed pump-turbine at a small GVO angle was investigated. The applied full RS model uses an improved coupling between the different components of the Reynoldsstresses. It could be shown that the stability and robustness of the model are beneficial for the challenging calculation of transient processes in the S-shap region. The transient simulation of the different operating points could be carried out with the same set-up which was also used for the calculations with standard two equation turbulence models. The additional computational effort resulting from the increased number of turbulent equations (5 additional eq. compared to standard two equation models) can thus be kept within limits.
Global simulation results and measured values show good agreement with regard to the 4Q-characteristic curve, as already observed with other anisotropic turbulence models such as EARSM. In the present case the unstable runner is recognized as such, since the calculated speed-no-load line is located in the S-shape region. The full RS model therefore seems to be suitable in the design process of pump-turbine runner geometries.
In addition to the 4Q-characteristic curve, the local phenomenon of the rotating stall at the GV could be clearly identified in the simulations for low flow rates. Its frequency corresponds to the measured data in the test stand (approx. 0.7x runner speed). The evaluation of the simulation results shows that several adjacent guide vane channels are simultaneously blocked by vortices during this process.
The results generated with the full RS model can now be used to extend the knowledge of the flow processes in the S-shape region. In this work the analysis of vortex structures have been extended to the guide vanes. The observations lead to the assumption, that in addition to the runner, the flow in the guide vane also has a clear influence on the S-shape characteristic. This assumption will be further investigated in subsequent studies. Of particular interest are dynamic processes of the plant such as load-rejection and turbine-start with the full RS model. The integration of the CFD simulation into a complete 1-D model of the test rig for the investigation of actively excited oscillation states, for which exp. measurements have already been carried out, is another interesting task for future research.