Non-Linear Harmonic flow simulations of a High-Head Francis Turbine test case

This work investigates the use of the non-linear harmonic (NLH) method for a high- head Francis turbine, the Francis99 workshop test case. The NLH method relies on a Fourier decomposition of the unsteady flow components in harmonics of Blade Passing Frequencies (BPF), which are the fundamentals of the periodic disturbances generated by the adjacent blade rows. The unsteady flow solution is obtained by marching in pseudo-time to a steady-state solution of the transport equations associated with the time-mean, the BPFs and their harmonics. Thanks to this transposition into frequency domain, meshing only one blade channel is sufficient, like for a steady flow simulation. Notable benefits in terms of computing costs and engineering time can therefore be obtained compared to classical time marching approach using sliding grid techniques. The method has been applied for three operating points of the Francis99 workshop high-head Francis turbine. Steady and NLH flow simulations have been carried out for these configurations. Impact of the grid size and near-wall refinement is analysed on all operating points for steady simulations and for Best Efficiency Point (BEP) for NLH simulations. Then, NLH results for a selected grid size are compared for the three different operating points, reproducing the tendencies observed in the experiment.


Introduction
Over the last years, the interest for unsteady CFD simulations in hydraulic machinery has been increasing for various motivations. One important objective is to accurately predict Rotor/Stator interactions, as hydro turbine and pump designers need to reduce the pressure fluctuations on hydraulic surfaces to minimum values. The final goal is to improve the behavior of their machines regarding to vibration damages.
For this purpose, NUMECA has developed an efficient modelling approach, the Non-Linear Harmonic (NLH) method, into his commercial software packages FINE™/Turbo and FINE™/OPEN. Compared to more "classical" time-marching unsteady simulations, the NLH method presents various advantages. Because of the transposition to the frequency domain, only one blade channel is required like for a steady flow simulation. As only the perturbations generated by the adjacent blade rows are computed, the computing cost is largely reduced, sometimes up to 2 orders of magnitude. The method also features an improved treatment that enhances the flow continuity across the Rotor/Stator (R/S) interfaces by a reconstruction of the harmonics and the time-averaged flow on both sides of the interface.
Since its introduction into FINE™/Turbo in 2006, the NLH method has been applied and validated against experimental data and classical Unsteady Reynolds Averaged Navier-Stokes simulations on a wide range of academic and industrial application, mostly for compressible applications [1] but also for incompressible test cases [2,3,4,5]. A specific validation for incompressible and low Mach number flows has been presented in [2]. The method has also been successfully compared to experimental and simulation results from a classical time-marching approach on a hydraulic turbine configuration in [4]. The Francis99 workshop [6] offered a new opportunity to validate and present the benefits of NLH method for hydraulic machinery. In the work presented here, we applied the method to the workshop test case, with specific interest in studying the consistency of the method regarding to grid size and nearwall resolution. Some complementary results for different operating points (best efficiency, partial and high load) are also presented and discussed for a selected grid configuration.

Test case: Francis99 turbine
The Francis-99 test case is a high-head Francis turbine model investigated at the Water Power Laboratory at the Norwegian University of Science and Technology (NTNU) in Norway. The turbine features a spiral casing, 14 stay vanes (SV), 28 guide vanes (GV), a runner (RU) with 15 blades and 15 splitter blades and an elbow-type draft tube. More details and the geometry of the complete turbine are available in [6].
Experimental data are available for 3 operating points. The complete set of experimental data features time dependent pressure measurements and velocity profiles in the draft tube as well as global performances. In this work, we have been comparing simulations results with all available pressure measurements provided (except draft tube), at locations illustrated in Figure 1: one in the vaneless space between the guide vanes and the runner blades (VL01), two on a blade pressure side (P42 and P71), one on the suction side of the adjacent blade (S51). Operating points considered are Best Efficiency Point (BEP), Partial load (PART) and High Load (HIGH). The conditions for the operating points are defined in table 1; for partial load, they differ from the proposed conditions of workshop 1 and actually correspond to conditions for workshop 2.

Numerical model and Non-Linear Harmonic method
The commercial package FINE TM /Turbo is used for this study. The flow solver is a three-dimensional, density-based, structured, multi-block Navier-Stokes code using a finite volume method. Centraldifference space discretization is employed for the spatial discretization with Jameson type artificial dissipation. A four-stage explicit Runge-Kutta scheme is applied for the temporal discretization. Multigrid method, local time-stepping and implicit residual smoothing are used in order to speed-up the convergence.
The base NLH, presented by Vilmin et al. [7], relies on the nonlinear harmonic introduced by He and Ning [8]. The instantaneous flow in a blade row channel is decomposed into a time-averaged value and a sum of periodic unsteady perturbations, those being provoked by the blade-passing frequencies of the adjacent blade rows. Thus, the unsteady conservative flow variable U reads at location r and time t: Each unsteady perturbation U' can be decomposed into a Fourier series with a finite number NH of time harmonics where ̌ and ωk are the complex amplitude and the frequency of the time harmonic k, respectively; ̌− is its complex conjugate.
The complex amplitudes ̌ are solved by marching in pseudo-time the first-order linearized Navier-Stokes equations cast in the frequency domain. Following the finite-volume formulation, it reads at each control volume i of volume Ωi where ̌, ̌ are the complex amplitudes of the linearized convective and viscous fluxes; ̌ is the complex amplitude of the linearized Coriolis and centrifugal terms, plus the frequency source term -iωǩ.
The frequency casting allows for keeping only one interblade channel per blade row. Thus, the connection of the time harmonics at the periodics simply follows the phase lag boundary condition characterized by the interblade phase angle.
The time-mean flow is solved by marching in pseudo-time the time-averaged Navier-Stokes equations, formulated following the finite-volume method Like Reynolds averaging, it includes additional stresses, under the form of time-averages of products of fluctuations. This contribution of the periodic unsteadiness to the time-mean flow is the so-called deterministic stress and represents non-linear effects in the time-mean equations. It is easily calculated from the in-phase and out-of-phase components of the solved harmonics, thus achieving the closure of the system.
For the NLH solver, time-mean and harmonic equations are solved using the spatial and temporal discretization schemes previously described. In order to avoid a too strong artificial dissipation and a too low convergence rate for low-speed flow applications in a density-based code, preconditioning has been applied to equations (3) and (4). Therefore, the left hand side of both equations is multiplied by a preconditioning matrix. The preconditioning matrix relies on a pseudo-sound speed that depends on the local velocity, the viscosity, and the pressure gradient. It is, therefore, just a function of local variables and does not depends on global parameters. Such a method allows to be able to solve purely incompressible fluid in a density-based flow solver and to have a convergence rate independent of the Mach number. More details on the preconditioning method that has been retained can be found in [2].

Computational Meshes
Various grids were generated for this work using NUMECA software AutoGrid5™. All grids are structured, featuring only hexahedral type cells, and are fully matching inside the mesh domain but also at periodic boundaries. For the stator domain, a passage mesh of 1 Stay Vane and 2 Guide Vanes was defined. For the rotor domain, a passage mesh including 1 main blade and 1 splitter blade was generated, till the end of the draft tube cone. Complex geometrical features were taken into account: guide vane earrings and volume under the runner cap till the axis. Some views of the mesh generated are shown by figures 2 and 3.  To analyse the grid size influence as well as the near-wall resolution, 3 grid points' distribution were defined, for the 3 operating points considered, resulting in a total of 9 grids used for computations. The different grids can be classified as follows:  Coarse High-Reynolds mesh (y + values > 20): "High Re".  Intermediate Low-Reynolds mesh (y + values < 10): "Low Re".  Fine Low-Reynolds mesh (y + values < 5): "Low Re fine". Characteristics of these meshes are given in table 2 and 3, and some comparison of mesh distribution are shown by figure 4. It can be seen from values in table 3 that the mesh quality is high-level and very uniform for all operating points and guide vane opening. As complementary characterization of all grids, average angles values are always over 75º, for and average expansion ratio under 1.2.  The strategy followed for refinement was: 1. Mesh was defined for High Re ensuring best possible grid quality from a first cell size, adapted to the use of a High Reynolds turbulence model (i.e. with wall functions). 2. From "High Re" to "Low Re": refinement of boundary layer mesh for use of Low Reynolds turbulence model (i.e. resolve the entire boundary layer, including the viscous sublayer). As an important number of cells are located in the boundary layer mesh, this modification leads to a more than 50% increase of cells count. 3. From "Low Re" to "Low Re fine": doubling the number of cells, by increasing the number of cells by about 25% in all directions (I, J and K as mesh is structured).

Figure 4.
Comparison of Blade to blade meshes in the inter-blade spacing between guide vane and runner: "High Re" (topleft), "Low Re" (top-right) and "Low Re fine" (bottom-left). A grid refinement can be seen upstream runner leading edge, due to the capture of seal leakage geometry at hub and shroud.

Simulations
First, standard steady computations have been run using mixing plane approach for all operating points and all grid types. For boundary conditions, inlet flow rate and rotation speed are defined as per table 1 and an outlet average static pressure is set for all computations. Inlet flow direction and turbulent variables were determined based on precursor spiral casing and distributor computations. Pitch-averaged flow variables are imposed at inlet surface of Stay-Vane + Guide-Vane + Runner configuration presented here. These steady computations were run for two different turbulence models, Spalart-Allmaras (SA) [9] or k-Omega « Shear Stress Transport » (SST) [10], both with extended wall-functions (i.e. wall functions are activated or not depending on the local y + values). Very little difference were observed between the two turbulence models. Indeed, the maximum discrepancy observed for the global torque is no more than 10 N.m on the fine grid (4 N.m on the intermediate grid and 2 N.m on the coarse grid). For the total pressure difference between the inlet and the outlet of the computational domain, the difference between turbulence models are at maximum 3,500 Pa, for finest grid, PART operating point. Those numbers have to be compared to the torque and the total pressure difference displayed on figure  5 and 6. Therefore only Spalart-Allmaras turbulence model was retained for the remaining NLH simulations of this study. Those harmonic simulations were run for BEP operating condition on the three different meshes already presented above. For the two others operating points, i.e. PART and HIGH, only the High Re mesh was used.
The computed frequencies are shown in Table 4 and are defined using the formulas 5 and 6. It should be mentioned that first frequency calculated is based on a periodicity defined by 2 blades (2 guide vanes or 1 main blade and 1 splitter blade), therefore the main effects are occurring for second and fourth computed frequencies that correspond to the blade passing frequency and its first harmonic.
Ω is the rotation speed in RPM and NbS number of stator blade, NbR number of rotor blades. All computations were analyzed after reaching full convergence: stable residuals (steady and/or harmonics), global values (for example torque) and local values (average static pressure at probe locations, real and imaginary parts of harmonic pressure).

Results
Results are presented in two sections. First the grid size impact is analyzed both for steady and NLH simulations based on global and local outputs from the various CFD simulations. Then, the results from the three operating conditions are compared focusing on Harmonic results in order to investigate the impact of Rotor/Stator interactions.

Grid size influence
To analyze the grid influence, steady computations were run for all 9 meshes and NLH computations for BEP conditions only.
Results are first compared looking at global values directly obtained from the solver for the 9 steady computations (figures 5 and 6): runner torque and total pressure difference between inlet and outlet. A typical convergent behavior of these quantities versus grid size is observed, and tendencies are well respected for all grids and all operating points. Maximum computed variations in torque and total pressure difference between the two finest meshes are observed for the partial load configuration and do not exceed 1.8 % and 1.2 %, respectively.  The sensitivity of local variables to mesh density is performed based on results from NLH simulations at Best Efficiency Point. Time-averaged static pressure at sensors locations together with amplitude of pressure fluctuations at the same places are presented on figures 7 and 8.
All the simulations exhibits almost the same time-mean static pressure and are able to produce a pressure distribution that is in close agreement with experimental data. Maximum difference in static pressure is 2.4 % and is observed for simulation on the High Re. Average difference for all grids and locations is 0.5 %.
On figure 8, it appears that when refining the mesh, the amplitude of the pressure fluctuations at the BPF are getting closer to experimental values for all the sensors. The major difference is observed for the most upstream sensor (VL01). At this location strong azimuthal variations of the amplitude of the pressure fluctuations at the BPF are observed along the turbine shroud in the vaneless space between the guide vane and the runner (figure 9). Indeed, just downstream of the guide vane trailing edge, moving the location of the numerical sensor by a few degrees (corresponding to a few grid cells) could leads to variations of amplitude of pressure fluctuations of up to several hundreds of Pascal. Nonetheless, whatever mesh is used, NLH results are able to reproduce the observed tendencies with higher pressure fluctuations in the middle of the runner blade pressure side, lower variations downstream and at the sensor location along the suction side.

NLH simulations results for the 3 operating points
According to the above mesh sensitivity analysis, it is assumed that trustable tendencies can be obtained using the "high Re" mesh type, i.e. the coarser one. It offers a more stable behaviour and lowest computing time for NLH simulations and thus should be more suitable for a design purpose. Therefore, all the remaining simulations at off design condition have been carried-out on this mesh. Simulations were run for the three operating points and are compared to experimental values for conditions BEP and HIGH. Measured amplitude of pressure fluctuations at the various sensors locations are available for those two operating points and are displayed on figures 10 and 11. Computed amplitude of pressure fluctuations at the same locations are also displayed on figure 12 and 13. From the comparison between experimental and numerical results it appears that the NLH method is able to accurately predict the observed tendencies.
At BPF, a strong increase of about 700 Pa in pressure amplitude is observed for the upstream most sensor VL01. Even if it is of a lesser extend (400 Pa) such a tendency is also observed in numerical results. This is also linked to a measured increase in the pressure fluctuations at the same location but at a frequency corresponding to half of the BPF. This can be attributed to the effect of a distinct flow pattern along the main runner and the splitter blades at higher flow rate. For all the sensors located along the runner blade a slight decrease in pressure fluctuation is measured at the blade passing frequency. CFD results are also in agreement with such a finding.
At the frequency corresponding to 2×BPF, the measured pressure fluctuations cannot be easily distinguished from the background noise. Nonetheless, both measurement and CFD results exhibit an increase in pressure fluctuations for VL01 sensor at such a high frequency. It should also be noted that the flow solver does not produce enough pressure fluctuations along the runner for 2×BPF as can be identified in figure 10. As already mentioned, strong azimuthal variations of pressure fluctuations are observed in the vanelees space around VL01 either for the BPF or twice the BPF. This can be identified on figures 14 to 17. Such a pattern has also been observed either numerically or experimentally for a centrifugal pump in [2]. Therefore, it should be of interest to have access to some measurements that can provide values for an array of sensors at these locations to further validate simulations results.

Conclusions
Steady and harmonic simulations have been carried out successfully on the Francis99 High-Head Francis Turbine test case. High quality structured grids have been generated for three operating points (three GV openings) and three grids were generated for each operating point, with increased mesh size and near-wall resolution. Steady computations have been run for all grids, and harmonic computations have been run for the three operating points on coarsest grid after ensuring consistency of the results on this mesh. Analysis of the steady and harmonic results focused primarily on grid size influence, for global and local values. Tendencies are well respected for all grid sizes. Comparisons with experimental data has also been performed and bring a high value of confidence in the numerical results.
On the overall, excellent consistency and stability of the modelling has been demonstrated and Hydraulic R&D/Design engineer can safely choose some parameters of their CFD harmonic computations, for example the grid size, to increase robustness and reduce even more computing time of their Rotor/Stator unsteady analysis.