Francis-99 turbine numerical flow simulation of steady state operation using RANS and RANS/LES turbulence model

We performed numerical simulation of flow in a laboratory model of a Francis hydroturbine at three regimes, using two eddy-viscosity- (EVM) and a Reynolds stress (RSM) RANS models (realizable k-ϵ, k-ω SST, LRR) and detached-eddy-simulations (DES), as well as large-eddy simulations (LES). Comparison of calculation results with the experimental data was carried out. Unlike the linear EVMs, the RSM, DES, and LES reproduced well the mean velocity components, and pressure pulsations in the diffusor draft tube. Despite relatively coarse meshes and insufficient resolution of the near-wall region, LES, DES also reproduced well the intrinsic flow unsteadiness and the dominant flow structures and the associated pressure pulsations in the draft tube.


Introduction
An important task of hydropower plant is the regulation of power in the energy system. During load changing, hydraulic units repeatedly undergo off-design modes of operation. In these modes, the flow remains essentially swirling after passing through the water turbine runner. The instability of swirling flow results in the appearance of intense low-frequency hydrodynamic pulsations, which threaten the safety and reliability of turbine structure.
The improvement of hydraulic machines efficiency and failure tolerance is impossible without studying the physical mechanisms of hydrodynamic processes, among which an important role is played by the transient phenomena associated with the formation of large-scale vortex structures. One of the mechanisms responsible for generating flow pulsations is the precessing vortex rope. It is formed downstream the runner under part-load or over-load conditions when the flow has a large residual swirl after passing through the runner [1][2][3][4][5].
The precession of vortex rope can be a serious danger for hydraulic equipment due to powerful flow pulsations that lead to strong vibrations of turbine structure. It the case of resonance, it can result in the turbine structural failure. Pressure pulsations generated by precessing vortex rope may also affect the cavitation process and enhance cavitation erosion. To predict the resonance phenomena and find the ways of suppressing instability, one needs the detailed information about the characteristics of pulsation regimes and flow structure. It should be noted that the working approaches must meet the requirement of minimizing energy losses (increasing the efficiency of water turbine). It can be realized The precession of vortex rope has been studied for quite a long time. These issues have long been in focus of the experimental research complemented by some simplified analytics, but a full account of the three-dimensional time dynamics has remained to a large extent beyond the reach of even most advanced laser-based measurements and diagnostics techniques. Computer simulations have been seen as a potential tool for capturing the subtleties of the development of instabilities and associated vortical and turbulence structures. Such information can be obtained by direct and large-eddy simulations (DNS, LES) over a range of important scales. However, the accuracy of DNS and LES depends much on the numerical resolution. For complex configurations the Reynolds-averaged Navier-Stokes (RANS) approach, being computationally less demanding, has been considered as a more rational option for industrial purposes. However, because if its empirical contents, the RANS have not been regarded as a trustful research tool, and moreover, their credibility depends much on the choice of the turbulence model to close the time-or ensemble-averaged conservation equations.
Despite the extensive literature (e.g. Doerfler [2], Alligne et al. [6], comprehensive experimental data for hydroturbines are still scarce, apart from mean velocity and second-moment turbulence statistics measured at selected locations in laboratory models of some popular turbine types [7][8]. On the other hand, there are still many open questions regarding the numerical simulation of flows in turbines of realistic configuration and practical relevance. Because of typically very high Re numbers, DNS and even LES, requiring enormous numerical mesh and computer resources, are usually beyond the reach of designers, but also of many researchers, leaving the URANS or hybrid RANS/LES as the only viable routes. Here the first challenge is related to selecting the appropriate RANS model that can reproduce satisfactorily the unsteady flow and the dynamics of at least the dominant large-scale vortex structures subjected to strong swirls, curvature and adverse pressure gradients in complex geometric shapes. The popular linear eddy-viscosity modes (k-ε, k-ω, and their many "improved" variants) lack on sufficient physics and have shown to perform poorly in describing such flows.
We report here the results of simulation of flow in a Francis-99 hydraulic turbine at three regimes using several RANS models, a DES (detached-eddy simulations) and LES. The models were benchmarked against the experimental data. We discuss the flow patterns and vortical structure in the turbine draft tube and its impact on flow stability and pressure pulsations, especially at part loads characterized by unsteady helix ropes.

Turbulence models
Flows in hydroturbine systems are characterized by several features which pose a challenge to the RANS models. The flows in the guide vanes and the rotating impeller are dominated by the bounding wall -blades and housing, and a strong favourable pressure gradient which for the design operating conditions (no flow separation) can be handled reasonably well with a well-tuned wall-integration (WIN) model that (for the runner) accounts properly for the rotation effects. Accurate predictions of the blade-tip leakage, the follow-up vortical wakes and flow separation on blades at part loads, require more advanced, anisotropy-accounting RANS models. The phenomena in a draft tube are different. Here the flow is governed by a strong swirl and a bluff-body recirculation zone behind the impeller hub, which is usually reinforced by the negative pressure due to the swirl. The 90 degree bend encountered in elbow draft tubes, with a subsequent circular-to-rectangular change of the cross-section further modify the flow especially at part loads creating very complex vortical patterns with additional recirculation and secondary flows. On the other hand, the walls are smooth, and despite a strong elbow curvature, the wall-adjacent azimuthal fluid velocity and turbulence seem to follow reasonably well the common wall scaling so that computations seem manageable by the common wall-integration schemes or, according to some publications, even with using the standard wall-function approach. The apparent success of relatively standard wall treatment could in part be explained by strong radial pressure gradients and by inviscid interactions that diminish the wall effects as well as the role of the employed turbulence model on the dominant vortical structures and the consequent mean flow pattern. The above brief overview of some major phenomena encountered in hydromachinery illustrate the importance of choosing the proper RANS model, which depends on which phenomena are in the focus of the computations. As stated above, our focus is on the draft tube and, specifically, on the strong and coherent swirling structures, primarily the precessing vortex cores and especially on the undesirable twin-ropes and the associated strong pressure pulsations at part loads. The linear eddy-viscosity (LEVM) turbulence models, widely spread in engineering calculations, are known to describe poorly most features listed above. Nevertheless, we tested some popular models, all validated against a wellresolved dynamic LES, aimed at scrutinizing the ability of various models to capture the above mentioned specific structures and pressure oscillation in the considered draft-tube model.
The computations here reported were performed with the ANSYS-FLUENT code using the built-in two popular eddy-viscosity (k-ω SST and the realizable k-ε), a Re-stress model (Launder-Reece and Rodi, LRR, with the linear pressure strain and wall-echo terms), and a DES (based on the k−ω SST Menter's model). In parallel, for reference, LES with the WALE subgrid-scale viscosity model was carried out on finer grids.
The high-Re-number realizable k-ε model and the basic RSM were solved using the standard wall functions. A test with the enhanced wall treatment and the non-equilibrium wall function available in the FLUENT code produced no significant differences in the results. For the k−ω SST we used the enhanced-wall-treatment ω-equation model (EWT-ω).
A specific issue in the computer simulation of hydraulic turbines is the treatment of the rotating impeller and the rotor-stator interaction. Several approaches can be found in the literature, i.e. dynamic, sliding and moving grid methods, and those based on a moving reference frame. The latter is the most common and the simplest way to model the runner rotation. It assumes that the runner is fixed and the equations are solved in a rotating reference frame. This formulation is often referred to as the "frozen rotor" approach. In this paper, the modelling of the runner rotation was performed in the rotated reference frame for the runner zone. The obtained results are then rotated with the runner rotation speed and as such imposed as the inflow field as the inlet into the draft tube. The earlier test calculations proved that this approach is credible for describing the integral flow characteristics including the dominant flow pulsations [4, 5, 8 -10]. A comparison with the computationally more demanding method using sliding meshes showed that for this type flows with a focus on draft tube dynamics the results are almost the same.
The computations were performed using the finite-volume method on unstructured grids. The coupling of the velocity and pressure fields for incompressible flow was ensured using the SIMPLE-C procedure. For the URANS models the convective terms in all equation were discretized by a secondorder upwind scheme, whereas the second-order central difference scheme was used for LES and DES (in the LES region). The time derivatives were approximated by an implicit second-order scheme. The time step was set from the condition CFL<2, but in the bulk of the draft tube the CFL number was less than 1.0. Some peaks exceeding the value of 2.0 appeared close to the guide-vane and impeller blades. The statistics was gathered for each run during about 10 sec of real time.

Computational domain, boundary conditions and grids
The computational domain included all the elements constituting the turbine system (spiral casing, guide vane, impeller and draft tube) see (Fig.1). The 3D geometry of the turbine has been taken on the website workshop Francis (2016). Since the intake channel was not considered in the calculations, the boundary conditions were set at the inlet to the turbine spiral casing by imposing a total head measured in the experiment and taken from Table 1. As the LES method requires initial timedependent forcing mimicking turbulence, random velocity fluctuations were imposed at the inlet, generated by the method proposed in [11]. Kinematic viscosity (m 2 s -1 ) 9.57E-7 9.57E-7 9.57E-7 Gravity (m s -2 ) 9.82 9.82 9.82 The computations have been carried out using unstructured grids with the total number of 5.12 and to 14.28 mln. nodes for the whole domain (Fig.2). The average distance of the wall-nearest grid nodes normalized by the wall units, 1 y + , for the 5.12 mln. grid was 150, for the 14.28M grid 1 y + was 89. As seen, both grids are coarse grid for the LES method.
To check further the grid resolution, the characteristic cell size production, was found by Picano and Hanjalić [12] to serve well as a criterion for mesh resolution in free flows: when the grid-cell size was kept smaller than the shear scale (nominally / 1.0 s L ∆ < , relevant primarily in high-shear regions), the LES of a free round jet reproduced well the mean flow and turbulence statistics at very high Re number using a very coarse grid despite the cell sizes being two order of magnitude larger than the Kolmogorov scale. Fig. 3 shows the distribution of Δ/L K and Δ/L s in the central section of the turbine calculated by RSM model for part load (PL) regime for grid 5.12M. As can be seen, the ratio / K L ∆ of around 40 to 80. Interestingly, the cell-size to shear-scale ration / S L ∆ is below 1.0, except in the core region, and of course very close to the wall. This is in agreement with the findings in [12], which shows that the LES method can return satisfactory predictions when / S L ∆ is less than 1.0. Comparison of the calculated and experimental radial velocity profiles are shown in Fig. 7, 9, 11. As can be seen, for all regimes the discrepancy between the calculation and the experiment is very great for all turbulence models. Perhaps, this is due to the fact, that the radial velocity has lower values compared to axial and tangential velocity and, therefore, radial velocity has high relative error of measurement.

The results of simulations flow in the Francis-99 turbine
Comparison of the calculated and experimental axial velocity profiles are shown in Fig. 6, 8, 10. As can be seen, for all regimes the calculated and experimental data agree well. The worst results are for the PL regime. For this regime all computational models give the deflection of axial velocity at the center of the draft tube. The experiment has not of this drop of axial velocity. Using of more detailed grid does not improve the result (see Fig. 6a). The results of LES calculations for the rough (5.12M) and the detailed mesh (14.28M) are similar. Also it should be noted, that the difference between the different models is not very large. For BEP and HL regimes the agreement between the calculation and experiment is a lot better. The best agreement between the calculation and experiment is obtained for the BEP regime. The results of the analysis of the axial velocity profiles showed that the best agreement with experiment has an LES model. The DES model also gives good results. The RANS models are close to each other. As mentioned above, the radial velocity profiles analysis does not make sense.  Another reason of disagreement between calculation and experiment is different value of discharge in the calculation and experiment. At given head the calculated value of discharge is about 10% lower than the experimental value (see Table 2 -4). For this reason, the calculated profiles are slightly higher than the experimental. The discrepancy between experiment and simulation is significantly reduced in case of fixing the discharge on the turbine input. An example of this is shown in the Fig.8a for LES calculation.
Comparison of the calculated and experimental integral characteristics of the turbine for all regimes are shown in the Tables 2-4. As can be seen, for all regimes calculated characteristics are 7-10% lower than the experimental value, except of the efficiency, which calculated value overestimated in comparison with the experimental one. The part load operation point shows the largest discrepancy. The calculations by different turbulence models give similar results. It is impossible to say which model is best for finding the integral parameters. Practically, it is very important to know the level of the pressure pulsation in the turbine. Therefore, it is interesting to compare the calculated and experimental values of the pressure fluctuations. These data are shown in Fig. 12. The plot shows the static pressure on the draft tube wall. The pressure were received at the point DT5, (see workshop web site and Fig. 1).
Obtaining of the maximum level of pressure fluctuations is very important for safe operation of the turbine. Analysis of the results shows that the magnitude of the experimental pressure fluctuations is 1.68 kPa. The calculations are in good agreement with experiment (see Fig.12). Thus, the LES calculation give maximum of the pressure fluctuations 1.07 kPa, DES model give 0.99 kPa, and RSM model give 1.22 kPa. This is interesting, because it was previously thought that the RSM model should give the lower level of the pressure fluctuations. Thus it is shown that the RSM model can be used to estimate the level of the pressure pulsations in the turbines. The behavior of the flow structure in the turbine is of great interest for analysis of the pressure pulsation of the flow. Fig. 13 shows vortex structure of the flow downstream the runner for different unsteady models at PL regime. Flow structure is shown by an instantaneous values of iso-surfaces of pressure. We see, that in this regime vortex rope structure is formed. In this case the LES and DES calculation gives two intertwined precessing vortex. Also it is interesting that, the Reynolds stress model captures well the phenomenon, although with smoothed fine-scale structures. It should be noted that the URANS model also provides some transient pressure pulsation in the draft tube.
Interesting details can be observed regarding to the pressure pulsations and axial velocity profile downstream the turbine runner. In the modes with the high pressure pulsations, the wide areas of counterflow are visible (Fig. 6). This again proves the fact that a zone of counterflow on the axial velocity profile can be an indicator of the intensity of transient phenomena in the turbine [4,5].  Figure. 13. Vortex structure in the PL regime for different turbulence models.
In addition to the pressure pulsation it is also very important to know the frequency of the pressure pulsations. It is important for estimating the potential resonance of the turbine design. A comparison of the experimental and calculated spectra is shown in Fig. 14. The flow after the runner has precessing vortex structures (Fig. 13). For this reason, the spectrum has a dominant frequency. In the experimental data the frequency of the pressure pulsations for this regime is 1.53Hz. This frequency is the frequency of the precession of intense vortex core.
As it is obvious from the spectrum of the pressure pulsations, the calculation results are well proved by the experiment, both for the frequency and intensity of the pulsations. The LES calculation provides frequency 1.57 Hz, DES provides frequency 1.22 Hz. And the RSM calculation provides frequency 1.47 Hz. Figure. 14. The spectrum of pressure pulsations at the DT5 point for PL regime.

Conclusions
The paper presents the results of numerical simulation of a flow in a laboratory model of a Francis hydroturbine at three regimes, using two eddy-viscosity-(EVM) and a Reynolds stress (RSM) RANS models (realizable k-ε, k-ω SST, LRR) and detached-eddy-simulation (DES), as well as large-eddy simulation (LES). Three operation point proposed for the Francis-99 workshop were calculated. Comparison of calculation results with the experimental data was carried out.
The results of the analysis of the axial velocity profiles showed that the best agreement with experiment has an LES model. The DES model also gives good results. The RANS models are close to each other. The worst results are for the PL regime. For BEP and HL regimes the agreement between the calculation and experiment is a lot better. The best agreement between the calculation and experiment is obtained for the BEP regime.
The reason of disagreement between calculation and experiment may be different value of discharge in the calculation and experiment. At given head the calculated value of discharge is about 10% lower than the experimental value (see Table 2 -4). For this reason, the calculated profiles are slightly higher than the experimental. The discrepancy between experiment and simulation is significantly reduced in case of fixing the discharge on the turbine input.
The structure of the flow behind the runner was analyzed by numerical simulation, which revealed the effect of the flow structures on the frequency and intensity of the pulsations. The calculations confirmed that a vortex rope behind the runner is the main cause of low-frequency pressure pulsations in hydraulic turbines. Its dynamics determines predominantly the behavior of turbine pulsations.