T RANSIENT 3D CFD SIMULATION OF A P ELTON TURBINE – A STATE - OF - THE - ART APPROACH FOR P ELTON DEVELOPMENT AND OPTIMISATION

,


NOMENCLATURE
cross-section (of the jet) [

INTRODUCTION
The application of computational fluid dynamics (CFD) has a long history in investigating, designing, and optimising hydraulic fluid machinery. This is especially true for Francis and Kaplan turbines, which can largely be transferred to a steady-state problem. Success and experience with the simulation of flow in Kaplan and Francis turbines subsequently led to the technique and know-how also being applied to components of Pelton turbines [1].
Due to the entirely different operating principle of Pelton turbines, the transfer of know-how took place in various stages. The research started with the components that are easier to simulate, namely the distributor and the nozzle, which have a primarily single-phased flow. The flow in the distributor, including the nozzles, was investigated in a singlephase flow setup in the early days of numerical simulations in Pelton hydraulics in [2]. Conclusions were made about possible losses and disturbing influences on the formation of the jet. Downstream of the nozzles, the flow is two-phased, which had to be considered in the modelling. Muggli and Zhang [3] carried out a two-phase flow simulation of the jet formation, based on a planar numerical grid using rotational periodicity, and compared the results with experimental data. However, the situation turns out to be much more complex when investigating the Pelton runner. As of the working principle, significantly more complex flow phenomena occurred than is the case with Kaplan or Francis turbines. For example, the transient behaviour has to be considered when simulating the jet-bucket interaction and the bucket flow. In addition, high gradients of pressure, velocity and the two phases to be considered, water and air, further complicate numerical investigations. Splash water and droplets in the turbine housing cause ventilation losses due to interaction with the runner and the jets, which can hardly be tracked by the grid-based numerical simulations. In the early stage, CFD simulations of the jet-bucket interaction were carried out using the impingement of a stationary bucket, as shown in [4]. With the increasing performance of modern computers and data storage, the complexity of the numerical investigations carried out on Pelton turbines also progressed. This is increasingly evident in numerical experiments including the entire flow path of the Pelton turbine, as shown in [5].
Progress in this area enables an increasing understanding of the flow processes. It also reduces the dependence on model tests in the development process, even if these cannot be entirely dispensed. As a result, numerical simulation of Pelton turbines is no longer only used for theoretical analysis. On the contrary, it has already taken the path to a development and optimisation tool, not only for the distributor and the nozzles but also the turbine runner, as Židonis has shown in [6].
Many published studies on numerical investigations of Pelton turbines use a single-nozzle turbine model, e.g. as in [6][7][8]. However, suppose the aim is to investigate the overall performance of multi-nozzle Pelton turbines. In that case, it is sometimes crucial to consider all nozzles in the simulation setup. Only in this way can the bucket emptying be assessed, and possible jet interference be detected, which considerably influences the overall operating behaviour.
When using a Eulerian CFD code, the number of necessary grid cells respectively nodes may increase considerably, and so does the simulation time. This detail is especially significant when combining the number of nozzles and the number of buckets of a given turbine does not allow for symmetry effects in the circumferential direction. The Pelton turbine investigated in the present paper represents such a case, with a six-nozzle turbine and a runner with 19 buckets.
As alternative to a Eulerian CFD code, Lagrangian CFD codes, such as Smoothed Particle Hydrodynamics (SPH) or Fast Lagrangian Solver (FLS), are used to analyse Pelton turbines. An overview of the application of SPH in Pelton turbines is given in [9].

NUMERICAL SIMULATION
This study's numerical modelling and simulation were based on the commercial code Ansys CFX, Release 20.1. This code is widely used in hydraulic fluid machinery, as it has provided good results in the simulation of Pelton turbines in terms of accuracy. Still, it is associated with long simulation durations and correspondingly high numerical costs [10].
The numerical study of the given Pelton turbine was split into two distinct parts. In the first part, an analysis of the Pelton distributor wa s conducted, whereby the hydraulic losses and the jet shape were of interest. Therefore, the calculations were performed in a single-and two-phase setup.
In a second step, the behaviour of the runner wa s investigated. Regarding the runner simulation, the focus was laid both on the overall performance in the form of the hydraulic efficiency as well as the generated torque and the evaluation of the jet-bucket interaction. The latter was based on the pressure distribution analysis on the runner surface and the flow situation between jet and bucket. Unlike in various publications [11][12][13], the numerical investigation did not cover a series of different operating points regarding a variation of the head or the flowrate . In this case, rather a numerically stable and most efficient procedure for the simulation of Pelton turbine runners should be developed. This procedure should enable the modelling of a wide range of turbines, regardless of the number of nozzles or buckets, and finally, it should allow for a numerically based, relative performance comparison of different runner geometries based on an identical load case.

Prototype Turbine Model
The study was based on the geometric model of a vertical six-nozzle turbine in prototype scale. The design of the distributor and the runner represent the result of a design optimisation process. As a result, the final distributor geometry incorporates changes of the general layout and the design of the bifurcations compared to the baseline design, as it is shown in Figure 1. The design of the Pelton runner was slightly modified in its main dimensions throughout the study, leading to three different runner geometries. Subject to changes were e.g. the length, width, depth of the bucket. Also, minor changes were applied to specific geometric details.

Figure 1. Tubine design, baseline (left) vs. final (right)
The nominal operating point of the turbine was given with a head of 382 m and a flowrate of 5.5 m³/s, both reached at a nozzle opening of 75%. The rotational speed was given with 600 rpm. The critical dimensions of the turbine are provided in Table 1, whereby indicates the pitch diameter of the runner, the inner width of the buckets, 0 the jet diameter and the number of buckets. The ratio / gives a number, which is widely used to characterise Pelton turbines.

Numerical Modelling of the Distributor
The modelling of turbulence in the steady-state simulations of the distributor was realised with the k-ω SST model with automatic wall function.
The high-resolution scheme in ANSYS CFX was used to model the advection terms. In addition, relaxation factors for gradient dependent advection terms were set to enhance numerical stability.
For modelling of the two-phase flow, there are two reasonable options, the inhomogeneous model and the homogeneous model. The first mentioned uses a set of conservation equations for each fluid with specific terms, which model the coupling of the different phases [14]. The higher complexity of this model with its higher number of equations results in a better performance compared to the homogeneous model [15]. On the other hand, the complexity causes an increased computational effort to solve the equations. The homogeneous model assumes that both fluids share the same flow field. Therefore, only one set of conservation equations needs to be solved with an additional equation, which determines the concentration of each fluid in the gird cells. Thus, a less demanding numerical scheme is necessary than with the inhomogeneous model. However, the homogeneous model was preferred over the more accurate inhomogeneous model due to higher numerical costs. The selected model also provided a good trade-off between the duration of the simulations and the numerical stability.
Surface tension and gravity were not considered in the modelling since the flow in Pelton turbines is determined by high Froude numbers and Weber numbers . Hence, influence of both is considerably weak on the flow in the distributor and the jet [5].
The timestep ∆ t was set to 5•10 -5 s, which provided a good balance between convergence speed and numerical stability. Additionally, double precision was selected to enhance the stability further.
The geometric model of the distributor was extended by two components. A cylinder of the length of five times the inlet diameter was added to the inlet of the distributor. Furthermore, another component was added at the outlets of the nozzles. These volumes correspond to the space in which the jets are formed. Only minor geometric simplifications were made to simplify the numerical modelling.

Figure 2. Geometry of the distributor simulation domain
The geometry shown in Fig. 2 represents the basis for the spacial discretisation. The final mesh was composed of several blocks of structured and unstructured meshes. Components including rather complex geometry elements, e.g. bifurcations, were meshed with tetrahedral and pyramidal elements. The same mesh topology was used for the leading and the trailing edge of the support ribs in the nozzles and the area of the needle tip. The other elements of the geometry were discretised with structured meshes containing primarily prismatic and hexagonal elements. Along all solid walls, thin layers of prismatic elements were placed. Fig. 3 shows the mesh of the distributor model. In regions of particular interest or with high gradients in relevant flow variables, such as the nozzle exit and the area where the phase boundary of the jet was to be expected, the grid was refined (see Fig. 4). Between the mesh components of the distributor and the nozzles general interfaces were placed. The remaining connections between structured and unstructured mesh topologies had a direct mesh connection.   Table 2. Regarding the boundary conditions, all solid walls of the model were treated as smooth, no-slip walls, without any specific surface roughness defined. The boundary surfaces of the jet domains were modelled as opening boundaries with the relative pressure rel set to zero, while the reference pressure ref is equal to the atmospheric pressure atm . At the inlet a total pressure tot,inlet was imposed, which corresponds to the given hydraulic head of the turbine.
In order to accelerate the simulation and to increase the numerical stability results of single-phase simulations of the distributor were used to initialise the flow field. At the start of the simulation, the jet domains were fully filled with air.

Numerical Modelling of the Runner
In recent publications [10] the k-ω SST model was used for turbulence modelling in Pelton runner simulations. Perrig showed in his thesis [15], that the k-ω SST model outperforms the k-ε model in terms of accuracy. Similar behaviour wa s summarised in [16]. However, due to significantly, less stable numerical behaviour in runner simulations compared to the distributor simulations, the k-ε turbulence model was applied for the runner simulations. The lack of accuracy of the chosen model is seen as a minor problem, since the overall degree of model simplifications is relatively high. Moreover, a relative comparison of different runner geometries, based on identical numerical settings, was not impaired by this.
In analogy to the distributor simulation the high resolution scheme was used to model the advection terms. The same relaxation factors as in the study of the distributor were applied. The Second Order Backward Euler scheme was selected for the time discretisation of the transient terms.
Analogous to the numerical modelling of the distributor, the homogenous model was used to model the two-phase flow in the runner. Concerning surface tension and gravity, the same assumption as in the numerical setup of the distributor were made.
The definition of the transient time step ∆ t is linked to a discrete angle of runner rotation. The size of the angle-step has been proven as crucial for some factors. A small step provides more stability to the numerical scheme, which leads also to a faster convergence behaviour and a lower overall computation time. Both effects have been observed in several simulations. Moreover, the size of the timestep influences the simulation results concerning the torque. The best simulation behaviour was found with a timestep ∆ t of approximately 4.63•10 -5 s, which corresponds to one sixth of a degree runner rotation. In the solver settings double precision was selected to further increase the numerical stability of the simulation.
The runner simulation was based on the geometric model given in Fig. 5. The simulation domain is composed of three main components: the runner, the casing and the nozzles. The former is defined as the rotating domain while the casing and the nozzles remain stationary. The runner is further split into two subdomains, which are indicated by different shadings in Fig. 5. Subsequently, these subdomains serve differently in the analysis of the runner performance. The medium grey shaded runner component is later used to analyse the performance, while the dark grey component represents some sort of a dummy.

Figure 5. Geometric model of the runner simulation domains
To simplify the model, the distributor as well as the real turbine housing were not considered for the simulation. Moreover, the symmetric layout of Pelton turbines was utalised to further reduce the compexity of the geometric model. This is a common strategy, which is seen in various publications, e.g. [8,11,14]. A non-idealised jet configuration was used. Thus, the last section of the nozzle, the nozzlehead, was included in the model to account for the effects of the fluid acceleration and the friction in the spear valve on the jet formation. The nozzle opening was set to 75%, what corresponds to the opening in the nominal operation point of the prototype turbine.
The numerical mesh of the geometry in Fig. 5 is composed of components with structured and unstructured meshes (see Fig. 6). The nozzles and the casing are represented by a structured mesh of prismatic elements. The two subdomains of the runner, with its complex geometry were split into two blocks. The first block contains the buckets and wa s therefore meshed with an unstructured, tetrahedral mesh. The second block of prismatic mesh elements adds up with the first one to complete the rotating runner subdomain. In areas with strictly aligned flow, such as in the nozzle and jet area, the mesh cells are primarily oriented according to the main flow direction. Grid refinements were placed in strategic areas, as shown for the distributor mesh. Near the walls of the runner as well as the nozzle layers of prismatic elements were placed.

Figure 6. Numerical mesh of casing and nozzles
The subdomains of the runner were meshed with two different approaches. The one, which incorporates only three consecutive buckets (medium grey shading in Fig. 5), was meshed with a comparably high mesh density. The mesh of the other runner subdomain is ra ther coarse. This approach allows to use the entire symmetric runner half to be taken into account and yet limits the size of the numerical grid to a manageable number of nodes and elements. For both subdomains, refinements were placed along the leading edge, the trailing edge as well as the cutout edge and its adjoining surfaces. A detailed view of the refined mesh of the buckets is given in Fig. 7.

Figure 7. Detail view of the mesh refinements in the cutout region of the bucket
The statistics of the mesh components of the runner simulation are given in Table 3. Quality characteristics of the runner simulation mesh are given in Table 4. Regarding the boundary conditions, the solid walls of the nozzle and the runner were treated in the same way as those of the distributor, which were modelled as no-slip walls, without any surface roughness defended. A symmetry condition was applied to the symmetry plane of the turbine model. The remaining boundary surfaces were treated with an opening boundary condition since the turbine housing was not considered. The relative pressure rel is set to zero, while the reference pressure ref was defined with the atmospheric pressure atm .
Each nozzle has an inlet surface, which is perpendicular to the corresponding nozzle axis. On these a mass flow rate inlet normal to the boundary was imposed. The link between the sliding meshes of the rotating runner and the stationary casing was treated with a transient rotor stator interface.
At the start of the first simulation a result of a distributor simulation with corresponding nozzle opening was used to initialise the stationary domain. The rotating domain was fully filled with air at the start of a simulation run. Subsequent simulations were then initialised with results of previous transient turbine simulations.

Results of the Distributor Analysis
The simulation of various Pelton distributor geometries in this study showed, that the applied changes to the geometry have led to a reduction of the hydraulic loss in the distributor and improved the flow situation upstream the nozzles. The latter is important for a good jet quality, which is indicated by a circular cross section and a homogeneous velocity profile of the jet.
A direct correlation between the jet shape and quantitative flow characteristics is difficult to determine. To overcome this, a swirl number was introduced in order to quantify the intensity of secondary flow. The definition of the is given in Eq. (1) whereby circ and ax are the velocities in circumferential and the axial direction respectively on a given plane of the cross section .
For each nozzle the velocities are referred to the corresponding nozzle axis (positive in jet direction). In the optimisation process of the distributor design, the SN was later used as an additional indicator for an improvement or a deterioration of the jet shape. If the SN decreases, an improvement of the jet shape can be expected, as smaller values indicate a lower intensity of the secondary flows, which primarily act negatively on the jet shape. Fig. 8 shows the cross section of jet number 4 on a plane about 4.25 jet diameters downstream the nozzle outlet, which was the jet with the most disturbed jet shape of all. The jet velocity, which is related to the theoretical velocity th , is plotted on the cross sections in a range of 0.75 to 1. The boundary of the cross section is determined by a water volume fraction of 50%. th is given in Eq. (2).
The left plot corresponds to the baseline geometry of the distributor, the right plot to the new design. The swirl numbers, which are evaluated at the nozzle inlet, are 0.12 and 0.036 respectively, which is equivalent to a reduction of about 70%.

Figure 8. Jet corsssection, baseline (left) vs. new design (right)
The hydraulic losses in the distributor are determined as a hydraulic head loss from the inlet to the outlet of the corresponding nozzles.
is given in Eq. 3 whereby ̅ , and ̅ , represent mass flow averaged total pressure values at the inlet and a plane, which was placed 30 mm upstream the corresponding nozzle mouthpiece.
These were lowered by roughly 16% compared to the baseline design. A comparison of the head losses, which are related to the nominal head , are given in Fig. 9. The largest share of the losses was caused by the nozzle, while the distributor contributed less than 20% to the head loss.

Results of the Runner Simulation
The results of the three runner simulations, based on different runner geometries, are analysed with a special procedure. For the analysis of the turbine performance in general, a reference bucket was defined. The reference bucket corresponds to the second bucket of the light grey shaded and comparably fine meshed runner domain in Fig. 5. Additionally, to the refence bucket, the leading and the successive bucket are considered to analyse possible, negative interactions on the runner outside surface. To calculate the runner torque, a "bucket load cycle" was defined. One load cycle corresponds to what a single bucket experiences during 360 timesteps or 60 ° of runner rotation, since the angular spacing between the nozzles in the case of the present six-nozzle turbine equals 60 °. In this period every bucket is once impinged by a jet. The arithmetic average of the torque ̅ ref over 360 timesteps or 60 ° of runner rotation, which is acting on the reference bucket, was finally used to calculate the runner torque. Additionally, ̅ ref was used as an indicator for the state of convergence of the simulation. A constant level of ̅ ref indicates a steady-state operation as well as a sufficiently converged solution. The total torque of the entire turbine runner corresponds to Eq. (4).
The efficiency is calculated from ̅ ref and the corresponding hydraulic head at the inlets as well as the total imposed flowrate at the inlets. The efficiency is given in Eq. (5), whereby represents the angular velocity of the turbine, the density and the gravitational acceleration.
The simulation of the baseline runner design in this study has shown an efficiency of 91.65%. The periodic curve of the generated torque, normalised by ̅ th for the inside, the outside and the total torque is shown in Fig. 10. The splitting of the bucket surface into an interior and an exterior surface is depicted in the upper left-hand corner of Fig. 10. The dark grey shaded surface indicates the interior, the light grey one the exterior surface. ̅ th represents a theoretical torque of a bucket based on the simulated head and the flow rate in the nominal operating point with an efficiency of 100%. ̅ th is given in Eq. (6).

Figure 10. Normalised torque (baseline design)
The torque of the primarily impinged interior surface and the exterior bucket surface add up to the torque contribution of the reference bucket. The external curve in Fig. 10 indicates the start of the load cycle with its sharp rise after roughly 390 timesteps. A low-pressure zone on the exterior surface causes the growth of the torque in this phase (see Fig. 11, left side). Perrig reported in his thesis the existence of such an external torque and explained it with the Coanda Effect [15].

Figure 11. Pressure distribution on the runner exterior surface (baseline design)
However, due to non-appropriate trailing edge angles, the same curve shows a loss following timestep 610. The impingement of the external bucket surface is caused by water, which is exiting the leading bucket. A zone of increased pressure on the exterior surface indicates the source of the loss (see Fig. 11, right side).
The bucket or the runner geometry design was modified to overcome adverse effects in an iterative process, as seen in the baseline geometry. Changes to the bucket width, the trailing edge angles, as well as the splitter and the cutout design have improved the efficiency compared to the baseline runner design. As a result, a n efficiency of 91.95% was reached.
Another set of geometry parameters of the Pelton bucket was changed consecutively in a second approach. As a result, the bucket width was slightly increased compared to the first design iteration process, and a modification of the bucket depth was realised. In addition, the positioning of the bucket on the pitch circle diameter was changed towards a smaller bucket position angle. Finally, the leading edge (splitter) and the cutout edge were further altered compared to the first design update. In the end, the numerically calculated efficiency was increased to a level of 92.30% by those measures. Fig. 11 shows the normalised torque of the second design update. Additionally, the total torque of the baseline design and the design update 1 is overlaid (dashed, green and orange curve respectively). In addition to the change of the bucket position angle, the increase of the bucket width and the bucket depth significantly contributed to a remarkable rise in the magnitude of the internal torque compared to the baseline and a small rise compared to the design update 1. On the other hand, the magnitude of the external torque in its maximum decreased compared to the baseline, which is most likely due to the changed bucket position. However, the external torque load cycle average is still higher because of the modified trailing edge angles.

CONCLUSION
While the simulation of the Pelton distributor for optimisation is quite common these days, the numerical investigation of the Pelton runner raises more questions.
As presented in this paper, the numerical approach to simulate Pelton runners, primarily of multi-nozzle turbines, should give a direction for an efficient way to analyse their hydraulic behaviour. By focusing on a single reference bucket, especially the numerical grid of the runner domain reduces to a manageable size. Furthermore, with a smaller mesh but still preserving the advantage of a high spatial discretisation in the selected regions, the computational costs and the required simulation time can be reduced significantly. Thus, numerical simulation has come one step closer to possible utilisation for the design and optimisation of Pelton runners.
Within the approach presented, three different Pelton runner designs were analysed with a consistent numerical setup, resulting in efficiencies of 91.65%, 91.95% and 92.30% for the baseline design and two design updates, respectively.
Nevertheless, the calculated efficiencies are partly subject to strong model assumptions and simplifications, such as the non-consideration of the turbine housing, the missing impact of the distributor, and the flow modelling. Therefore, it is essential to point out that the efficiencies given in this paper should not be confused with experimentally determined efficiencies.
In various publications, see e.g. [5,6,11], discrepancies between numerically simulated and experimentally measured efficiencies of Pelton turbines were reported. Whereby over-and underprediction have been shown, primarily depending on the numerical modelling, including mesh density, and physical modelling. Based on this knowledge, the reported efficiency enhancement of 0.65% points between the baseline design and the design update 2 is primarily regarded as an indicator for the actual improvement rather than an absolute value for the efficiency enhancement. Verification of the simulation results still has to be done based on a model test.