Experimental and Numerical Simulations Predictions Comparison of Power and Efficiency in Hydraulic Turbine

Abstract

On-site power and mass flow rate measurements were conducted in a hydroelectric power plant (Mexico). Mass flow rate was obtained using Gibson's water hammer-based method. A numerical counterpart was carried out by using the commercial CFD software, and flow simulations were performed to principal components of a hydraulic turbine: runner and draft tube. Inlet boundary conditions for the runner were obtained from a previous simulation conducted in the spiral case. The computed results at the runner's outlet were used to conduct the subsequent draft tube simulation. The numerical results from the runner's flow simulation provided data to compute the torque and the turbine's power. Power-versus-efficiency curves were built, and very good agreement was found between experimental and numerical data.

1. Introduction

Nowadays, there is a clear need for low-degrading and low-polluting energy transformation processes. In this respect, wind and hydraulic turbines have taken an important role due to the intrinsic absence of combustion in the electricity production. Along this line of thinking, this work illustrates an effort taking place in Mexico which addresses the performance of hydropower stations and, in particular, the enhancement of the efficiency of the turbine components.

In Mexico, the energy generation through hydroelectric power plants corresponds to 22.14% of the total (11,094.90 MW) [1]. Most of these plants use turbines from the 1950s, and their current operating conditions are somehow different from those specified by the original design. It is, thus, necessary to assess the present conditions to look for possible modifications including optimization and corresponding refurbishing. In order to do this, it is essential to measure the water discharge. This task is not easy to perform in field conditions and is particularly complex in water power plants with no access to the penstock from outside and with no suitable measurement installations prepared during construction. To solve this problem, the technique known as Gibson’s method, by means of which the flow rate is obtained by integration in time of the pressure difference variation between two cross sections of the penstock, can be used.

Although Gibson’s method has been recently improved [2], the handling of practical aspects still requires the installation of experimental apparatus (pressure transducer, gauges, holes, etc.) and, obviously, to operate the turbine in test rather than in production mode for a period of time. Considering this type of constraint, it is useful to apply computational fluid dynamics as an attempt to lower the cost of gathering experimental data.

In this work, both aspects have been considered. First, on-site measurements of flow rate and power were conducted at the hydropower plant in the state of Oaxaca in Mexico, then, simulations were performed using the commercial software ANSYS CFX to reproduce aspects of the real-life situation.

2. Gibson’s Method for Flow Measurements

This method, devised to measure the rate flow in a hydraulic turbine, is based on the water hammer phenomenon taking place in a closed pipe. It was introduced by Gibson [3], who benefitted from the work conducted by Jukowsky in 1898 concerning the water-hammer theory. Gibson’s method measures a static pressure difference between two cross-sections of the penstock as a result of momentum variation caused by a quick closing of the wicked gates of the turbine. The flow rate is then obtained by integration, within a proper time interval (~10 s) [4–6]. Figure 1 illustrates the scheme of measuring in the penstock.

Figure 1 

Penstock and measuring sections in Gibson’s method.

In order to derive a relationship for computing the flow rate 𝑄, the following hypothesis are considered:(i)a closed pipe with a flow section area 𝐴 that may change along its length, (ii)initial constant velocity and pressure fields between two given sections along of the penstock, (iii)that the water flow be completely stopped when the water hammer occurs, (iv)constants density and constant flow section during the water hammer.

Based in these assumptions, the relation between the parameters of the one-dimension unsteady flow between two selected sections of the pipe can be described using the energy balance equation𝜌𝑣²₁2+𝑝₁+𝜌𝑔𝑧₁=𝜌𝑣²₂2+𝑝₂+𝜌𝑔𝑧₂+Δ𝑝_𝑓+𝜌^𝐿₀𝜕𝑣𝜕𝑡𝑑𝑥. (1) Considering 𝑄=𝑣𝐴, the above equation becomes [6]:𝜌𝑄²2𝐴²₁+𝑝₁+𝜌𝑔𝑧₁=𝜌𝑄²2𝐴²₂+𝑝₂+𝜌𝑔𝑧₂+Δ𝑝_𝑓+𝜌^𝐿₀𝑑𝑄𝑑𝑡𝑑𝑥𝐴(𝑥). (2) Integrating this equation between time interval (𝑡𝑜,𝑡𝑘), and grouping terms yields:𝑄₀=1𝜌𝐶^𝑡𝑘_𝑡0Δ𝑝(𝑡)+Δ𝑝_𝑑(𝑡)+Δ𝑝_𝑓(𝑡)𝑑𝑡+𝑄_𝑘, (3) where 𝑄₀=𝑣₀𝐴 represents the flow rate at steady state before closing the wicked gates, Δ𝑝=𝑝₂+𝜌𝑔𝑧₂−𝑝₁−𝜌𝑔𝑧₁ indicates pressure difference between sections 1 and 2, Δ𝑝_𝑑=𝜌𝑄²/2𝐴²₂−𝜌𝑄²/2𝐴²₁ Denotes dynamic pressure difference between sections 1 and 2, 𝑄_𝑘: Flow rate measured in the leaking test (ideally should be 0), and 𝐶 = specifies the geometrical modulus of the penstock segment of a length 𝐿 given by𝐶=^𝐿₀𝑑𝑥𝐴(𝑥). (4)

For a penstock with constant diameter 𝐶=𝐿/𝐴.

The leakage flow 𝑄_𝑘 is determined by the empirical equation𝑄_𝑘=𝜇𝐴_𝑠2Δ𝑝_𝑘𝜌. (5)

3. Measurement Instrumentation

To carry out the measurements, two separate systems of signal acquisition were used [4]; the first consisted of a laptop computer and a 16 bits signal acquisition conditioning module card. The second consisted of a National Instrument PXI1010 system with signal acquisition module SCXI 16-bit card. The two pressure transducers used for testing were PS-2p (Vist, Gdañsk, Poland): range 0–700 kPa (relative) with accuracy of 0.1% and temperature compensated output current of 4–20 mA calibrated by the manufacturer and on site. In the case of the transducers of section  2-2, they were checked for calibration before the test. Figure 2 shows the scheme of the sensing points on penstock, and Figure 3, the waterproof manifold installed, with the pressure transducer, for measurements.

Figure 2 

Scheme of sensing points.

Figure 3 

A waterproof manifold with details of the absolute pressure transducer location.

The signal acquisition system, in addition to capturing the behavior of the pressure, also recorded the opening of the wicked gates, the active power, and the level of tailwater. The level in the tailwater was measured manually and compared with indications from the pressure measured in Section  2-2; these measurements were necessary to calculate the efficiency of the turbine.

The data recording was made with the sampling frequency of 500 Hz, and the files are prepared in ASCII format with 100 Hz frequency. The flow rate was calculated with the program GIB-ADAM developed in the Szewalski Institute of Fluid-Flow Machinery in Poland [2, 6]. Each measurement required a quick closing of the distributor and the generator being connected to the network and in stable condition.

4. Measurements Results

The flow rate was determined from the recorded pressure time histories using a GIB-ADAM [4] program. Figure 4 shows a sample of record of the time history of turbine guide vane closure and the differential pressure between the measuring sections of the penstock.

(a)

(b)

Figure 4 

Guide vane closure percentage and change in pressure between section  2-2 depending on time 𝑡.

Once made, the 8 rejects of charge and recorded all the data. Flow measurements were conducted in five loading conditions: 25%, 50%, 75%, 85%, and 100% of load. A summary of measurements is show on Table 1.

Table 1 

Summary of measurements.

5. Governing Equations for Flow Simulation

For simulations a three-dimensional incompressible flow was considered, with constant properties and isothermal at 25°C·k-ε turbulence model standard. The flow inside the suction tube is described by the Navier-Stokes-Reynolds averaging for the k-ε.

Continuity:
𝜕𝜌𝜕𝑡+∇⋅𝜌𝐔=0, (6) where 𝜕∇=,𝜕𝜕𝑥,𝜕𝜕𝑦,𝜕𝑧𝐔=𝐔+𝐮, (7) where 𝐔=(1/Δ𝑡)∫_𝑡^𝑡+Δ𝑡𝐔𝑑𝑡.

Momentum:
𝜕𝜌𝐔𝜇𝜕𝑡+∇⋅𝜌𝐔⊗𝐔−∇⋅_eff∇𝐔=∇𝑝^′+∇⋅𝜇_eff∇𝐔^𝑇+𝐁, (8) where B is the sum of efforts, 𝜇_eff is the effective viscosity for turbulence, and 𝑝^′ is the modified pressure gives by 𝑝^′2=𝑝+3𝜇𝜌𝑘,_eff=𝜇+𝜇_𝑡, (9)𝑘-𝜀 model assumes that the turbulent viscosity is related to the turbulent kinetic energy and dissipation via the relation 𝜇_𝑡=𝐶_𝜇𝜌𝑘²𝜀. (10) The values of 𝑘 and 𝜀 come directly from the transport equations for turbulent kinetic energy and dissipation range of turbulence 𝜕𝜌𝑘𝜇𝜕𝑡+∇⋅𝜌𝐔𝑘=∇⋅𝜇+^𝑡𝜎_𝑘∇𝑘+𝑃_𝑘𝜕−𝜌𝜀,𝜌𝜀𝜇𝜕𝑡+∇⋅𝜌𝐔𝜀=∇⋅𝜇+^𝑡+𝜀𝜎𝜀∇𝜀𝑘𝐶_𝜀1𝑃_𝑘−𝐶_𝜀2𝜌𝜀, (11) where 𝐶_𝜀1,𝐶_𝜀2,𝜎_𝑘, and 𝜎_𝜀 are constants and 𝑝_𝑘 is the production of turbulence due to the forces of viscosity which is modeled as follows: 𝑃_𝑘=𝜇_𝑡∇𝑈⋅∇U+∇𝑈^𝑇−23∇⋅𝑈3𝜇_𝑡∇⋅𝑈+𝜌𝑘+𝑃_𝑘𝑏. (12) To incompressible flow, ∇⋅𝑈 is small, and the second term at the right from equation does not contribute significantly in production.

5.1. Numerical Analysis

The commercial software used employs the finite volume numerical method for solving the RANS equations and 𝑘-𝜀 standard model. To perform the numerical simulations, we used an unstructured mesh of 1,759,261 elements for the runner and a structured mesh of 1,024,956 elements for draft tube. Simulations were run on an HP Pavilion PC with two Intel 1.66 GHz each, with 2 GB of RAM.

To perform the numerical simulations, measurement results in Tests 1, 3, 5, 6, and 8 showed in Table 1 were used.

The selected convergence criteria of the numerical simulations was RMS (root mean square) normalized values of the equation residuals with a value of E-006.

5.2. Boundary Conditions

In the runner, the boundary conditions were stipulated as follows:(i)inlet: defined as Mass flow inlet,(ii)turbulence intensity: 𝑘-𝜀 from spiral case,(iii)outlet: defined as Static Pressure (measured at the inlet of draft tube),(iv)turbulence intensity: zero gradient,(v)blade, shroud and hub: were defined as wall with no movement and as a smooth surface.

6. Results of Simulation

Before obtaining the final results of the numerical simulation, we carried out a mesh independence analysis called grid convergence index (GCI for its acronym) to estimate the percentage of error in the solution. According to the work done by Roache [7], this is calculated by the following equation:3|𝜀|GCI=𝑟^𝑝−1, (13) whereℎ𝑟=²ℎ₁,𝑓𝜀=²−𝑓₁𝑓₁, (14) where 𝜀 is the relative difference between the results of torque (a parameter selected, in this case, the most important to calculate the power of the runner) for the finest 𝑓₁ and 𝑓₂ for the coarsest, 𝑟 is the ratio between the size of the element of finest mesh (ℎ₁) and coarsest (ℎ₂), and 𝑝 is the order of the method used for the reference solutions; in this case, it was used second order. Table 2 shows the results of the analysis.

Table 2 

Analysis GCI.

Considering the processing time and storage volume of the files generated by numerical simulation results, to obtain the final results mesh no. 3 was selected, which has a convergence index that is reasonable (according to references, it is ideal when it is less than 1.5)

6.1. Runner

From the simulations on the runner, torque was obtained, with which it proceeds to the calculation of mechanical power and compare it with the measurements. The results are presented in Table 3 and Figure 5.

Table 3 

Comparison between numerical and experimental data.

Figure 5 

Efficiency curves.

At the point of 91.16% efficiency, special attention was taken on the results of simulations, which should be analyzed carefully, as these results represent the maximum efficiency. One of the most important issues is the cavitation; in Figure 6, pressure contours do not show low pressure points over the blade that could produce it.

(a)

(b)

Figure 6 

(a) Pressure contours at 76.14 m3/s, (b) Pressure contours at one blade.

Figure 7 shows the blade loading, and it can be seen that there are no sudden changes in pressure around the blade, as well as negative pressure, which provides evidence that there is no possibility that air bubbles could form and cause cavitation.

Figure 7 

Blade loading in midspan, at 76.14 m3/s.

In Figure 8, the velocity vectors show that the flow does not produce eddies when passing through the runner.

Figure 8 

Velocity vectors at midspan.

6.2. Draft Tube

From the results of simulations on the runner, the boundary condition at the draft tube inlet for different measured flow rates was obtained; thereby, the behavior in this component was discovered, which is significant because it is very difficult to make direct measurements.

Figures 9 and 10 show streamlines and pressure contours in the draft tube at different flow rates. The behavior of the flow in this component is expected as observed in the literature [1, 8–10]. The draft tube has the function of discharging water from a turbine, in addition to acting as a device for energy recovery, helping to improve the overall performance of the unit. It can also allow the water level downstream to be lower or higher than the equatorial plane of the turbine, according to the installation needs. The draft tube, because of its divergent shape, produces a slowdown velocity of the water leaving the turbine, converting fluid kinetic energy into pressure energy (Figure 10).

(a)

(b)

Figure 9 

Streamlines at different flow rates (a) 76.14 m3/s and (b) 89 m3/s.

(a)

(b)

Figure 10 

Pressure contours in cross-section at different flow rates (a) 76.14 m3/s and (b) 89 m3/s.

In Figure 9, it can be observed that as the flow moves into the tube, their behavior is changing and the speed at the elbow down to about half the initial velocity. Besides, because of the pier on the draft tube, there is a separation of the flow, and the velocity decreases further in the side where the vortex produced by the runner is directed.

With these results, we will propose amendments to this component geometry to increase efficiency.

7. Conclusions

The mechanical power and efficiency obtained from the simulations have good agreement with those obtained on-site measurements.

From the results at the maximum efficiency, adverse conditions in the flow will not appear, so the turbine can work properly at this head, and thereby improve the performance of flow in the unit of study.

Therefore, these results can be considered reliable enough for the speeds at the exit of the runner and could be used in the boundary condition at the entrance of the draft tube.

The results show that the flow in draft tube is sufficiently approximate to the actual flow, according to revised bibliography.

Notation

𝜌:	Density
𝑧1,2:	Hydrometric levels
Δ𝑝𝑓:	Pressure drop caused by friction losses
𝑄0:	Flow before the closing of the wicked gates
𝑄𝑘:	Leakage flow
𝐶:	Geometrical modulus of the penstock segment
𝐿:	Length of the segment of the penstock
𝐴:	Transversal area of the penstock
𝜇:	Ratio of slack flow (0.65–0.7)
𝐴𝑠:	Clear area of the wicked gate
Δ𝑝𝑘:	pressure difference between inlet and outlet of the distributor
𝜂𝑥:	Recovery ratio
G:	Acceleration of gravity
Re:	Reynolds Number
𝐷:	Diameter
𝜈:	Kinematic viscosity
𝑚̇:	Mass flow
𝜇:	Dynamic viscosity
𝐶𝜇:	Empirical constant
𝐾:	Kinetic turbulent energy
𝜀:	Turbulence dissipation rate
𝑈:	Axial velocity
𝑉:	Radial velocity
𝑊:	Tangential velocity
𝑃:	Static pressure
𝜇eff:	Effective viscosity
𝜇𝑇:	Turbulent viscosity
𝐶𝜀1:	Turbulence constant
𝐶𝜀2:	Turbulence constant
𝜎𝑘:	Turbulence constant
𝜎𝜀:	Turbulence constant
𝑃𝑘:	Turbulence production by viscosity forces
𝜏𝜔:	Shear stress
𝜅:	Von Karman constant
GCI:	Grid Convergence Index
𝜀:	Relative error
𝑅:	Element size ratio
𝐶𝑝:	Pressure coefficient
𝑝:	Average static pressure
V:	Volume
𝜔:	Turbulence frequency
𝜂𝑡:	Efficiency
𝑄𝑡:	Flow
𝑃𝑚:	Mechanical power.

Acknowledgments

The authors are grateful to the Federal Electricity Commission (CFE) for the facilities provided for the present study and to the National Council of Science and Technology (CONACYT) for the scholarship granted to the Ph.D. student.