Investigation of the fluid-structure interaction of a high head Francis turbine using OpenFOAM and Code_Aster

The increasing energy consumption and highly stressed power grids influence the operating conditions of turbines and pump turbines in the present situation. To provide or use energy as quick as possible, hydraulic turbines are operated more frequent and over longer periods of time in lower part load at off-design conditions. This leads to a more turbulent behavior and to higher requirements of the strength of stressed components (e.g. runner, guide or stay vanes). The modern advantages of computational capabilities regarding numerical investigations allow a precise prediction of appearing flow conditions and thereby induced strains in hydraulic machines. This paper focuses on the calculation of the unsteady pressure field of a high head Francis turbine with a specific speed of nq ≈ 24 min-1 and its impact on the structure at different operating conditions. In the first step, unsteady numerical flow simulations are performed with the open-source CFD software OpenFOAM. To obtain the appearing dynamic flow phenomena, the entire machine, consisting of the spiral casing, the stay vanes, the wicket gate, the runner and the draft tube, is taken into account. Additionally, a reduced model without the spiral casing and with a simplified inlet boundary is used. To evaluate the accuracy of the CFD simulations, operating parameters such as head and torque are compared with the results of site measurements carried out on the corresponding prototype machine. In the second part, the obtained pressure fields are used for a fluid-structure analysis with the open-source Finite Element software Code_Aster, to predict the static loads on the runner.


Introduction
The transition of the energy system all over Europe of nation states has broken up an old history of electricity supply and grid control. New suppliers based on renewable energy sources have appeared and made the electricity production more flexible. Forecasts about the gap between supply and demand became even more uncertain and future trends about energy storage technologies and their role in the energy system are under investigation and research. Different studies with various approaches cannot deliver satisfying answers to the needs of our energy system within the next 20 years and beyond. The crucial topic of energy storage is common to all studies: without short-term and long-term storage possibilities as well as decentralized storage capabilities no transition of the energy system is possible. Short-term storage possibilities are characterized by their short charging and discharging characteristic and are therefore ideal for a quick demand response and grid stabilization. In 2008 the study of the Energy Storage Task Force [1] predicted a new era of increased importance and usage for the battery technology. Today we see that these predictions are behind the expectations of 2008. Reality showed that especially the new storage technologies could not be successful due to the deterioration of the market sales price after 2009. In 2013 the roadmap of the EASE/EERA Core Working Group [2] assumed that the short term electricity balancing market will be there, where energy storage is based on commercial business cases. This means more research and development activities for recent technologies but also new challenges for the established ones. The need for research and adaption of existing hydropower plants as well as the development of decentralized pumped storage technologies are summarized in the Hydro Equipment Association roadmap [3]. The increasing amount of volatile energy sources (e.g. wind and solar) forces the operators of large pumped storage plants to deliver more regulating power to the electric grid and to decrease the response time of the hydraulic machines. Francis turbines, which can be used within a wide range of discharge and head, are therefore operated over longer periods of time at off-design conditions. This involves the appearance of unstable flow phenomena, depending on the specific speed and the operating point of the machine (see Magnoli et al. [4]). The dynamic flow effects have an important influence on the expected lifetime of the Francis turbine. In full load and higher part load operation, the rotor-stator interaction is mainly important, especially for high head machines. In lower part load, draft tube vortex ropes and interblade vortices appear, which leads to pressure pulsations in the turbine. Start-stop cycles and speed-no-load operations have a more stochastic influence on the runner and a significant impact on the life expectancy, as Coutu et al. [5] or Gagnon et al. [6] showed. Several studies have been published to examine the influence of the static and dynamic pressure effects on the structural behavior of the Francis runner by using fluid-structure interaction (FSI) (see [7][8][9][10][11]). In order to investigate the lifetime of a high head Francis turbine, an approach has been applied by Doujak et al. [12] within a large research project, using numerical flow simulations and Finite Element computations compared with prototype site measurements (see figure 1).

Prototype site measurements
In order to evaluate the accuracy of the numerical approach used for the fluid-structure interaction, measurements on a prototype Francis turbine with a specific speed of n q ≈ 24 min −1 were performed. They were done in different operating points at nearly constant head in steps of 10 MW or 20 MW from the machine start-up to the maximum power of 180 MW and back again. The Francis machine consists of the spiral casing, 10 stay vanes, the wicket gate with 20 guide vanes, the runner with 15 blades and the draft tube. The operating parameters such as head, output power, guide vane opening, the pressure in front of the spiral casing and at the draft tube outlet were recorded during the measurements. The stresses on the runner were obtained by eight strain gauges attached to the suction and pressure side of one runner blade (see figure 2 left). This blade's trailing edge region was measured with a 3D measurement arm and the geometry of the original CAD model has been adapted to increase the accuracy of the FEM simulations. The data of the strain gauge measurements were recorded with a sampling rate of 1 kHz using a data logger, which was mounted in a casing on the rotating shaft inside the draft tube (see figure 2 right).  [13] and Lenarcic et al. [14]. Two different types of physical models are used for the CFD calculations: A full model adjusted according to the prototype machine, including the spiral casing (SC), the stay vanes (SV) and guide vanes (GV), the runner (RN) and the draft tube (DT). In comparison, a reduced model with a cyclic domain for the stay vanes excluding the spiral casing is used. The inlet region is therefore extended in radial direction to avoid impacts from these boundary conditions. A similar approach has been performed by Lenarcic et al. [15]. The particular domains of both models are spatially discretized in multiblock hexahedral grids using Ansys ICEM 15.0. The non-conformal domains are connected at the interfaces (IF) using the arbitrary mesh interface (AMI) implemented in OpenFOAM (see figure 3). The whole mesh consists of about 5.9 million cells for the full model and 5.6 million cells for the reduced model. The quality of the meshes are evaluated regarding the minimum determinant, angle and aspect ratio. The number of cells and the quality of the meshes are summarized in table 1.   In order to estimate the uncertainties due to the discretization of the CFD model, a grid independence study has been performed with the full model according to the approach of Celik et al. [16]. The procedure for the mesh refinement and the evaluation of the results have been already published by Doujak et al. [12]. The results reveal a monotonic convergence for the three global parameters head, runner torque and efficiency with a decreasing grid size h i (see figure 4). Considering the computational costs and the accuracy of the CFD simulations, the medium grid (as described in table 1) with a normalized size of h 2,n = 1.32 is used for the further investigations.
The unsteady CFD simulations are performed in seven operating points with varying output power P , guide vane angles Φ and discharge Q from full load at 180 MW to low load at 20 MW. The meshes for the wicket gate domain are adjusted to the according positions of the site measurements. The operating parameters for the CFD simulations are listed in table 2. Due to unsteady flow conditions at off-design points, a transient rotor-stator analysis with a rotating runner domain is used. To improve the convergence, the unsteady simulations are initialized using the results of steady-state calculations. These are performed with the multiple frame of reference (MRF) solver, which couples the rotating and stationary interfaces with the frozen-rotor approach. The pressure-velocity coupling is done with the SIMPLE algorithm using moderate under-relaxation factors. The unsteady simulations are performed with the PIMPLE algorithm, which couples the SIMPLE and PISO procedure. Therefore, a stable solution can be reached even at increased physical time steps and higher Courant numbers due to the advantage of inner iteration loops together with under-relaxation. The time step is chosen according to a runner rotation of 0.36 • per time step, which yields a maximum Courant number of Co max = 15.

CFD Results
The results of the relative efficiency, compared to the best measured efficiency, at the different operating points for the two physical models are displayed in figure 5 on the left side. The characteristics show a good agreement between the site measurements and the CFD simulations. At full load operation the relative values for the efficiency are slightly lower than the measured ones, even the full model reveals more accurate results. At part load operation, the values are also in a good agreement to the site measurements due to the adjusted CFD setup. The relative  In figure 6 the streamlines in the runner and draft tube domain for the operating points at 180 MW and at 60 MW are illustrated. The displayed velocity is calculated in the relative frame of reference. At full load there are hardly any vortices visible in the runner and draft tube. In part load operation flow separation appears at the inlet of the runner due to the smaller discharge and the inadequate incident flow. The circumferential velocity at the draft tube inlet is also higher than at full load, which decreases the efficiency of the turbine. For the structural analysis the pressure distribution on the runner surfaces is necessary to define the boundary conditions in the different operating points. Therefore, the static pressure on one runner blade and the according hub and shroud region is recorded during the unsteady CFD simulations and averaged afterwards over ten runner rotations, to obtain a steady distribution. The normalized pressure regimes on the blade for two operating points at full load (180 MW) and part load (60 MW) are displayed in figure 7. In order to increase the accuracy of the numerical simulations, second order elements with a quadratic interpolation function are used. The mesh is refined in critical regions like the fillet between the blade an the hub respectively shroud or at the trailing edge of the blade, where the strain gauges have been applied. The sector model is connected at the interfaces using cyclic symmetry conditions. The runner is fixed at the bolt circle, centrifugal and gravitational forces are considered too. The flow in the side chambers of the runner has not been considered in the CFD simulations, but the pressure distribution is calculated analytically and applied to the appropriate FE surfaces, in order to increase the accuracy of the model. On the blade, hub and shroud surfaces, the pressure distribution from the CFD simulations is interpolated between the fluid and structural meshes and applied to the FEM model. The material of the runner is a stainless steel with an E-modulus of E = 216 kN/mm 2 , a Poisson ratio of ν = 0.30 and a density of ρ R = 7.700 kg/m 3 . The setup for the FEM simulations is summarized in table 4.   Figure 9 shows the mean stresses σ m from the site measurements and the FEM simulations for the strain gauge positions PS4 (left) and SS4 (right). The values are normalized to the yield strength σ y of the runner material. According to the CFD analysis, there is hardly any difference visible for the two physical models (full and cyclic). The stress curve at the strain gauge location PS4 reveals a larger gradient along the power curve compared to the measurements. At full load conditions, the stresses are higher, whereas at part and low load operation the values are below the measured ones. For the strain gauge location SS4 the results for the simulations are uniformly higher along the load conditions but the tendencies are in a good agreement. In figure 10 a linear regression analysis is used to determine the accuracy between the measurements and the simulations. The comparison for all operating points and strain gauges (left) shows a good agreement, even if the slope of about 1.27 is higher than 1. At 180 MW power (middle) there is a shift in the vertical direction but the gradient is close to 1. At 60 MW part load (right) there is also an accurate agreement. The equivalent von Mises stresses σ e on the Francis runner in two operating points for 180 MW and 60 MW are displayed in figure 11. The maximum values at 180 MW (left picture) appear with about 20 % of the yield strength σ y at the trailing edge close to the hub of the runner. At 60 MW (right picture) the maximum stresses are lower than at full load and occur in the middle of the blade at the fillet to the hub. However, the yield strength for the runner material is not reached in any operating point and the appearing stresses are all in the elastic regime of the stress-strain curve.

Conclusion
Unsteady CFD simulations and static FEM simulations have been performed for a high head prototype Francis turbine in different operating points from full load to part load conditions according to site measurements. The numerical flow analysis have been carried out using a full model including the spiral casing and a reduced model with a cyclic inlet boundary condition. Further, a grid convergence study has been performed to estimate the numerical uncertainties due to the discretization and to reduce the computational costs. The CFD results reveal a good agreement for the global parameters and small deviations from the measured head, torque and efficiency with the used OpenFOAM setup and both physical models. In order to determine the strains on the Francis runner at different operating conditions, static FEM simulations have been performed with the open-source software Code Aster using a cyclic sector model. The results reveal a varying behavior regarding the magnitude and location of the maximum stresses due to different pressure distributions in the corresponding operating points. The accuracy of the static FEM simulations has been determined by a comparison with strain gauge measurements performed on the related prototype machine, revealing an appropriate agreement. The achieved results of the numerical investigation are further used to assess the lifetime of the Francis runner due to unsteady pressure effects like the rotor-stator interaction or draft tube instabilities, which will be discussed in future publications.