Challenges in Dynamic Pressure and Stress Predictions at No-Load Operation in Hydraulic Turbines

Some of the potentially most damaging continuous operating conditions for hydraulic turbines are the no-load (NL) conditions. At NL conditions the flow passes through the turbine without power generation, but with non-negligible flow rate, and therefore all the potential energy in the flow has to be dissipated. This takes place through a mechanism where the runner channels are partially pumping, thus generating large scale unsteady vortex structures which, by their nature, break down into smaller and smaller vortices until energy dissipation occurs at the smallest scales. This type of flow, dominated by its turbulent character, is inherently difficult to simulate by means of numerical methods since turbulence model and numerical dissipation have a major influence. The resulting dynamic loads on the runner are largely of stochastic nature, exciting a broad band of frequencies and thus, almost always interact with at least one deformation mode. The presented investigations are aimed at predicting the effect of the unsteady NL pressure loads on the fatigue life of a Francis turbine runner. A combination of computational fluid dynamics (CFD) and finite element analysis (FEA) methods has been employed. The results from transient CFD simulations are presented. Comparison of the results with prototype strain gauge measurements at no load conditions shows that the stochastic nature and the approximate range of the dynamic stresses can be predicted.


Introduction
No-Load (NL) operating conditions occur under various circumstances in the operation of Francis turbines. Most notable are speed-no-load, runaway and transients such as load rejections. At speed-noload the turbine is operated at synchronous speed without electricity production typically in stand-by ready to be connected to the grid. This operation can occur over extended periods of time. Runaway is typically an exceptional condition which occurs if a generator is rejected from the electrical grid and the guide vane mechanism fails to close. However, runaway can also be used on purpose for example in cases where water levels need to be maintained after disconnection of a machine from the electrical grid often referred to as sluicing or ramping. During transient events such as load rejections the machine can temporarily pass through a NL condition.
As reported in [1] even the mean contribution of speed-no-load to fatigue life of Francis runners can be significant. However, NL operation is typically highly dynamic. Pressures and strains strongly fluctuate as can be seen in figure 1. Clearly the unsteadiness of this condition with its peak-to-peak pressure and strain variations is relevant for fatigue life as has been reported in the literature [2]- [6]. Some recent studies on conditions around NL have been performed in the field of reversible pump turbines, e.g. [7] [8]. Primary concern in those investigations is the prediction of the characteristic instability in the S-shaped region of the characteristic curve as well as the radial force imbalance on the machine due to rotating stall. Computational Fluid Dynamics (CFD) based prediction methods of some unsteady -mostly periodic -phenomena in hydraulic turbines such as rotor-stator interaction (RSI) and von Karman vortex shedding is well-established standard practice [9] [10]. However, to the authors' knowledge, no studies on the prediction of the stochastic loads that act on the blades or other parts of radial turbines during part load or no load operation have so far been published. b In order to assess the ability of a machine to operate under NL condition before fatigue damage will lead to cracking, a reasonable prediction method for the pressure loads is required. Due to the stochastic nature of the loads, an appropriate and efficient conversion of the dynamic pressure loads into deformations and further stresses is equally required. A CFD based method for the prediction of the unsteady stochastic pressure loads is presented here and results are compared to unsteady pressure measurements on a prototype Francis runner. A Finite Element Method (FEM) based approach is used for the conversion of the stochastic pressure loads into stochastic dynamic stresses, the results of which are also presented and compared to strain gauge measurements on the same prototype.

Setup and meshing
The work presented here is the summary of the collaboration and experience of two teams, one from Andritz Hydro and one from Hydro Québec. A number of simulations of NL conditions are presented for which also prototype measurements are available. In order to test a larger number of parameters of influence and therefore to converge faster onto a suitable method, two different general setups are used. They are shown in figure 2.
The first setup consists of the distributor, the runner and the draft tube resulting in a total mesh size of about 5M nodes. All meshes are generated with automatic in-house mesh generators. The runner and draft tube meshes consist of hexahedral elements, while the distributor mesh is a hybrid mesh made up out of hexahedral and prismatic elements. As CFD solver the commercial code Ansys CFX is used. The combination of total pressure at the inlet and static pressure at the outlet is found to be suitable because the known head can be imposed directly. Additional hydraulic losses in the casing are considered to be relatively small at NL conditions and are therefore neglected in this setup. A turbulence intensity of 5% is specified at the inlet. The runner mesh is connected to the meshes of the adjacent components by means of transient rotor stator interfaces. In this study the guide vane opening as well as the runner rotational speed as measured on the prototype is used. If not known, both can be found through an iterative process such as the one described in [1]. The time-dependent flow is analyzed with a second order backward Euler time marching scheme. While this is not the case in reality, the flow is considered single phase, i.e. neither cavitation nor possible aeration are considered directly in the flow solution. This is necessary in a first step in the interest of limiting computation times. As advection scheme, Ansys CFX's "high resolution" scheme is used. All unsteady processes are potentially susceptible to the initial conditions. In this first setup a steady state solution with "Frozen Rotor" interfaces is chosen. A more detailed study of the effect of the initial condition is planned for the future. While the second setup is similar, some significant differences exist. First, the complete spiral casing is included to provide a less uniform flow at the distributor entrance. Second, the runner hub cap is meshed according to the actual geometry rather than being simplified. Fully hexahedral meshes are used in the distributor, runner (Numeca Autogrid) and draft tube (Ansys ICEM-CFD) while a hybrid hexahedral-dominant mesh (Numeca Hexpress Hybrid) is used in the spiral casing (which includes the stay vanes). This last mesh accounts for about 14M elements, while the total number of nodes is just over 20.5M. Also, in this second setup, water is considered as slightly compressible through an equation of state relating pressure and density. For the boundary conditions, the same prescribed head was used, but the inlet is located at the entrance of the spiral casing, while the outlet is an opening with entrainment and a static pressure. A value of 1% turbulence intensity is set at the inlet. The advection scheme used is also Ansys CFXs' "high resolution", but with a minimal value of the blending factor equal to 0.5, to prevent the solver from falling back onto a less accurate first order approach. The so-called "high resolution" scheme is also used for the turbulent equations.

Significant numerical and physical parameters of influence
Ansys CFX is a general purpose solver for the unsteady Reynolds averaged Navier-Stokes (URANS) equations. Consequently most of the turbulent scales are accounted for by means of a turbulence model. The choice of turbulence model can be expected to have a relatively large influence on a highly dynamic flow condition such as the NL condition. Figure 3a) shows the pressures at a measurement equivalent monitor point P2 (see figure 5 for location) calculated with two different turbulence models. Clearly the SAS turbulence model brings out more of the stochastic elements that are characteristic for the investigated NL condition. In particular the maximum pressure range is larger with the SAS model. Where the mesh is fine enough, the SAS model enables the resolution -rather than modelling -of much smaller turbulent scales than simpler turbulence models such as the k-ε model. Figure 3b) illustrates the general degree of disorganization of the flow. Highly variable channel vortices are present, a number of which in this example extend into the upper part of the draft tube resembling multiple ropes. Clearly, many more of the vortices which are responsible for the dynamics of the flow are brought out by the simulation with the SAS than with the k-ε turbulence model. Details on the differences in modeling turbulence with the two models can be found in Ansys CFXs' documentation. As figure 3 shows, for flow under NL when the irregular unsteadiness is of interest, the SAS turbulence model can be considered the better choice. Other turbulence modelling options have also been investigated. The DES model of Ansys CFX did not show any improvement over the SAS model. Reynolds stress models should be a good choice for the type of flow at hand. However, they are computationally significantly more expensive and proved to be unstable since numerical divergence rapidly became a concern.
Another factor that most certainly has an effect on the runner pressure loads is cavitation. Under NL large volumes of vapor are present in the runner and the upper part of the draft tube. Since our analysis is single phase, a first step in assessing the influence of cavitation can be taken by clipping the pressure at vapor pressure for post-processing. Figure 4a) shows in blue the areas on the blades that at one particular time step would have had pressure below vapor pressure. Clearly, this way of accounting for cavitation neglects significant parts of the actual physics, such as the displacement of liquid flow by the vapor and the resulting changes in the complete flow field. Nevertheless, it gives us a first order estimate of the influence cavitation might have. In figure 4b) the effect of using clipped versus unclipped pressure to determine the instantaneous blade torque is shown. It can be seen that the blade torque range based on the clipped pressure is smaller than that of the unclipped pressure.
Although the former should be more physical, we prefer to use the unclipped pressure in order to have pressure load estimates tending towards higher values for conservative estimates regarding dynamic stresses and resulting fatigue life. Preliminary investigations using cavitation modelling within the CFD analysis did not show significant increases of the pressure variations. However more in-depth investigations are needed to for conclusive evidence on this. These are planned and may be reported in a future publication. Given that a large part of the turbulence is modelled rather than resolved, it is interesting to determine the amount of turbulent pressure fluctuations that are "hidden" inside the turbulence model. Preliminary investigations indicated that the pressure fluctuation amplitudes hidden in the turbulence model are relatively small compared to the large scale pressure pulsation amplitudes that are resolved by the URANS code. More detailed investigations on this aspect of CFD modelling of the NL conditions are needed and may be reported in the future. Cleary, even if the pressure pulsations "hidden" in the modelled turbulence quantity (turbulent kinetic energy) are relatively small, this does not mean that the way of modelling turbulence is not important for NL conditions. As figure 3 above clearly indicates and as the results presented below confirm, the turbulence modelling has a large effect on how many of the large vortex scales that drive the large scale pressure fluctuations are resolved and therefore appear in the results of the simulation.
An additional factor that should be considered in the validation is aeration. The investigated machine has an aeration system installed that was most likely operational while the measurements were taken. However, the air admission flow rate was not recorded with the unsteady pressure and strain gauge measurements. In the single phase CFD analysis air admission is necessarily neglected, but it may have to be considered in future measurement campaigns and validations.

Pressure measurements
It cannot be expected that standard URANS based CFD will be able to predict all aspects of the flow at NL condition in full accuracy. The main question is rather to what extent CFD can predict the largest amplitudes that will occur on a real machine and, in a general manner, the stochastic nature of the flow. Of actual interest are the stress amplitudes which can be measured with reasonable accuracy on a Francis runner prototype [6]. For these stress amplitudes to occur, pressure fluctuations have to be present. Measurements of pressure fluctuations can directly be used for validation of the CFD predictions, while strain gauge measurements serve as an indirect means of validating analysis results. Nevertheless, strain gauge measurements are a powerful method of validation because their results represent area integral measures of the effect of pressure fields onto the structure.
For the present study, unsteady pressure measurements from commissioning tests on a prototype Francis runner at NL condition are available. ENTRAN EPN-B13-35B-/X/Z1 pressure transducers are used. The sensors are calibrated to 15V and 35 bars. The acquisition system of the sensors on the rotating components is an EDAQ Plus by Somat. Pressures are sampled at a frequency equivalent to 0.65° runner rotation per time step (°/time step). Further details of the measurement approach are confidential Andritz Hydro information and beyond the scope of this paper. Measurements under NL are challenging for the measurement equipment. Out of 12 measured pressure sensors, after standard initial signal quality checks, 4 were considered to have recorded valid data. The pressure sensor locations are shown in figure 5 while the recorded unprocessed time signals are shown in figure 6a). The time signals look reasonable for all four sensors, containing a large amount of stochastic with a differing range for the four sensors. The Fourier Transforms (FT) for those time signals in figure 6b) show broad spectra. A distinct peak at 20 times runner rotation frequency corresponding to the guide vane passing frequency (RSI) is only present in one of the sensor signals. RSI can be viewed as an indicator signal component. It should be present as it is in the CFD calculations; therefore indicating that probably only sensor P4 delivered valid data in this measurement series. We use the results from sensor P4 for our direct pressure amplitude validations.

CFD results
The time dependent CFD simulations have pressure monitoring points at the locations of the pressure sensors on the prototype. Figure 7a) shows the time signals from the P4 sensor location for the measured values (blue) in comparison with the calculated results. Four calculated results, considering two different setups with two different numerical settings each are shown. One can see that the two setups using the SAS turbulence model can bring out the stochastic nature in the flow. The general pressure range is similar to the measured one in all three calculations using the SAS model (green, red, magenta). In contrast, the calculation using the k-ε turbulence model (light blue) shows a much Computation times are long (in the order of 3 to 7 days); therefore only about 10 runner rotations were calculated with 1° runner rotation per time step (°/time step). The red line figure 7a) represents the continuation of the magenta line but with a change in time step size corresponding to 5°/time step. In this way about 45 runner rotations can be computed in less computation time than needed for the 1°/time step (magenta line). Fourier transforms of the pressure monitor results in figure 7b) show that a strong RSI signal is present in all of them. Taking a closer look at the RSI signal on the right of that figure, it becomes clear that with 5°/time step (red) some of the RSI pressure amplitude is lost. All other settings, even the one with the k-ε model -as expected from previous experience [10] -give RSI amplitudes similar to the measured one. The k-ε results show much less spectral density in the low frequency range than the measurements, which is the consequence of the practically nonstochastic appearance of the time signal. The other calculation settings bring out a reasonable amount of low frequency amplitudes in comparison with the measurements when inspecting the Fourier transforms in figure 7.
A point of interrogation is the best time step size to use for NL computations. If RSI is to be reproduced correctly, the best accuracy is achieved with 1° runner rotation per time step. (Smaller time steps may yield even better results for RSI, e.g. [10], but RSI is not the focus of this study.) However, it appears that if the main interest is the maximum range under NL, a larger time step corresponding to 5° runner rotation per time step is equally suitable and more time economical. In the example of the red line in figure 7a) the maximum range in this particular sequence actually exceeds the measured range. This indicates the importance to simulate long time sequences in order to have a higher probability of including in the calculated time a superposition of the various vortices and flow structures that lead to the high pressure peak events which will dominate the fatigue assessment.
It is important to note that Fourier transforms are not the most suitable way of analyzing primarily random or stochastic signals (Instead, FT's are well suited for signals with periodic content). One suitable method for the analysis of stochastic processes, especially in the context of fatigue analysis, is that of rainflow counting, e.g. [2]. Figure 7c) shows a rain flow diagram comparing the measured (blue) with two calculated pressure signals at P4 (red: setup one, SAS, 5°/time step; green: setup two, SAS, 2°/time step). In order to compare rainflows from measurements and from different CFD calculations, time signals are extrapolated to a unique length using the method presented by Poirier [11]. Both CFD calculated curves show a relatively similar tendency over most of the number of cycles but the green curve from the 2°/time step analysis overlaps more closely with the curve from the measurements. At the largest pressure ranges, which occur with the lowest number of cycles, the prediction -in this example -exceeds the measurements. These large range cycles are of highest significance for fatigue life assessment. Therefore it is desirable to have an over-rather than an underprediction in order to obtain conservative estimates for the fatigue life. In summary, when the stochastic is present in the CFD calculation, such as in the red and green curves of figure 7, the quality of the results is quite satisfactory in terms of range of pressure fluctuations. In the interest of a fatigue evaluation based on rainflow cycle counting, the large peak events are of greater interest than a good resolution of RSI. Therefore the calculation of a large number of runner rotations is important, and in this example the use of 5°/time step is a reasonable choice. With even larger time steps too much information is lost. However, if the computational resources are available, a 2°/time step is a better overall choice.
An additional way of presenting the calculation results is in form of the dynamic runner blade torque as shown in figure 8. An advantage of the blade torque is the fact that it represents an integral value, i.e. is the result of the integration of the instantaneous pressure over the blade surface (weighted with the local radius). It is therefore more closely related to what will load the blades mechanically rather than the pressure measured at a single location. The results shown in figure 8 confirm the findings from the pressure monitors: The stochastic nature of the flow resulting in large irregular blade torque amplitudes is well predicted with both setups using the SAS turbulence model. In contrast, the k-ε model predicts much lower blade torque variation in the time signal while in the Fourier transform practically only RSI and one lower frequency are present. The quality of the prediction of the blade torque amplitude at guide vane passing frequency is primarily a function of the time step size. Both investigated turbulence models are suitable for the prediction of the amplitudes at this frequency.

Structural modelling and stress predictions
Once CFD pressure load time signals are available for a sufficient number of runner rotations, they can be converted into deformations and stresses. Flow structures and pressure fluctuations due to the NL flow are relatively large (see vortex core region visualization in figure 3). Consequently, it is assumed that the feedback of the deforming mechanical structure (e.g. blades: deformation in the order of 1mm over characteristic blade span length of 1 to 3m) onto the flow can be neglected. This leaves one-way fluid structure coupling as the method of choice. At the current point in time, this is also necessary in the interest of computational resources and numerical stability for a flow condition as complex as the NL condition. The pressure fields are transformed into loads for mechanical analysis and converted into stresses as a function of time by means of a numerical model. (A detailed description of the method of converting CFD pressures into blade stresses would exceed the available space for this paper. It will be reported in the future.) Once the stresses are available as a function of time, rainflow counting can be performed for comparison with strain gauge measurements. a b Figure 9. a) and b) No-load condition rainflow comparisons of strain gauge positions 1 and 3, measured blue on two different blades (one curve per blade) and predicted red (setup1 SAS 5°/time step). c) Statistical evaluation of all available strain gauge measurements. Measured over calculated values with slope 1.09 of regression line (black), 95% confidence interval of regression curve (red) and 95% confidence interval of the data (green).
As mentioned above, using rainflow diagrams is a convenient way of representing the dynamics of stochastic processes in particular with respect to fatigue assessment. Figure 9a) and b) show the results of dynamic stress ranges as a function of number of cycles for two strain gauge locations. Lines in blue represent the measured values and lines in red the predicted values. The general tendency is predicted by the analysis method. For fatigue considerations, the largest stress ranges, which typically occur with the lowest number of cycles, are the most critical. One can see that the maximum value of the range is reasonably well predicted. On the graph of figure 9c), the slope of the linear regression of the measured over the calculated strains is 1.09. The correlation between the measured and simulated data is similar to what is obtained for standard static stress analysis. This is an encouraging accuracy regarding future systematic predictions of dynamic stress under no-load operating conditions. However due to the small number of data points and the relative scatter, the confidence interval is rather big. Therefore, in order to conclude that the method is reliable, a larger number of validation cases are required.

Discussion and conclusion
Inevitably, no-load operating conditions will occur in the life time of a hydraulic turbine, either at machine start, during load rejection or in exceptional cases as a result of ramping or sluicing operation. As has been seen from previous experimental studies (e.g. [2]) the stochastic pressure loads that occur under such no-load operating conditions are significant under structural fatigue considerations. Therefore it is important to be able to predict these stochastic pressure loads in order to assess turbine fatigue life during the design phase where there is still the possibility to take corrective measures, rather than based on strain gauge measurements on an existing prototype. As is apparent from the limited number of prior studies, the prediction of stochastic loads such as the ones that occur under no-load operation, poses difficulties using standard computational fluid dynamics methods. In order to obtain a relatively good prediction of the apparently random large scale pressure fluctuations, long enough real times need to be calculated and computational meshes need to be relatively large. This results in high resource requirements. However, our study shows that the prediction of significant aspects, such as the maximum pressure and stress ranges, as well as the general stochastic nature of the flow, can be predicted to a reasonable accuracy using a suitable setup within a commercial CFD code. These preliminary results are encouraging, although a significant number of additional validation cases will be needed to assess the reliability of the approach before a standard prediction of this kind of load can be included in the design process. Therefore it is planned to perform additional validation with already existing measurements. Additional physics such as cavitation and aeration may also be included.