Hydrodynamic damping and stiffness prediction in Francis turbine runners using CFD

Fluid-structure interaction (FSI) has a major impact on the dynamic response of the structural components of hydroelectric turbines. On mid- to high-head Francis runners, the rotor-stator interaction (RSI) phenomenon has to be considered carefully during the design phase to avoid operational issues on the prototype machine. The RSI dynamic response amplitudes of the runner are driven by three main factors: (1) pressure forcing amplitudes, (2) excitation frequencies in relation to natural frequencies and (3) damping. All three of the above factors are significantly influenced by both mechanical and hydraulic parameters. The prediction of the first two factors has been largely documented in the literature. However, the prediction of hydro-dynamic damping has only recently and only partially been treated. Two mode-based approaches (modal work and coupled single degree of freedom) for the prediction of flow-added dynamic parameters using separate finite element analyses (FEA) in still water and unsteady computational fluid dynamic (CFD) analyses are presented. The modal motion is connected to the time resolved CFD calculation by means of dynamic mesh deformation. This approach has partially been presented in a previous paper applied to a simplified hydrofoil. The present work extends the approach to Francis runners under RSI loading. In particular the travelling wave mode shapes of turbine runners are considered. Reasonable agreement with experimental results is obtained in parts of the operating range.


Introduction
Fluid-structure interaction (FSI) has a major impact on the dynamic response of the structural components of hydroelectric turbines. On mid-to high-head Francis runners, the rotor-stator interaction (RSI) phenomenon has to be carefully considered during the design phase to avoid operational issues on the prototype machine [1], as illustrated in Figure 1. The RSI dynamic response amplitudes of the runner are driven by three main factors: (1) pressure forcing amplitudes, (2) excitation frequencies in relation to natural frequencies and (3) damping. All three are significantly influenced by both mechanical and hydraulic parameters. The prediction of the first two factors has been largely documented in the literature, e.g. [2] [3]. However, the prediction of hydro-dynamic damping and added stiffness due to the flow has only recently and only partially been treated in the domain of hydraulic turbines.
Some authors have employed full two-way fluid structure interaction approaches. While these are conceptually straight forward and accessible through commercial software packages, in practice they prove to be difficult to tune [4] when applied to hydro-dynamic structures as well as requiring very high computational resources [5]. Satisfactory results can however be obtained as the two mentioned references show. A large body of literature exists on fluid-structure interaction in general, summarized e.g. in monographs such as [6], treating a variety of topics including galloping, flutter, valve and gate selfexcitation to name but a few. In thermal turbomachines, such as gas turbines or aircraft engines, FSI predictions have been part of the design process for many years. As an example, Moffat [7] describes a simplified approach to calculate the forced response for a transonic fan blade. In addition to forced vibrations, blade flutter has been of particular interest in thermal turbomachinery. Vasanthakumar [8] for example describes a method to compute aerodynamic damping on the blades of a transonic low bypass fan to assess flutter risk using a modal work based approach. The methods employed in the field of thermal turbomachinery are also applicable to hydraulic turbines, but some adjustments may be required.
Regarding the flow-added dynamic parameters on hydro-foils there have been a number of investigations in recent years. Monette [9] derives the theoretical equations governing added dynamic parameters and applies them to a plate-like hydrofoil with good comparison to experimental results. Münch [10] presents a CFD approach to determine the flow-added dynamic parameters for a rotating 2D hydrofoil. A frequency sweep and complex transfer functions are used for that purpose. More recently CFD modal work based studies to determine damping on turbine blades have also been performed [11]. A single degree of freedom rotating oscillator model coupled to a 2D CFD computation is used in [12] to assess pump turbine guide vane stability regarding self-excitation.
A number of experimental studies, some quite dated [13] and some relatively recent [14]- [16] are also available. These mostly show investigations of geometrically simple structures such as plate-like or 2D hydrofoils and represent good cases for the validation of prediction methods.
Other than the highly resource-intensive full two-way FSI calculation of a Francis runner, no studies are known to the authors that deal in detail with the specifics of predicting flow-added dynamic parameters of Francis runners.

General
Due to the above mentioned high resource requirements of full 2-way coupled FSI approaches [5], in particular in terms of computation times, two simplified approaches of capturing the interaction between the flowing water and the vibrating structure of Francis runner blades have been investigated. They are both based on the reduction of the dynamic behavior to individual vibration modes. These are determined by means of modal analysis using FEA of the runner in still water. In Francis runners the modes are associated with nodal diameters. In the case of RSI excitation the nodal diameter ND that can be excited is known. It is a function of the number of guide vanes N gv and runner blades N bl (ND=±m•N bl -n•N gv , m,n=0,1,2…). The normalized travelling wave mode shape vector ̃ imposed by the rotating RSI pressure pattern and resulting from the linear combination of two orthogonal standing wave mode shapes with the same natural frequency, can be expressed as a function of the location vector x of the individual blades l: By normalization we understand ̃. to be dimensionless with its maximum local norm max (‖̃. ‖) equal to 1. The modal deformation of a considered blade is part of a nodal diameter runner deformation and thus, the motion of a single blade is not perfectly in phase as this is often assumed for motion of axial thermal turbomachinery blades. However, the large deformations and even more the reaction of the fluid will have a clearly dominant phase on the blade. The complex nature of the vector for the modal deformation per blade is apparent by the fact that points on the blade surface do not reach their local maximum value at the same time. The real travelling wave modal motion of the blades in a CFD model can then be imposed as real displacement Here ̂ scales the travelling wave modal motion and is of dimension length. In order to extract the equivalent flow-added dynamic modal parameters, the reaction of the unsteady fluid flow to the modal motion can be approximated -considering the dominance of the real part of the mode in typical Francis blade modes -as a single degree of freedom (1dof) oscillator by Here is the time dependent reference displacement which is chosen to be the instantaneous local displacement on the surface of interest where max(‖̃. ‖) = 1. The modal mass of the structure m s and of still water m w is obtained by FEA. The damping coefficient contains contributions from structural damping B s , damping that occurs in still water B w resulting from acoustic radiation and viscous effects, and damping resulting from momentum exchange between the flow and the structure as a result of the oscillation B f . The first two, B s and B w , can typically be neglected as they are 2-3 orders of magnitude smaller than B f under flow velocities in hydraulic turbines. Modal stiffness includes contributions from the structure k s as well as the water k w (compressibility effects) both obtained from FEA, and the flow k f . Figure 2 illustrates the concept and terminology for a 1dof motion of a plate. a b Figure 2. Concept and terminology for modal representation of hydro-foils (example damping plate H0 [14]): a) Example of mode shape (spatial distribution) b) Modal representation.
The modal work is the energy that is exchanged between the vibrating mode and the flow during any given period of time (usually used as work per cycle, e.g. [8]): Here A is the surface area of the blade. In a CFD calculation ̇ can directly be interpreted as the mesh velocity. The shear only contributes at free edges where it is aligned with the motion. On Francis turbine blades its contribution is very small (2-3 orders of magnitude less than pressure contribution) and can generally be neglected. In the following sections two CFD based methods are described that permit the extraction or consideration of the effect of flow-added parameters onto the dynamic behavior of blade modes under RSI loading. The fundamental assumption underlying these methods is that the shape of the mode as determined in still water is not changed by the flow.

Modal work approach
The travelling wave modal motion of the blades according to equation (2), reacting to RSI pressure pattern, is imposed during an unsteady CFD calculation. Since the CFD is incompressible and does not contain the structure, the part of equation (3) that is represented in the CFD solution reduces to This gives the following equation for modal work in modal space, the value of which needs to be equal to that of equation (4) obtained from CFD: This means that the exchanged energy is represented by equivalent 1dof parameters. Using u r =û r e it as ansatz for the harmonic motion we get the following expression for the modal work exchanged in an arbitrary time period:  s is the imposed angular oscillation frequency. One can see that the expression inside the integral contains a constant term and two terms with the sin and the cos of 2 s . This implies that in the time derivative of the modal work dW m /dt there will be an oscillation of twice the frequency of the displacement u r (t). Transforming dW m /dt into the frequency domain by means of a Fourier transform permits the extraction of values for the determination of the flow-added parameters: The flow-added parameters can also be obtained by integrating equation (7). Integration over the oscillation period T will yield the damping and integration over quarter periods synchronized with the displacement (0…T/4, T/4…T/2 etc.) will yield expressions to determine the added stiffness. The damping coefficient can be converted into the more convenient non-dimensional damping ratio (ratio of damping to critical damping ). It is important to use the correct modal mass and natural frequency reference. Using those in flowing water should be the most appropriate: Here  n.w and  n.wf are the angular natural frequencies in still and flowing water respectively. The sign of W m.T (=2Re) is negative when the structure transfers energy to the flow and would therefore see declining amplitudes in a free vibration. To make  positive for this situation the minus sign is introduced in accordance with [7]. For the frequency in flowing water follows:  Second part: interface between rotor and stator has been switched to "Stage"  free oscillation.
An advantage of the modal work approach is that all flow-added parameters can be obtained from a single CFD computation. An alternative approach to compute the added stiffness is the use of the static force over displacement characteristic by k f =F m /u r [11]. This however requires at least two CFD computations and does not consider any secondary dynamic effects.

Coupled 1dof approach
In the coupled 1dof approach the modal motion is not imposed in the CFD calculation. Instead the mechanical dynamic parameters are coupled to the flow by means of a 1dof oscillator equation which is integrated numerically together with the unsteady CFD solution. Again the equation to be coupled to the flow is derived from equation (3): The modal dynamic parameters included in the equation stems from the mode in still water. Modal mass of water m w , damping B f and added stiffness k f are included in the solution and enter the equation by means of the modal fluid force F m : The modal force for the coupled 1dof approach is calculated with the real part of the mode shape only. This is an approximation, which in practice works reasonably well if the mode shape is turned such that . ( ) is at its maximum displacement at the reference location. Equation (11) can then be integrated numerically with the following approach   2 1 . .
A real valued approximation of the displacement imposed at every time step then follows from a single standing wave mode shape by: (14) One can see that here the characteristic of points on the mode shape reaching their maximum at different points in time, i.e. the travelling wave nature expressed by a variable phase, is lost. This simplification and approximation is justified if in one section of the runner the phase difference compared to the maximum displacement point is small for all locations with significant displacement. In Francis turbine runners this is generally the case as the example of a typical mode shape in Figure 4 shows. Using two 1dof equations does not really solve the problem of losing the wave nature of the modal displacement because both modes are only coupled by the fluid and are therefore synchronized with the excitation. To correctly impose the wave-like modal motion one would have to know the response frequency, but this is part of the solution. Determining instantaneous values numerically is relatively difficult. Figure 3b shows an example of the two types of results that can be calculated by means of the 1dof coupling to CFD: When coupled to an RSI calculation the forced response is obtained. If the rotor-stator interface is switched to "mixing plane" (circumferential averaging), a free oscillation results, from which damping and response frequency can be extracted. The advantage of this approach over the modal work approach lies in the fact that results can be obtained on the actual unsteady flow rather than on the mean flow. a) b) Figure 4. Example of typical blade mode in a) with amplitude (left) and phase (right). Same mode represented in b) as real (left) and imaginary (right) parts. Mode is turned such that Re is at its maximum.

Setup
The commercial CFD software Ansys CFX v16 is used for this study. Both described approaches to account for flow-added dynamic parameters are based on unsteady CFD computations. Figure 5 shows the typical setup for a modal work computation on the left. Three runner channels and one guide vane channel are used. Guide vane and runner domains are coupled by means of a stage (mixing plane) interface. Therefore the computation represents the mean flow with the unsteadiness added by modal motion. Stability based unsteady phenomena such as von Karman vortices can also form but RSI is eliminated due to the stage interface. For the coupled 1dof approach a typical setup as needed for an RSI computation is used as a base such as the one shown on the right in Figure 5. The number of guide vane and runner channels will depend on their numbers in the real machine.

Numerical and physical parameters of influence
The influence of numerical parameters onto the results has been studied using the damping plate [9] and are shown in Figure 6. Standard size meshes in the order of 200k nodes for a runner channel appear to be sufficient. The time step size has a significant influence. About 200 time steps per period seem to be a reasonable choice. The vibration amplitude that is imposed during a modal work computation also has a significant effect. Results from the damping plate show that a minimum is reached near 0.5mm amplitude. Computation times obviously depend on the number of processors used but are generally high. A modal work calculation will take about 12 hours while a coupled 1dof calculation requires twice that.  Figure 5. CFD setup for a modal work computation (left) and coupled 1dof computation (right). It is sufficient to use three channels with the correct inter-blade phase angle for the modal work approach as shown in the middle: using 5 runner channels gives same result for all three middle channels.

Damping plates
The first type of test case is based on the damping plates introduced in detail in [9] and [14]. Good quality measurements exist for these damping plates. They exhibit straight forward real mode shapes, well suited for the validation of the basic concepts. Comparisons of experimental and CFD results for the plates H0 and H3 with two mode shapes M1 (1 st bending) and M2 (1 st torsion) are shown in Figure  7. Very good results are achieved for H0M1 both for damping and frequency prediction. For H3M1 and H3M2 the comparison is less good. Damping is over-predicted in the order of 20%. However, the tendency as a function of flow velocity and between modes is well captured. The measured data exhibit a certain degree of scatter and the predicted values come very close to the upper limit of these data points. Using the 1dof coupling gives similar results to the modal work approach. The modal work approach appears to give slightly better comparison for the frequency shift.

Runner blades
Experimental data for the damping values on runner blades are inherently less precise since they can only be indirectly calculated based on measured strains, thus including uncertainties from the natural frequency in water, the geometry, the damping and the flow-added stiffness. Comparisons here are made for operating points where the vibration response is near its peak and often close to resonance. That eliminates uncertainties regarding the natural frequency and allows more accurate assessment of the damping alone. The ratios of CFD-predicted to measured values for damping and maximum displacements are shown in Figure 8 for a specific Francis runner. At best efficiency both the modal work and the coupled 1dof approach achieve good results. At part load the modal work approach gives a significant over-prediction of the damping of about 30%. The coupled 1dof approach only slightly over-predicts the damping and to a similar extent under-predicts the maximum amplitude under forced response. This may be due to the fact that the underlying CFD solution is RSI rather than the mean flow. Even in free oscillation, the initial solution is still the one from RSI. The ratio of the CFD-predicted natural frequency in flowing water to the FEA-predicted natural frequency in still water is also shown. It is only slightly above 1 and increases slightly with power. No accurate measurement data are available for the natural frequency in flowing water at this point. a) b) c) Figure 7. Damping plates. a) H0M1 damping ratio from CFD modal work (blue) and coupled 1dof approach (green) compare very well with experiment. b) The increase in natural frequency from modal work compares slightly better than from 1dof. c) H3 results for the damping of modes M1 (1 st bending) and M2 (1 st torsion) show less good comparison than for H0 (overprediction ~20%) measured values. Figure 8. Ratio of the predicted over measured damping  CFD / Exp , of the predicted natural frequency in flowing water over the FEA predicted natural frequency in still water f nwf /f nw , and predicted over measured reference amplitude u r.CFD /u r.Exp . Overall good comparison is achieved near the best efficiency point. Towards part load, the deviations increase.
Using predicted damping obtained by the modal work approach, RSI forced response of three different Francis runners was performed using standard procedure [2]. Figure 9 shows a comparison of the predicted and measured strains that where obtained for strain gauge measurements on those Francis runners. Generally a very satisfactory comparison is achieved. In particular the critical maximum strains are well predicted. Figure 9. Correlation of predicted and measured strains at various strain gauge locations for measurements on three runners using the damping values calculated with the modal work approach. Satisfactory agreement is achieved.

Discussion and conclusion
Two approaches for the prediction of flow-added parameters based on CFD are presented. Both rely on the concept that the primary dynamic aspects can be represented by individual modes which have to be determined by means of FEA. The shape of these modes is assumed to remain unchanged by the flow. This of course is not fully correct. However, the damping results both for the simple plate and the runner blades show a reasonable comparison with experimental data. This leads us to the conclusion that both approaches are generally appropriate for the prediction of hydro-dynamic damping in Francis runners. Regarding flow-added stiffness, appropriate measurement data will be required in the future for full validation of this specific parameter.