Most of turbine flow simulations are performed in a steady state operating
regime, since the engineers are mainly interested in mean flows and the
effort required to perform such a simulation is significantly lower
compared to an unsteady flow simulation. However, the turbine has to
be designed to withstand extreme loads under normal, emergency or
catastrophic transient regimes and the designers have to predict the
turbine behaviour under such conditions and design
protective measures against them. In present research, we used CFD method
to simulate the transient of the sudden turbine shut down. The example was
the "Vrhovo" HPP on Sava river, Slovenia (Fig. 1). The turbine was a bulb type,
double regulated, where the term refers to ability to control the angle
of distributor vanes and runner blades during the turbine operation (Fig. 2).
Fig 1: The Vrhovo HPP flow passage geometry
Fig 2: Double regulation
During sudden load rejection, the runner is rotating, the guide vanes
of a distributor as well as the runner blades relatively to the hub
are moving simultaneously during the simulation period. If one wants to
simulate these movements by CFD method,
the computational grid has to follow the boundaries and stay valid. Complex
distributor vanes' and runner blades' movements prevent application of
simple computational grid transformations for description of boundary
movements. Therefore, a method
was developed which preserves the consistency of the numerical grid
throughout the interval of complex geometrical changes.
Numerical method
The water flow in a turbine flow-passage is considered viscous and turbulent.
It may be fully described by the 3D time-averaged Navier-Stokes equations,
which are accompanied by the two-equation turbulence model. In this study
we have utilized a CFD flow-solver code ICCM COMET. The transient term
of the NS equation was discretized using the first order Euler
implicit scheme. The approximation of the convective term of the
NS equations has an important impact on the accuracy and stability of
the numerical solution. The upwind differencing scheme (UD) and central
differencing scheme (CD) were alternatively used.
The turbine flow-passage calculation domain is bounded by walls, where
standard no-slip boundary condition was applied. The surfaces of distributor vanes
move according to the prescribed law, whereas the computational domain
comprising runner is rotating about the main turbine axis. During the steady
operating turbine regime the runner rotational speed is held
constant by the turbine governor according to the electrical power consumption.
In the case of sudden load rejection, the runner speed becomes dependent solely
on the water flow field in the turbine flow passage and mechanical properties
of the system (mass moment of inertia of rotating parts). The runner
acceleration can be derived from the Torque acting on runner and the mass moment
of inertia. Torque is generated by fluid-flow-induced hydraulic forces and may
be extracted from the CFD solver by integrating the torque contributions at
each control volume surface on the whole runner and hub surface.
Equation was coupled with CFD solver within a code which automatically
generated new position of computational grids for each time step.
Fig 3: Boundary conditions
At the inlet and at the outlet of the turbine the total and static pressure
were prescribed respectively. The gross head was selected as the inlet total
pressure. The static pressure at the outlet was set to zero. Consequently
the turbine discharge can be calculated as a function of geometrical position
of distributor vanes and runner blades. Fig. 3 shows the boundary
condition setup.
Grid generation
Grid generation for all turbine components has been acomplished by the
proprietary tool "gridedit" (see references). Whole turbine conduit represented
the computational domain. Three levels of grid density were genereated. The most
dense grid comprised approx 400.00 cells. Fig 4 shows grid of density level 2
with approx 150.000 cells.
Fig 4: The grid of the whole water passage
The geometric changes of the
domain were approximated by generating a number of intermediate positions
in the preprocessing step ("click" positions, see Fig 5). In the course
of calculation, the grids corresponding to intermediate geometric positions
were switched on in appropriate moment according to a prescribed scenario,
which drove the temporal events (Fig. 6).
Fig 5: Scenario of grid changes of the distributor
Fig 6: The temporal events
The (linear) interpolation was used to generate different
grid in every time step. Every node coordinate was interpolated between
successive "click" grids for an appropriate time step prior to flow
calculation. With this technique, a "continuous" movement of computational
domain boundaries was ensured in all temporal steps of the numerical procedure.
Fig. 7 shows part of the distributor grid at selected temporal points.
Fig 7: The temporal changes of the distributor
grid (82%, 40%, 2.5% opening, respectively)
The results
To evaluate our method the results obtained by simulation were compared
with the experimental data on the prototype. Measurements were performed
on the Vrhovo HPP, Slovenia. It is a run-off river type power plant and
consists of three bulb turbine units (T1, T2 and T3), each of rated power
P = 11 MW at rated head H = 7.3 m and rated discharge Qr = 166 m3/s.
The turbine rotational speed is n = 1.786 s-1. Turbine sudden load
rejection experiment was carried out on turbine T2. The following
characteristics were measured: turbine rotational speed n, axial force
Fa and pressure p at a selected point MP (Fig. 8). The positions of
distributor and runner were measured indirectly by potentiometers, mounted
on servomotors and are denoted as y_gv and y_rb, respectively.
Fig. 8: The location of sensors and experimental setup
The steady operating regime
The initial distributor vane and runner blade servomotor strokes were
y_dv = 82 % and y_rb = 48 %, respectively. For the sake of simplicity
the whole flow field was initialized with constant speed in direction,
parallel to the main turbine axis, while the runner immediately began to
rotate with constant speed. At such initial conditions, it took approximately
1500 time steps (4.2 revolutions of the runner) for the discharge and
torque to become nearly steady. Fig. 9 shows this
initial period of time dependent simulation. Before "triggering" the transient
event, the simulation was running for about 10,000 time steps or 28
revolutions of the runner.
Fig 9: Turbine characteristics in the initial period of simulation
The transient operating regime (sudden load rejection)
Fig. 10 shows comparison between simulated rotational speeds for all
computational cases and the experiment. We have performed a number of
calculations variing several numerical parameters, such as grid density and
the discretuzation scheme. g_UD_1 denotes upwind differencing with low grid
density, whereas g_CD_3 denotes central differencing with the highest grid
density. It can be noticed that UD simulations
yield lower rotational speeds during a transient event. The deviation
of "the best" simulation from the experiment is lower than 5% in the whole
observed time interval (30 seconds of real time).
Fig 10: Comparison between simulated rotating speeds and the experiment
(y_dv = distributor vane servomotor stroke; y_rb = runner blade servomotor
stroke; g* = computational case).
Fig. 11 shows comparison between simulated and experimentally obtained axial
forces. The deviation of simulated axial force from the experiment is lower
than 5% in nearly whole observed time interval.
Pressure at the selected point MP (depicted in Fig. 8) is compared in Fig. 12.
Because of large measured pressure oscillations (passing of the runner
blades near the sensor) the diagram also shows the averaged (measured)
pressure.
Fig. 11: Comparison between simulated and experimentally obtained axial forces
Fig. 12: Comparison between simulated and measured pressure at the
point MP (between distributor and runner).
There is a good qualitative agreement between all simulated and measured
pressures at the point MP. The best agreement with experiment
can be found in the computational case g3_CD. The simulated characteristics
approach the measured one by increasing the numerical accuracy
(g1 ... g3, UD ... CD). The deviation may partially be attributed to the
position of the sensor in simulation, which is very close to sliding interface.
Animation of the turbine pressure during shut down
Flash player 8.0 required.
Additional documents
[1] KOLSEK T, DUHOVNIK J, BERGANT A. Simulation of unsteady flow and runner
rotation during shut-down of an axial water turbine. J. Hydraul. Res.,
2006, vol. 44, nr. 1, pp 129-137.
[2] BERGANT A, KOLSEK T. Developments in bulb turbine three-dimensional
water hammer modelling. V: Proceedings of the 21st IAHR Symposium on Hydraulic
Machinery and Systems. Lausanne, 2002, pp. 565-571.
[3] KOLSEK T, SUBELJ M, DUHOVNIK J. Generation of block-structured grids in
complex computational domains using templates. Finite elem. anal. des..
[Print ed.], 2003, vol. 39, pp. 1139-1154.
[4] KOLSEK T. Design of specific shapes in hydraulic conduits, PhD thesis,
Faculty of Mechanical Engineering, UNI-Ljubljana, June 2002
(available in
pdf in slovene (17 MB) and
pdf in english (13 MB))
[5] SUBELJ M. Grid generator for complex 3d domains, BSc thesis,
Faculty of Mechanical Engineering, UNI-Ljubljana, 2001 (in slovene only)
