CFD Investigation of a High Head Francis Turbine at Speed No-Load Using Advanced URANS Models

Featured Application: The present numerical investigation gives some insights regarding the challenges to compute the ﬂow in hydraulic Francis turbines at off-design operating points. However, such simulations are useful to design hydraulic turbines that will be operated with more ﬂexibility. Therefore, the present results can be used to determine which turbulence models are useful to compute off-design operating points, keeping in mind limitations and future improvements. Abstract: Due to the integration of new renewable energies, the electrical grid undergoes instabilities. Hydroelectric power plants are key players for grid control thanks to pumped storage power plants. However, this objective requires extending the operating range of the machines and increasing the number of start-up, stand-by, and shut-down procedures, which reduces the lifespan of the machines. CFD based on standard URANS turbulence modeling is currently able to predict accurately the performances of the hydraulic turbines for operating points close to the Best Efﬁciency Point (BEP). However, far from the BEP, the standard URANS approach is less efﬁcient to capture the dynamics of 3D ﬂows. The current study focuses on a hydraulic turbine, which has been investigated at the BEP and at the Speed-No-Load (SNL) operating conditions. Several “advanced” URANS models such as the Scale-Adaptive Simulation (SAS) SST k − ω and the BSL- EARSM have been considered and compared with the SST k − ω model. The main conclusion of this study is that, at the SNL operating condition, the prediction of the topology and the dynamics of the ﬂow on the suction side of the runner blade channels close to the trailing edge are inﬂuenced by the turbulence model.


Introduction
Hydropower is the main renewable source of energy generation in the world with 3.7% of the total final energy consumption [1]. In order to increase the part of the hydropower in the energy mix, efforts are made to develop micro-hydropower [2], to recover energy from existing water networks [3], or to improve the efficiency of the existing power plants during refurbishment [4]. Beyond these efforts and due to the development and integration of renewable energy resources, hydraulic turbines and pump-turbines are key technical components for load shifting and frequency control of the grid [5,6]. However, in order to act as a grid stabilizer, the hydraulic machines have to extend their operating range [7].
Computational Fluid Dynamics (CFD) is one of the main tools used to investigate the dynamics of the flow in various hydraulic turbines at both stable and unstable operating points [8,9]. Most of the 2 of 26 numerical investigations are carried out by solving the Unsteady Reynolds-Averaged Navier-Stokes (URANS ) equations. The Reynolds stresses are often modeled using Boussinesq's assumption, which requires computing the eddy viscosity, usually calculated using a two-equation turbulence model. Such an approach has been used to compute full head Francis turbines at the Best Efficiency Point (BEP) including also full load [10][11][12][13], partial load [13][14][15], and Speed-No-Load (SNL) [16][17][18] operating points or transition from pump to turbine mode [19], start-up [20], or the S-shaped region of pump-turbines [21]. The Francis-99 workshop, whose summary is available in [22,23], gave an overview of the state of the art for the simulation of high head Francis turbines at design and off-design operating points. Regarding the turbulence modeling, the authors mentioned that the BEP is well predicted by both URANS and hybrid turbulence models. The part load operating point is more challenging, and results show a larger discrepancy compared to the measurements. This discrepancy is explained by the influence of the mesh near the wall boundaries. Hybrid turbulence models improve the results at part load operating conditions compared to the standard URANS models.
Regarding more precisely the SNL operating conditions, several authors put in evidence the inability of standard two-equation turbulence models to capture the pressure fluctuations [17,18]. Indeed, these models are built to compute steady flows in absence of streamline curvature, stagnation points and large vortices, which are characteristic of the flow at SNL. Therefore, in order to overcome the limitations of the standard turbulence models, Scale-Adaptive Simulation (SAS) [17] based on the Shear Stress Transport (SST) k − ω model or Large Eddy Simulation (LES) [18] has been carried out. In both cases, the content in pressure fluctuations is enhanced and shows an improved agreement with the experimental measurements.
LES in a complete hydraulic machine requires a very fine mesh with around 100 million nodes [18,24]. In order to avoid the use of such large meshes, which require large Central Processing Unit (CPU) resources, some authors intended to use hybrid RANS/LES models. Hybrid RANS/LES models are designed to adjust the level of the eddy viscosity in relation to the grid spacing, at least in the framework of Detached Eddy Simulation (DES) [25,26]. Krappel et al. [27] used the Improved Delayed-Detached Eddy Simulation (IDDES) model to simulate the flow in a Francis pump turbine with a mesh containing 20 million elements. Compared to the results provided by the SST k − ω turbulence model, the IDDES model improved the prediction of the standard deviation of the velocity components and of the turbulent kinetic energy in the draft tube cone. Minakov et al. computed the part load and full load operating conditions in a reduced domain of the Francis-99 turbine using a DES SST k − ω model and a mesh with only six million control volumes [28]. The results show an improvement of the accuracy regarding the r.m.s. velocity components compared to the results provided by a standard Reynolds stress model. The LES behavior can also be obtained without using explicitly the grid spacing, as is the case in the SAS framework proposed by Menter and Egorov [29,30]. The SAS formulation has been used by several authors to compute the flow through hydraulic turbines [17,31]. Neto et al. [31] showed that the SAS model on a grid with nine million control volumes improves the prediction of the velocity profiles in the draft tube diffuser for the part load operating condition compared to the standard SST model. Excluding the resolution of the turbulent eddies, Explicit Algebraic Reynolds Stress Models (EARSM) and Reynolds Stress Models (RSM) are able to capture the flows characterized by stagnation points and streamline curvature. Such models have been used for instance by Mössinger et al. [32] to compute the flow in a high head Francis turbine including part load, BEP, and full load operating conditions. They concluded that the EARSM or RSM model does not improve the results compared to the SST model. Moreover, at part load condition, the RSM model was not able to converge.
Based on the previous literature review, the linear eddy viscosity models are able to compute accurately the BEP of a high head Francis turbine, but have some difficulties in computing the flow at part load and SNL operating conditions. During the last 15 years, several improvements have been achieved in the development of URANS and hybrid RANS/LES models [33][34][35]. Consequently, the present paper focuses on the influence of some "advanced" URANS turbulence models regarding the dynamics of the flow at the SNL operating condition of a high head Francis turbine already investigating at part load by Müller et al. [36] and at various operating points by the current authors [37][38][39].
The paper deals first with the description of the case study, the governing equations, and the numerical setup. In the BEP results, a detailed analysis regarding the influence of the mesh, the boundary conditions, and the numerical setup is carried out. Finally, the SNL operating point is computed with four "advanced" URANS models: the SST k − ω model with production limiter and curvature correction, the SAS SST k − ω model, the Organised Eddy Simulation (OES) SST k − ω model, and the Baseline k − ω model (BSL) EARSM model. The results are compared with the ones provided by the standard SST k − ω model and with experimental data.

Case Study
The considered case study is the Grimsel 2 pumped-storage power plant operated by Kraftwerke Oberhasli AG (KWO) [40] and built between 1974 and 1980 in Switzerland. This power plant is equipped with four 100-MW ternary groups: a horizontal-axis motor/generator, a Francis turbine, and a pump. The number of guide vanes is z g = 24, and the number of runner blades is z b = 17. The speed, discharge, and power factors along with the specific speed at the BEP are respectively: At the BEP, the Reynolds number is equal to 21.6 × 10 6 based on the discharge velocity and the diameter at the runner outlet section (these two quantities are used to define the dimensionless variables in the paper. The dimensionless pressure is defined with the head at the BEP.).
As the electricity market requires more and more flexibility, the number of start-up, stand-by, and shut-down procedures per year has increased strongly during the last decade. Consequently, the turbine spends more and more time at operating points for which the hydraulic conditions are far from the ones at the BEP. After several operation cycles, some cracks have been observed on the runner at the junction between the trailing edge of the blades and the hub. The phenomenon responsible for the development of the cracks is not yet clearly identified. However, experimental measurements on the prototype showed that the maximum amplitude of strain occurs at the SNL operating point during the synchronization of the generator with the electrical network, as well as during the shutdown phase (see Figure 1). The frequency associated with this phenomenon is close to 49 * f n (with f n the runner frequency).

Governing Equations
The incompressible flow is modeled by the URANS equations that are expressed in Cartesian coordinates as follows: with: • ρ the fluid density. • p the mean pressure. • u i the mean velocity vector. • σ ij the viscous stress tensor computed using a Newtonian linear constitutive equation. • τ ij = −ρu i u j the Reynolds stress tensor (an instantaneous variable a is decomposed following the Reynolds decomposition in a mean part a and a fluctuating part a ).
Two types of closure models are considered to compute the Reynolds stress tensor: models based on Boussinesq's assumption and models based on an explicit algebraic formulation.
Boussinesq's assumption considers a closure formulation of the Reynolds stress tensor (see Equation (4)) by analogy with the Newtonian linear constitutive equation. This assumption reduces the six unknown components of the Reynolds stress tensor to one scalar µ t , classically called the eddy viscosity. In order to compute the eddy viscosity, various approaches based on an algebraic formulation or one-or two-equation models have been developed [41].
The explicit algebraic models require formulating the Reynolds stresses as a function of the anisotropy tensor a ij (see Equation (5)). In the general 3D case, the anisotropy tensor can be expressed as the combination of ten tensors and five invariants [42].
In the present paper, three models based on Boussinesq's assumption have been considered. Each of them requires solving two additional equations in order to compute the eddy viscosity. The first one is the SST k − ω model proposed by Menter [43]. This model solved a transport equation for the turbulent kinetic energy k (Equation (6)) and a transport equation for the turbulent eddy frequency ω (Equation (7)), which is in fact a blending between the ω-and -equations (in order to limit the paper length, the values of the coefficients used in the turbulence models are not given in the paper, except the case when the value differs from the standard one. The readers are kindly requested to refer to the cited papers.).

∂ρk ∂t
The eddy viscosity is then computed as: The formulation of the SST model can include: • A production limiter for the turbulence production term P k (see Equation (9)), which avoids the build-up of turbulent kinetic energy in stagnation regions [44].
• A rotation-curvature correction term, which allows the model to be sensitive to streamline curvature and system rotation [45].
For ten years, an advanced formulation of the SST model has been developed by Menter and Egorov [30,46,47]. This advanced formulation, named Scale-Adaptive Simulation (SAS), is based on a revisited version of the k − kl model of Rotta [48], which makes an additional production term Q SAS (see Equation (10)) appear in the transport equation of ω (see Equation (7)). This term makes the von Karman length scale L vK appear, which allows the model to adjust to the already resolved scales in the flow field.
with ξ 2 = 3.51, C SAS = 2, and σ Φ = 2/3. For non-equilibrium turbulent flow with massive separations, the group of M. Braza have reconsidered the value of the coefficient C µ in the formulation of two-equation isotropic eddy-viscosity models [49,50] and anisotropic eddy-viscosity models [51]. They have shown that the value of the coefficient C µ = 0.09 derived from an equilibrium assumption is overestimated for unsteady detached flow. They proposed in the framework of the OES models a value of C µ = 0.02. In the present work, the OES modeling is applied by changing the value of the coefficient β * (which corresponds to the coefficient C µ in the k − two-equation model formulation) from 0.09 to 0.02.
The EARSM model used in the present paper is the one developed by Wallin [52] and adjusted to the BSL k − ω model [53]. The differences are the use of a value of 1.245 for A 1 instead of 1.2, the neglect of the fourth order tensor polynomial contribution β 9 , the change of the tensor basis according to [54], and the use of the standard eddy viscosity instead of the effective eddy viscosity to compute the diffusion term in the k and ω transport equations. Furthermore, the correction of the diffusion term is not used in the present calculation.
For all the models considered in the study, a wall law is used to compute the turbulent stresses in the first cell [55].

Mesh
Two computational domains have been considered taking into account the whole turbine from the spiral case to the standard draft tube (see Figure 2 on the left) or an artificially-extended draft tube (see Figure 2 on the right). Two types of mesh have been considered: two structured hexahedral meshes (see Figure 3 at the top) and one unstructured mesh (see Figure 3 at the bottom) with tetrahedral and prism elements. All the meshes have been built with the ICEM CFD v17.2 software [53]. Considering the hexahedral meshes, five sub-domains are meshed separately using a blocking method. These five sub-domains are the spiral case, the stay vanes, the guide vanes, the runner, and the draft tube. The two structured meshes differ only by the height of the first cell layer at the walls in the guide vanes, runner, and draft tube domains. On the contrary, for the unstructured mesh, the whole turbine is meshed including only three fluid domains and two internal interfaces: one between the guide vanes and the runner and another between the runner and the draft tube, in order to delimit the fixed and the rotating part of the computational domain. Along the solid walls, five prism layers are generated to better resolve the boundary layer. The number of nodes and elements for each mesh are shown in Table 1. The minimum angle and the maximum aspect ratio for each mesh and each sub-domain are provided in Table 2. For the structured meshes, the minimum angle in the runner domain is low compared to the recommended value of 20 degrees; however, the percentage of the mesh with a poor minimum angle represents less than 1% of the runner mesh. The average y + on several solid walls for the BEP is given in Table 3. It is obvious that the y + in the runner for the structured mesh is high compared to some standard requirements. This point will be further discussed in the Results Section.

Numerical Setup
The flow inside the turbine is computed using the Ansys R CFX v17.2 software [53]. The CFX solver uses an element-based finite volume method with a co-located grid. The Rhie-Chow discretization is used for the pressure-velocity coupling [56]. The system of equations is solved using a coupled formulation and a multigrid accelerated Incomplete Lower Upper (ILU) factorization technique. For the steady state simulations, a pseudo-time step is used to reach convergence, whereas for the unsteady simulations, the second order backward Euler scheme has been chosen. Regarding the convection terms, a high order scheme [57] is used for the mean transport equation and a first order scheme for the turbulence transport equations, except for the simulations with the SAS SST k − ω model, for which a high order scheme is specified. The fluxes at the interfaces between the stationary and rotating domains are computed using a General Grid Interface (GGI) algorithm. A frozen formulation is used for the steady simulations, whereas a fully-transient approach is used for the unsteady simulations.
Regarding the boundary conditions: • no-slip walls are specified at the solid walls. • the mass flow rate or the total pressure is set at the inlet of the spiral case. • the surface averaged pressure is set at the outlet through an opening boundary condition. • the rotational speed is imposed in the rotating domain.

BEP Results
The simulations at the BEP aim at checking and validating the CFD setup using only the SST k − ω turbulence model. Several simulations (see Table 4) have been performed in order to investigate the influence of the mesh, the size of the computational domain, the algorithm, and the inlet boundary conditions. Regarding the unsteady simulation CFD6, the time step was set to a 0.9 • revolution per time step, and three runner revolutions have been computed. The results compare the discharge Q ED , the speed n ED , and power P ED factors and the dimensionless efficiency (the dimensionless efficiency was computed as the ratio between the CFD and measured efficiencies) η * with the measured values (see Table 5). Regarding the discharge factor Q ED , the discrepancy between the CFD and the measurements was between +1.4% for the computation CFD5 with the imposed total pressure at the inlet and +3.8% for the computation CFD2 using the refined structured mesh. For the speed factor n ED , the discrepancy was between −1.5% for the computation CFD4 with the extended draft tube and +1.5% for the computation CFD5 with the imposed total pressure at the inlet. The power factor P ED was overestimated at least by +1.8% and at most by +4%. Finally, the dimensionless efficiency was predicted with 1% of precision. Keeping in mind the uncertainties on the input and measured values and the fact that the leakage flow through the labyrinths is neglected in the CFD, the accuracy of the numerical results was good and not strongly influenced by the numerical setup or the mesh.  Figure 4 shows the dimensionless pressure and axial velocity contours in the meridional plane provided by the simulations CFD1 (structured mesh) and CFD2 (refined structured mesh). No differences were observed on the pressure contours between the two meshes. Regarding the axial velocity contours, both simulations gave the same pattern, with spots of high magnitude at the trailing edge of the runner blades close to the shroud. The small differences in the draft tube cone can be related to the use of a frozen interface. The y + contours on one of the blades are displayed in Figure 5 for both simulations CFD1 and CFD2. The refined mesh led to a reduction of the y + by a factor of three. The y + were very high in two regions: one on the pressure side at the trailing edge close to the shroud and one on the suction side at the leading edge close to the shroud. In order to investigate the influence of the y + values on the flow in the turbine, the contours of the dimensionless magnitude of the relative velocity, pressure, eddy viscosity, and wall shear stress on the suction side of one blade are shown in Figure 6 for both meshes. For each quantity, the pattern of the contours was the same. The differences between the two meshes were observed mainly for the dimensionless eddy viscosity, which is lower on the refined mesh, as expected. The wall shear stress was slightly higher on the refined mesh at the shroud close to the leading edge. The relative velocity differed along the shroud with a lower value predicted by the refined mesh, which is in agreement with the refinement of the mesh at the walls. Finally, the contours of the dimensionless magnitude of the vorticity and eddy viscosity in a cross-section of a blade channel are displayed in Figure 7. The pattern for the two quantities was rather the same for the two meshes, even if close to the hub, the simulation CFD1 made a region with a low magnitude of the vorticity appear, which is not observed on the refined mesh (simulation CFD2). From the above analysis of the flow, it can be concluded that the simulation CFD1, with a slightly coarse mesh at the walls, allowed capturing the main features of the flow with an acceptable accuracy. Consequently, this mesh was used for the computations at the SNL operating point.

SNL Results
The SNL operating point was computed using the different turbulence models described in the governing equations subsection. The procedure for the computations was as follows: first, the SST model was used to compute a steady state simulation, then 6.25 runner revolutions with a time step set to 0.9 • of revolution per time step were computed, and finally, 6.25 additional runner revolutions with a time step reduced to 0.36 • of revolution per time step were computed. This last result was used as an initial guess to compute 6.25 additional runner revolutions simulated using: the SST model with the activation of the Curvature Correction and the Production limiter (SST-CC-Plim), the SST-OES model with a value of β * = 0.02, the SST-SAS model, and the BSL-EARSM model. The CPU time required for 6.25 runner revolutions was around 5 × 10 7 s, which corresponds to approximately one month on 24 CPUs.
The time history of the discharge and speed factors are displayed in Figure 8. A low frequency was observed on both factors for each turbulence model except for the SST-OES model, which seemed to converge to a stable operating point. The four other turbulence models presented a similar time evolution of the global parameters of the turbine. The time-averaged values of the Q ED and n ED over 6.25 runner revolutions are given in Table 6 for each simulation. However, for the SST-OES model, the final values of the Q ED = 0.0168 and n ED = 0.269 were more valuable. Comparing with the experimental data, the SST-OES model was the one that provided the closest values of the Q ED and n ED within an accuracy of 1%. The others turbulence models overestimated the Q ED and n ED , whereas, for the BSL-EARSM model, a maximum overestimation of the Q ED value by 4.7% and of the n ED value by 5.2%.  The contours of the dimensionless axial velocity ( Figure 9) and pressure ( Figure 10) on the meridional plane did not show any difference between the turbulence models. Comparing to the BEP, the pressure field was not strongly affected, whereas the axial flow was mainly concentrated close to the wall of the draft tube with a recirculation in the center of the draft tube around the tripod.
Despite the weak influence of the turbulence model on the mean pressure and axial velocity, the turbulence field was strongly affected by the turbulence model, as shown in Figure 11, which displays the instantaneous contours of the dimensionless eddy viscosity on the meridional plane. Compared to the SST model, the BSL-EARSM and the SST-OES models led to an increase of the eddy viscosity in the runner and the draft tube, whereas the SST-CC-Plim and the SST-SAS models led to a decrease of the eddy viscosity. The large magnitude of the eddy viscosity predicted by the SST-OES model was the explanation for the stable behavior observed on the Q ED and n ED factors. Several probes have been located in one guide vane and blade channel (see Figure 13). The probe GV1 was located upstream of the passage reduction due to the low opening of the guide vanes, whereas the probe GV3 was located just downstream of the passage reduction. The probes Ru1 and Ru19 were located in the middle of the blade channel, whereas the probes Ru26 (close to the trailing edge of the blade) and Ru35 (close to the hub) were located on the suction side downstream of the large spanwise vortex described previously. The time history of the dimensionless pressure and velocity magnitude in the guide vane channel is shown on Figure 14. At position GV1, the turbulence model influenced the pressure level, which was higher for the SST-OES model and lower for the SST-SAS and BSL-EARSM models compared with the SST model. The SST-CC-Plim model did not change the pressure level compared to the SST model. At position GV3, the pressure level was the same for all the models, even if the fluctuations were damped by the SST-OES model. Therefore, the pressure drop at the guide vane passage differed between the different models, explaining the differences in the predicted speed and discharge factors. Indeed, since the SST-OES model predicted a higher pressure drop at the guide vanes, then the predicted specific energy was higher, and the speed and discharge factors were lower. Regarding the velocity magnitude, at position GV1, all the models provided the same behavior characterized by a positive slope compared to the SST model. This feature was absent at position GV3, where large velocity fluctuations were observed compared to position GV1, except for the SST-OES model.
In the runner blade channel (see Figure 15), at location Ru1, the models predicted more or less the same dynamics, except the SST-OES model, which damped the pressure and velocity fluctuations. This feature was more pronounced downstream at the locations Ru19. At this position, the SST-SAS showed larger fluctuations than the others turbulence models.
Finally, at locations Ru26 and Ru35 (see Figure 16), the SST-OES once again damped the fluctuations and provided a stable signal. The SST-CC-Plim model followed the behavior of the SST model. The BSL-EARSM model presented a similar behavior, but with a shift regarding the time-averaged pressure and velocity. The SST-SAS model predicted a higher amplitude of the pressure and velocity signals. The pressure signals at the locations Ru26 and Ru35 have been transformed into the Fourier space, and the corresponding dimensionless amplitude is plotted against the scale frequency in Figure 17. At location Ru26, the SST model showed several peaks at frequencies lower than the scale frequency f / f n = 24, which corresponds to the guide vane passage frequency. Above this frequency, the amplitude decreased by one order of magnitude. Apart from the SST-SAS model, all the others models exhibited a peak at the scale frequency f / f n = 24. Beyond this frequency, the decrease in amplitude disappeared, which means that contrary to the SST model, these models increased the fluctuations at frequencies higher than the guide vane passage frequency. In the particular case of the SST-SAS model, no specific frequency was put in evidence due to the higher amplitudes predicted for all the frequencies, masking the guide vane passage frequency. At location Ru35, no specific frequency was captured by the models; however, SST-SAS still provided higher amplitudes for all the frequencies compared to the other models.
Regarding the large strain fluctuations that were observed experimentally at a scale frequency f / f n = 49, none of the simulations made such a frequency appear. Therefore, it seemed that the flow did not provide a specific excitation source for the structure. Nevertheless, as shown by the pressure spectrum provided by the SST-SAS model, the level of pressure fluctuations was increased around the scale frequency f / f n = 49 by the advanced turbulence models compared to the SST model.

Conclusions
The SNL operating point of a high head Francis turbine was investigated using several "advanced" RANS turbulence models. The unsteady simulations have been compared with the results obtained using the standard SST k − ω turbulence model and the available experimental data. The mean flow was slightly influenced by the turbulence model, which means that the flow pattern and the time-averaged global parameters of the turbine such as the speed or discharge factors were computed within an accuracy of less than 5%. This range of 5% was due to a difference in the pressure drop at the restriction passage imposed by the small guide vane opening. On the contrary, the turbulence field and the dynamics of the flow were strongly influenced by the turbulence models. For instance, the SST-OES model led to a quasi-stable behavior of the flow due to an increase of the eddy viscosity in the runner and the draft tube. Conversely, the SST-SAS model increased the level of pressure and velocity fluctuations compared with the standard SST model. The BSL EARSM model and the SST model with the curvature correction and production limiter shared the same behavior as the SST model, but they provided a little bit higher pressure and velocity fluctuations, mainly for frequencies beyond the guide vane passage frequency.
However, none of the models showed a specific excitation at the scale frequency f / f n = 49, whereas experimentally, the strain gauges exhibited a strong excitation. As this frequency could correspond to an eigenfrequency of the turbine, further investigations are required. It could be already suggested that the SST-SAS should provide a higher excitation source due to the higher pressure fluctuations predicted.

Conflicts of Interest:
The authors declare no conflict of interest.