Numerical Estimation of Torsional Dynamic Coefficients of a Hydraulic Turbine

Abstract

The rotordynamic behavior of a hydraulic turbine is influenced by fluid-rotor interactions at the turbine runner. In this paper computational fluid dynamics (CFDs) are used to numerically predict the torsional dynamic coefficients due to added polar inertia, damping, and stiffness of a Kaplan turbine runner. The simulations are carried out for three operating conditions, one at about 35% load, one at about 60% load (near best efficiency), and one at about 70% load. The runner rotational speed is perturbed with a sinusoidal function with different frequencies in order to estimate the coefficients of added polar inertia and damping. It is shown that the added coefficients are dependent of the load and the oscillation frequency of the runner. This affect the system's eigenfrequencies and damping. The eigenfrequency is reduced with up to 65% compared to the eigenfrequency of the mechanical system without the fluid interaction. The contribution to the damping ratio varies between 30–80% depending on the load. Hence, it is important to consider these added coefficients while carrying out dynamic analysis of the mechanical system.

1. Introduction

Thomas [1] initiated the research on fluid-rotor interactions on turbines in 1958. He suggested an analytical model of destabilising forces due to nonsymmetric clearance in steam turbines. Alford [2] developed a similar model for compressors, where the forces are obtained as a function of the change in efficiency due to increased eccentricity. Ulrich [3] carried out the first research in a test rig and suggested corrections to Thomas and Alford's models. At the same time Iversen et al. [4], Agostinelli et al. [5], and Csanady [6] introduced models of hydraulic unbalance forces due to asymmetry of the flow channel geometry in centrifugal pumps. Hergt and Krieger [7] studied the influence of radial forces during off-design operating conditions. Colding-Jorgensen [8] used potential flow theory to determine damping and stiffness coefficients. Adkins [9] were the first to introduce an analytical model of both mass, damping and stiffness coefficients and harmonic forces. Adkins and Brennen [10], and Bolleter [11, 12] used test rigs to continue the development of models for fluid-rotor interactions of pump impellers. Childs [13] used bulk flow theory to determine rotordynamical coefficients at the pump-impeller-shroud surface.

The use of computational fluid dynamics (CFD) has recently increased within the area of fluid-rotor interactions. It was introduced by Dietzen and Nordmann [14] in 1987, but has due to the computational cost not been widely used in the past. The first applications of CFD within rotordynamics have been in the area of hydrodynamic bearings and seals. Recently, CFD has entered into the research of fluid-rotor interactions in centrifugal pumps [15]. CFD has been more common in research and development of hydraulic machinery. Ruprecht et al. [16, 17] used CFD to calculate forces and pressure pulsations on axial and Francis turbines. However, the results were not used in rotordynamical analysis. Liang et al. [18] carried out finite-element fluid-structure interactions of a turbine runner in still water and showed a reduction of the nonrotating eigenfrequencies compared to a runner in vacuum. The result had good agreement with the experimental results presented by Rodriguez et al. [19]. Karlsson et al. [20] analyzed the influence of different inlet boundary conditions on the resulting rotordynamic forces and moments for a hydraulic turbine runner. The benefits of using CFD to calculate rotordynamical forces and coefficients of hydraulic turbines have not yet been fully explored. In the present work CFD is used for the determination of the torsional dynamic coefficients due to the flow through the turbine.

2. Modelling and Simulation

2.1. Fluid-Dynamical Model

2.1.1. The OpenFOAM CFD Tool

In the present work the OpenFOAM (www.openfoam.org) open source CFD tool is used for the simulations of the fluid flow through the Hölleforsen water turbine runner. The simpleFoam OpenFOAM application is used as a base, which is a steady-state solver for incompressible and turbulent flow. It is a finite volume solver using the SIMPLE algorithm for pressure-velocity coupling. It has been validated for the flow in the Hölleforsen turbine by Nilsson [21]. New versions of the simpleFoam application have been developed in the present work, including Coriolis and centrifugal terms and unsteady RANS. All the computations use wall-function grids and turbulence is modelled using the standard 𝑘−𝜖 turbulence model. The computations have been run in parallel on 12 CPUs on a Linux cluster, using the automatic decomposition methods in OpenFOAM. The version number used for the present computations is OpenFOAM 1.4.

2.1.2. Operating Conditions

All the computations are made for the Hölleforsen Kaplan turbine model runner, shown in Figure 1. The computational grid is obtained from earlier calculations by Nilsson [21]. The operating conditions used for the present investigations are for runner rotational speeds of 52rad/s, 62rad/s, and 72rad/s, which correspond to loads of about 70%, 60%, and 35%, respectively. The boundary conditions are kept the same for all operating conditions (in the inertial frame of reference). The change in the load due to the rotational speed is explained by the fact that the pressure drop (or head of the system) needed to drive the same flow through the turbine will change with different rotational speed. The runner rotational speed is finally perturbed with a sinusoidal function in order to identify added coefficients for the torsional dynamic system. This is described below.

Figure 1 

The computational domain.

2.1.3. Boundary Conditions and Computational Grid

The inlet boundary condition was obtained by taking the circumferential average of a separate guide vane calculation, yielding an axisymmetric inlet flow [22]. This corresponds to a perfect distribution from the spiral casing and without any disturbance from the guide vane wakes.

Wall-functions and rotating wall velocities were used at the walls, and at the outlet the homogeneous Neumann boundary condition was used for all quantities. Recirculating flow was thus allowed at the outlet, and did occur. The turbulence quantities of the recirculating flow at the outlet are unknown, but to set a relevant turbulence level for the present case the back-flow values for 𝑘 and 𝜖 were assumed to be similar to the average of those quantities at the inlet. The background of this assumption is that the turbulence level is high already at the inlet due to the wakes of the stay vanes and the guide vanes. It is thus assumed that the increase in turbulence level is small compared with that at the inlet. It is further believed that the chosen values are of minor importance for the overall flow. For the pressure the homogeneous Neumann boundary condition is used at all boundaries. The computations are made for a complete runner with five blades. The computational domain is shown by Figure 1. A block-structured hexahedral wall-function grid was used, consisting of approximately 2 200 000 grid points.

2.1.4. Discretization Schemes

For the convection divergence terms in the turbulence equations the Gamma discretization scheme by Jasak et al. [23] was used. For the convection divergence terms in the velocity equations the GammaV scheme was used, which is an improved version of the Gamma scheme formulated to take into account the direction of the flow field. The Gamma scheme is a smooth and bounded blend between the second-order central differencing (CD) scheme and the first-order upwind differencing (UD) scheme. CD is used wherever it satisfies the boundedness requirements, and wherever CD is unbounded UD is used. For numerical stability reasons, however, a smooth and continous blending between CD and UD is used as CD approaches unboundedness. The smooth transition between the CD and UD schemes is controlled by a blending coefficient 𝛽_𝑚, which is chosen by the user. This coefficient should have a value in the range 0.2≤𝛽_𝑚≤1, the smaller value the sharper switch and the larger value the smoother switch between the schemes. For good resolution, this value should theoretically be kept as low as possible, while higher values are more numerically stable. Studies of different 𝛽_𝑚 values have been made, and the results are however more or less unaffected by the choice of 𝛽_𝑚. In the present work a value of 𝛽_𝑚=1.0 has been used. The time derivative is discretized using the Euler implicit method.

2.2. Identification of Dynamic Coefficients

To describe how the eigenfrequencies and damping of a torsional dynamic system change due to the flow, the model illustrated in Figure 2 is used. In the model, the generator is assumed to be stiff due to the connection to a rigid electric grid, and hence only the torsional motion of the turbine runner is considered. The equation of motion for this system is given by 𝐽_𝑃𝜃̈+𝐶𝜃̇+𝐾𝜃=𝑀(𝑡), (1) where 𝐽_𝑃 is the polar inertia, 𝐶 is the damping, 𝐾 is the stiffness, 𝑀(𝑡) an external moment, 𝑡 is the time, 𝜃 is the angular displacement, 𝜃̇ is the angular velocity, and 𝜃̈ is the angular acceleration. It is further assumed that the flow through a turbine will give additional inertia, damping, and stiffness to the system. With these additional coefficients the equation of motion becomes (𝐽_𝑃+𝐽_𝑃,Fluid)𝜃̈+(𝐶+𝐶_Fluid)𝜃̇+(𝐾+𝐾_Fluid)𝜃=𝑀(𝑡), (2) where 𝐽_𝑃,Fluid is the added polar inertia, 𝐶_Fluid is the added damping, and 𝐾_Fluid is the added stiffness. External moments are negligible (𝑀(𝑡)=0) in the present work. CFD is used to identify the added coefficients from the torque of the turbine runner. Rewriting the moments due to the flow to 𝑇^′(𝑡)=−𝐽_𝑃,Fluid𝜃̈−𝐶_Fluid𝜃̇−𝐾_Fluid𝜃, (3) where 𝑇^′(𝑡) is the total torsional moment due to the flow, and inserting this into (2) yields 𝐽_𝑃𝜃̈+𝐶𝜃̇+𝐾𝜃=𝑇^′(𝑡). (4) To solve 𝑇^′(𝑡), the forces and moments from the CFD-simulations are calculated at each time step. The force on a control volume boundary face is given by 𝐹_face,𝑖=𝑝_face,𝑖𝐴_face,𝑖𝑛_face,𝑖, (5) where 𝑝_face,𝑖 is the pressure of the face, 𝐴_face,𝑖 is the area of the face, and 𝑛_face,𝑖 is the normal vector of the face. The moment of the centre of gravity of the runner at a face is 𝑀_face,𝑖=𝐹_face,𝑖𝑟_face,𝑖, (6) where 𝑟_face is the radius from the centre of gravity to the face. The total moment is calculated as 𝑀=^𝑛_𝑖=1𝑀_face,𝑖, (7) where 𝑛 is the number of faces. The torque is obtained as a scalar product of the moment and the direction vector of the shaft 𝑇(𝑡)=𝑀𝑛_𝑦. (8) During steady conditions the torque is constant in order to provide a constant power to the generator. In case of unsteady conditions, the torque can be written as 𝑇(𝑡)=𝑇_mean+𝑇^′(𝑡), (9) where 𝑇_mean is the constant part of the torque. In the present work the rotational speed of the turbine runner is prescribed in order to determine the dynamical coefficients of the turbine runner due to the flow. The angular displacement of the runner is given by 𝜃=Ω𝑡+acos(𝜗𝑡)=Ω𝑡+𝜃^′, (10) where Ω is the constant angular velocity, 𝑡 is the time, 𝑎 is an amplitude, 𝜗 is a frequency of the prescribed runner oscillation, and 𝜃^′ is the oscillating part of 𝜃. Below, we are only interested in the oscillating part, where 𝜃^′=acos(𝜗𝑡), (11) gives the velocity

Figure 2 

The mechanical model of a torsional dynamic system.

𝜃̇^′=−𝑎𝜗sin(𝜗𝑡), (12) and the acceleration 𝜃̈^′=−𝑎𝜗²sin(𝜗𝑡). (13)

Inserting (11), (12), and (13) into (3) results in an equation for the fluctuation of the torque 𝑇^′(𝑡)=𝑎𝜗²𝐽_𝑃,Fluidcos(𝜗𝑡)+𝑎𝜗𝐶_Fluidsin(𝜗𝑡)−𝑎𝐾_Fluidcos(𝜗𝑡). (14) This can be written as 𝑇^′(𝑡)=𝑇_Ampcos(𝜗𝑡−𝜙)=𝑇₁cos(𝜗𝑡)+𝑇₂sin(𝜗𝑡), (15)

where 𝑇_Amp is the amplitude of the torque, 𝜙 is the phase angle, and 𝑇₁ and 𝑇₂ are the cosine and sine components of the amplitude. Then the additional damping due to the fluid can be identified as

𝐶_Fluid=𝑇₂𝑎𝜗, (16) and the additional stiffness and polar inertia due to the fluid can be identified by solving 𝑎𝜗²𝐽_𝑃,Fluid−𝑎𝐾_Fluid=𝑇₁ (17) for two simulations with different values of 𝜗.

The eigenfrequency of (2) can now be solved as

𝜔_𝐷=𝐾+𝐾_Fluid𝐽_𝑃+𝐽_𝑃,Fluid−(𝐶+𝐶_Fluid)²4(𝐽_𝑃+𝐽_𝑃,Fluid)², (18) and the corresponding damping ratio is

𝜁=𝐶+𝐶_Fluid2(𝐽_𝑃+𝐽_𝑃,Fluid)(𝐾+𝐾_Fluid)/(𝐽_𝑃+𝐽_𝑃,Fluid). (19)

3. Results

In Figure 3 the torque is shown as a function of time for one of the simulated cases. The amplitude of 𝑇₁/𝑎 in (17) is presented as a function of perturbation frequency in Figure 4. The perturbation amplitude is 𝑎=4.0×10⁻⁶rad for all simulations and is selected in the area where torque/angular velocity is linear and the value is selected in order to separate the response from numerical noise. One can see that it is difficult to identify the coefficients as stated in (17). There are two possible explanations to this: the coefficients depend on frequency and the stiffness is probably small due to the incompressible fluid. The stiffness is therefore assumed to be negligible (𝐾_Fluid=0 in (17) in the analysis below. The added polar inertia is presented in Figure 5 and the added damping in Figure 6.

Figure 3 

The torque as a function of time for one of the simulated cases (rotational speed is 72rad/s and the oscillating frequency is 1809rad/s).

Figure 4 

Identification of the coefficients of (17), together with a curve-fitted function (rotational speed is 52rad/s).

Figure 5 

Additional polar inertia as a function of perturbing frequency and operating condition.

Figure 6 

Additional damping as a function of perturbing frequency and operating condition.

The later coefficients are added to the mechanical system, that is, (2). The polar inertia of the mechanical system is 𝐽_𝑃=1.57Nms², the damping is 𝐶=0Nms, and the stiffness is 𝐾=49000Nm. In Figure 7 the reduced eigenfrequencies (18) and in Figure 8 the damping ratio (19) due to the flow for such a fluid-mechanical system are presented and the influence of the different coefficients is illustrated.

Figure 7 

Reduction of the eigenfrequency (the eigenfrequency of the mechanical system is 1) due to the flow through the turbine. The “undamped” markers represent the effect of an added polar inertia alone.

Figure 8 

Additional damping due to the flow through the turbine (the damping of the mechanical system is zero).

4. Discussion

Both added polar inertia and damping have a significant effect on the eigenfrequency of the mechanical system. The added polar inertia decreases the eigenfrequency 3-5% for all cases (see Figure 7). Concerning the damping, an additionally decrease of the eigenfrequency of 5-60% is observed (see Figure 7). One can see that both damping and polar inertia increases for off-nominal speed and with frequency. Recent research by Liang et al. [18], and Rodriguez et al. [19] has shown that the eigenfrequencies are reduced by 10-39% for a nonrotating Francis runner in still water. The effect of added inertia in these papers are significantly higher than the case of nominal operating condition in the present work and the authors observe no strong effect of damping. An explanation to the difference between the present study and the earlier work is the dependency of frequency for both added inertia and damping and that the present work includes the turbine flow.

Iso-surfaces are here used to illustrate the difference between the different operating conditions. Figure 9 shows iso-surfaces of regions where the turbulent kinetic energy is high. In Figures 10, 11, and 12 smearlines at the blades are presented in order to see the details of the flow.

(a)

(b)

(c)

Figure 9 

Iso-surface of turbulent kinetic energy, (a) 52rad/s, (b) 62rad/s, and (c) 72rad/s.

Figure 10 

Smearlines and velocity vectors for 52rad/s.

Figure 11 

Smearlines and velocity vectors for 62rad/s.

Figure 12 

Smearlines and velocity vectors for 72rad/s.

The difference in the rotating speed results in different flow conditions for the different operating conditions. The guide vane angle is equal for all cases. Hence, the angle of attack at the leading edge of the runner blades is changed when changing the rotational speed. The tipclearance flow from the pressure side to the suction side is increased when the rotational speed is reduced. For high rotational speeds there is also a tip vortex at the runner blade pressure side due to the unfavorable angle of attack close to the tip. The tip vortex flow is the reason to the high turbulent kinetic energy near the tipclearence, which is shown in Figure 9. Figure 9 also shows high turbulence kinetic energy in the flow stagnation at the leading edges of the runner blades, and in separation regions. A major difference in the level of turbulence kinetic energy can be found below the runner cone in the recirculation region. The significant differences of the flow field for the different cases are also illustrated by the smearlines in Figures 10, 11, and 12. Figures 10 and 12 show a large non-axisymmetric recirculation area below the cone. The wakes below the runner vanes are also shown on the cone as well as the tipvortex flow. Figure 11 shows a small axisymmetric recirculation area below the cone.

Recent research of added mass of a cylinder by Wang et al. [24] has shown that the added mass is dependent on the velocity around a cylinder. The same effect is suspected in the present study, where the flow velocity differs between the cases.

5. Conclusions

The added polar inertia and damping due to the hydraulic system significantly affect the mechanical system. This results in a reduced eigenfrequency of 5-65% and an increase in the damping of 30-80%. It is further concluded that the added coefficients are dependent on the turbine load and oscillating frequency. A change in the system properties of the mechanical system is important to consider in design and operation. Future studies should include experimental verification of the results in the present work.

Nomenclature

𝜃:	Angular displacement (rad)
𝜔𝐷:	Damped natural frequency (rad/s)
𝜁:	Damping ratio (−)
𝜗:	Prescribed frequency (rad/s)
Ω:	Rotational speed (rad/s)
𝑛face,𝑖:	Normal vector at one face (−)
𝑝face,𝑖:	Pressure one face (N/m2)
𝑡:	Time (s)
𝐴face,𝑖:	Area of one face (Nm)
𝐶:	Damping (Nms/rad)
𝐶Fluid:	Added damping (Nms/rad)
𝐹face,𝑖:	Force on one face (N)
𝐽𝑃:	Polar moment of inertia (kgm2)
𝐽𝑃,Fluid:	Added Polar moment of inertia (kgm2)
𝐾:	Stiffness (Nm/rad)
𝐾Fluid:	Added stiffness (Nm/rad)
𝑀(𝑡):	External moment (Nm)
𝑀face,𝑖:	Moment at one face (Nm)
𝑇′(𝑡):	Total torsional torque due to flow (Nm)
𝑇1,2:	Sine and cosine components of the torque (Nm)
𝑇Amp:	Amplitude of the oscillating part of the torque (Nm)
𝑇Mean:	Constant part of the torque (Nm)

Acknowledgments

The research presented in this paper has been carried out with funding by Elforsk AB and the Swedish Energy Agency through their joint Elektra programme and as a part of the Swedish Hydropower Centre (SVC) (http://www.svc.nu/). SVC has been established by the Swedish Energy Agency, Elforsk, and Svenska Kraftnät together with Luleå University of Technology, The Royal Institute of Technology, Chalmers University of Technology, and Uppsala University. Computations have been carried out with support from the Swedish National Infrastructure for Computing on the Hive and Ada clusters at C3SE, Chalmers University of Technology.