Prediction of S-shaped characteristics in reversible pump-turbines using different numerical approaches

The prediction of characteristics and flow phenomena in reversible pump-turbines becomes increasingly important, since operations under off-design conditions are required to respond to frequency fluctuations within the electrical grid as fast as possible. Fulfilling the requirements of a stable and reliable operation under continuously expanding operating ranges challenges the hydraulic design and requires ambitious developments. Beyond that, precise estimations of occurring flow phenomena combined with a detailed understanding of their causes and mechanisms are essential. This study aims at predicting the S-shaped characteristics of two reversible pump-turbines by using different numerical approaches. Therefore, measurements at a constant wicket-gate opening of Δγ = 10° were done. Based on these experimental data, unsteady flow simulations are performed under steady and transient operating conditions respectively: Starting from the best efficiency point in generating mode, through the runaway, along the S-curve, down to operation in reverse pump mode. The hydraulic machines are spatially discretized in model size with a near-wall refinement of y+mean ≤ 5 and y+mean ≥ 30. The application of two different solvers discloses deviations in underlying methods. The turbulence modeling is basically executed by the k-ω-SST and the standard k-ϵ model. Focusing on higher order numerics, the Explicit Algebraic Reynolds Stress Model (EARSM) is selected in the commercial code and extended with an approach for curvature correction (EARSM- CC). In the open-source software, the four-equation v2-f model assumes the role of higher order numerics. The temporal discretization errors are observed using three different time-step sizes. As a supplement, experimental data obtained from the HydroDyna pump-turbine are used as additional validation, providing integral quantities and local pressure distributions at an operating point set on the S-curve. To sum this work up, a methodology is developed to approximate S-shaped characteristics and local effects within the hydraulic machinery properly.


Introduction
The varying energy consumption and the usage of volatile energy sources require temporary operation of hydraulic machinery at part-load or even at runaway to rapidly respond to fluctuating loads and to balance frequency shifts in the electric grid. At those off-design conditions, the flow field is dominated by unsteady flow phenomena which stress individual components and adversely affect life cycles of hydraulic systems. The operation of reversible pump-turbines might additionally become unstable in generating mode if so-called S-shape instabilities occur. This means that the current operating point (OP) switches back and forth between turbine (brake) and reverse pump mode. Their occurrences manifest as an unstable branch in characteristics, indicated by a positive slope of dQ ED /dn ED > 0 and dT ED /dn ED > 0. Such an S-shaped characteristic is depicted in Fig. 1. The definitions of the dimensionless coefficients read: Unstable conditions cause serious operational problems, particularly during the startup procedure where significant hydrodynamic effects delay or even prevent a successful synchronization. Another issue represents the case of load-rejections, at which large amounts of energy have to be transferred and dissipated within the hydraulic machine. As a consequence of all these facts, an accurate numerical prediction of occurring instabilities becomes crucial and numerical accuracy plays a major role in modern pump turbine design. It is of course an advantage to determine onset and intensity of S-shape related phenomena properly before proceeding with the model tests. Concerning this matter, several works have already been published that aim at precisely computing and investigating arising flow phenomena at offdesign conditions in generating mode [4,5,6,7,10,11]. The HydroDyna project -for example -was introduced to analyze the unstable behavior along S-shaped characteristics (S-curves) and provides a variety of experimental data and numerical results [5,7,11]. Apart from the expected RSI frequency (including its harmonics), the performed investigations reveal a dominant low frequency in turbine brake mode that represents a rotating stall in the runner passages. Yan et al. [11] achieve an appropriate approximation of those existing physical effects using the k-ε model and a time step size that corresponds to 1 • runner rotation (N ∆t /ς HD W G =18). Supplementing to that, this study focuses on the prediction of S-shaped characteristics and related local phenomena in reversible pump-turbines depending on different numerical methods. Based on experimental data of two reduced-scale model machines at a constant wicket gate opening (∆γ = 10 • ), the numerical investigations are executed by applying unsteady simulations under steady or transient operating conditions from best efficiency point, through runaway to reverse pump mode. The computed integral quantities are compared with model test data as well. The chosen settings are additionally validated with experimental data obtained from the HydroDyna pump-turbine scale model, serving integral quantities and static pressure distributions at runaway (OP-69) and at an operating point located on the S-curve (OP-73) [5]. Concerning the numerical setup, all physical models represent a spatial discretization of the entire hydraulic machines, see Fig. 2. Nevertheless, the presence of discretization errors is evident as geometrical simplifications and spatial constraints cause deviations to real conditions. Errors in grid consistency are estimated by a modified method of the Richardson Extrapolation [8] contributed by Celik [2]. The physics in the boundary layer is mainly affected by quantity gradients and shear stresses. The wall-function universally approximates near-wall effects and, therefore, causes inaccuracies that have to be quantified. In order to capture related deviations, the near-wall treatment is executed at two different near-wall refinement stages, namely y + mean ≤ 5 and y + mean ≥ 30. Apart from applying common two-equation turbulence models like the standard k-ε (k-ε) and the k-ω-SST (SST) model, the EARSM as higher order turbulence model reveals correlations between anisotropy of turbulent quantities and numerical accuracy [1]. In order to take streamline curvature and system rotation into account, the EARSM is extended by an approach for curvature correction. This implementation based on the work of Spalart and Shur [9] modifies the production term and sensitizes the turbulence models for these mentioned effects. Since the EARSM model has not been implemented in OpenFOAM yet, the four-equation v 2 -f model [3] is applied in the open-source software package instead. Temporal discretization errors are observed using time-step sizes of 1.2, 4 and 17/18 (depending on the considered model machine j) time-steps per wicket-gate pitch (N ∆t /ς j W G ). A comparison of two solvers (Ansys CFX 15.0 and OpenFOAM 2.3.x) discloses the influence of underlying numerics in either code.

Reversible pump-turbine scale models and discretization
Numerical investigations are carried out based on the reversible pump-turbine scale models RPT 1 , RPT 2 and HydroDyna (HD). Their computational domains represent the spatial discretization of the entire hydraulic machines and consist of a spiral case (SC), stay vanes (SV), a wicket gate (WG), a runner (RU) and a draft tube (DT). SV and WG grids are Table 1. Grid quality of RPT 1 and RPT 2 . The y + mean values are computed at runaway. in between. The related grids are illustrated in Fig. 2. Grid quality and machine specifications are given in Tab. 1 and Tab. 2. A determinant value of 1 indicates a fully orthogonal grid. Table 2.
Main specifications of RPT 1 , RPT 2 and HydroDyna. The specific speed n q is stated in turbine mode.
For the HydroDyna case the same model discretization is used as described in Yan et al. [11], whose contribution includes all relevant information about its attributes. The extensions at inand outlet restrain impacts of boundary conditions on the actual flow field.

Boundary conditions and evaluation settings
The boundary conditions set at the inlet comprise a discharge, a turbulence intensity of I T = 0.05 and a turbulent length scale of L t = 0.01m. The explicitly computed operating points are highlighted in Fig. 4. The following sensitivity study deals with a comparison of flow quantities at best efficiency point (BEP), runaway (RA) and zero discharge (Q0) on RPT 1 and at OP-73 on HydroDyna (labeled in Fig. 4). An area-averaged static pressure of p mean = 0 is prescribed at the outflow. A so-called opening condition allows a re-entering flow at the outlet and thus reduces numerically oscillations. The walls are defined according to the no-slip formulation. Considering a rotating runner domain, a constant rotational speed n   interactions which are needed to connect the rotating runner with its surrounding stationary domains. The interpolation methods used to describe convective movement is of second-order. Related implementations are the High-Resolution scheme (HRS) in CFX and a massflowweighted approach in OpenFOAM. The turbulence transport is approximated by an upwind difference scheme (UDS) to avoid divergence during the computations. In order to reduce computational costs, the maximum number of inner iteration loops is limited to 8-10 steps by an abortion criterion for the averaged values of the residuals. Moderate under-relaxations factors are chosen to stabilize numerical oscillations during the iterative process. Concerning the unsteady simulations under steady operating conditions, the maximum simulation duration corresponds to 40 runner revolutions. The thereby obtained solutions are used as initialization for the unsteady simulations under transient operating conditions. All evaluated quantities are averaged over a period of the last ten runner rotations. The elaboration of outstanding settings is part of the following sensitivity study.

Sensitivity study 4.1. Grid independence
In order to quantify the discretization error due to grid inconsistency, a grid independence study is done on RPT 1 at three different grid refinements (h i = 1, 1.3, 1.7) using CFX and the EARSM-CC turbulence model. Steady-state solutions are computed for the dimensionless coefficients Q ED and n ED at RA, extrapolated according to the Richardson Extrapolation [8] and compared with its estimated value f 0 : The quantities f 1,1.3 represent the computed values of Q ED and n ED at the different refinement stages h i and r ≥ 1.3 is the grid refinement ratio prescribed by Celik [2]. The order of accuracy p is specified by reaching values of p > 2 which denotes a sufficiently accurate method for the given case. The results show a good accordance with the asymptotic estimation using the fine and medium size grid h 1 and h 1.3 (see Fig. 5). However, the coarse refinement stage h 1.7 is not capable of accurately predicting the unknown quantities and yield increased deviations from f 0 . Focusing on a compromise of accuracy and computational costs, the hexahedral medium size grid with about 7.3 · 10 6 nodes has revealed itself as the best variant for all further investigations.

Near-wall refinements
The type of near-wall treatment might significantly affect the flow distribution in hydraulic machines since effects in the boundary layer are resolved differently. In order to investigate the influence of different near-wall refinement stages on integral quantities of RPT 1 , the medium size grid is adapted to reach y + mean ≤ 5 and y + mean ≥ 30. The resulting deviations r,a of Q ED , n ED and T ED from the measurements at BEP, RA and Q0 are shown in Fig. 6. As assumed, resolving the boundary layer reaches a higher accuracy than approximating physics with the universal wall-function -particularly at off-design conditions where gradients of flow quantities mainly affect the onset of flow separations and recirculation zones. However, the coarse resolution predicts S-shaped characteristics very well but causes slight offsets to the experimental data. Since the fine near-wall resolution requires about 50% more nodes than the coarse one, a final decision about the preferred option, therefore, has to reflect computational costs in addition to numerical accuracy. Due to the large number of CFD-simulations performed in this work, those numerical efforts are of major importance, which in turn requires the usage of the coarse refinement stage for all further computations.

Turbulence models and solvers
This section describes the dependence of calculation results on applied turbulence models and solvers. Therefore, a comparison of integral quantities is done in CFX and OpenFOAM using the k-ε and the SST model. Higher order numerics are considered by applying the EARSM W G = 17 to ensure convergence. The deviations of Q ED , n ED and T ED from the measurements are given in Fig. 7 and Fig. 8. Fig. 9   Exp.
Exp. overestimates the amplitudes at the RSI frequency and its second harmonic. Surprisingly, the k-ε model reaches acceptable results but cannot compete with the higher order numerics. Concerning the EARSM, the curvature correction tends to improve numerical accuracy at unstable operating conditions. Taking computational costs and stability issues into account, OpenFOAM requires, compared to CFX, significantly more computation time to reach similar results. As the EARSM-CC is proven to be most suitable to reproduce real flow conditions in the present reversible pump-turbines, it is used for all further investigations.

Temporal discretization
The choice of a proper temporal discretization determines the possibility of numerical methods to detect physics within a particular frequency band and to achieve convergence at each time step. However, certain flow phenomena need several runner revolutions of computation to fully develop. Here too, a balance between computational costs and accuracy should be achieved. The variation of time step sizes is intended to provide information about their impact on numerical results, considering values of 1.2 time steps per WG-pitch (N ∆t /ς j W G = 1.2), N ∆t /ς j W G = 4 and N ∆t /ς j W G = 17 or 18 (depending on the investigated model machine j). Fig. 10 shows the corresponding deviations of Q ED , n ED and T ED from the measurements at RA and OP-73, whereas the static pressure distribution, recorded at MP66 in the vaneless space of HydroDyna, is depicted in Fig. 11. The observations reveal a sufficiently accurate solution using the medium time step size (N ∆t /ς j W G = 4). The coarse discretization generates increased deviations of integral quantities and is not capable of predicting the dominant low frequency at MP66. Its temporal resolution restricts the upper limit for the frequency to a value that lies below the second harmonic of the RSI frequency. Since the fine discretization does not improve outcome compared to the medium time step size, the latter is selected as best variant to discretize the time-dependent terms of the underlying partial equations.  Figure 11. FFT of static pressure observed at MP66 using the coarse time step size N ∆t /ς HD W G = 1.2 (left), the medium time step size N ∆t /ς HD W G = 4 (center) or the fine time step size N ∆t /ς HD W G = 18 (right) for HydroDyna.

Prediction of S-shaped characteristics
Based on the final numerical setup depicted in Fig. 13, specific operating points of RPT 1 and RPT 2 are numerically predicted under steady or even transient operating conditions, starting from BEP, through RA, along the S-curve and down to reverse pump mode. The associated results are illustrated in Fig. 12  quantities occur in reverse pump mode and amount to r ≤ 3.1% for RPT 1 and to r ≤ 3.4% for RPT 2 . The corresponding discrepancies in turbine and turbine brake mode reach maximum values of r ≤ 2.3%.

Conclusion
Based on experimental data from two different reversible pump-turbine model machines, this work aims at developing a methodology to accurately predict S-shaped characteristics and occurring phenomena under the aspects of acceptable computational costs. A sensitivity study reveals the influence of specific input parameters on computed results. An additional validation is given by the HydroDyna's model test data [5], providing integral quantities and static pressure distributions in the vaneless space. Derived from these comparisons, a final numerical setup is developed (see Fig. 13) and S-shaped characteristics approximated under steady and transient operating conditions. The highest accuracy is yielded with Ansys CFX, combined with the spatial discretization of the entire hydraulic machine, a near-wall refinement stage of y + mean ≤ 5, an extended version of the EARSM considering curvature correction (EARSM-CC) and a temporal discretization of 4 time steps per wicket-gate pitch (N ∆t /ς W G = 4). As a compromise between accuracy and numerical stability, the high-resolution scheme (HRS) is used to interpolate convective terms, whereas the upwind difference scheme (UDS) describes turbulence transport through the cell boundaries. The selection of the type of near-wall refinement stage depends mainly on the self-defined limits for the computational costs since both techniques gain acceptable results. Applying the final numerical setup together with the coarse near-wall resolution, a good accordance between numerics and measurements is apparent for both reversible pump-turbines, reaching maximum relative deviations of r ≤ 3.4% for Q ED and n ED observed in reverse pump mode. Related discrepancies in turbine and turbine brake mode amount to r ≤ 2.3%. The continuous development in numerical methods manifests itself in a significant improvement in prediction quality compared to previous studies. The unsteady simulations under transient operating conditions reproduce the measured S-shaped characteristics with acceptable accuracy and supply beneficial information about occurring phenomena along the S-curve.

Acknowledgment
Measurements on the HydroDyna pump-turbine were performed within the frame of the HydroDyna project (Eureka N • 3246) which represents a cooperation of Alstom Hydro, Andritz, Voith Hydro and UPC-CDIF. The authors would like to thank all involved partners for providing model test data.