Flow Simulations in a Small Bulb Turbine Using Two-Equation Turbulence Models

Numerical flow simulations in a small bulb turbine by the use of Ansys CFX software were performed in this paper. Simulations were performed for three different RANS-based models: k   , k   and SST. For each of these models, seven different operating regimes were considered. In order to reduce computational effort, steady state simulations were performed. In all cases, good numerical stability and convergence of solution were obtained. Based on the obtained results, performance curves for each of selected turbulence models in different operating regimes are formed. In turbine’s optimal operating regime, velocity profiles in a selected cross section of the draft tube were calculated, as well as the static pressure distribution on runner blades. Comparison of obtained performance curves was performed. Analysis of the velocity profiles and distribution of static pressure are given for each of the selected turbulence model.


INTRODUCTION
Fluid flow in hydraulic turbines is very complex, due to the geometry of the flow field.In order to understand the flow phenomena that are present in turbomachinery, one can perform a series of experiments that are quite time consuming [1], or use CFD.Reynolds averaged Navier-Stokes (RANS) equations in combination with Boussinesq hypothesis are mainly used for industrial simulations of this kind of flows.Performances of hydraulic turbines can be estimated according to numerical simulations based on Boussinesq hypotesis [2], which assumes isotropic eddy viscosity.Modeling of eddy viscosity can be done by algebraic, oneequation or two-equation turbulence models [3].
Standard k   model, proposed by Launder and Spalding in 1974, is a semi-empiric two-equation turbulence model [4].Empiric equation for dissipation in the k   model describes dissipation of the smallest turbulence scales in the homogeneous turbulence.In this case, eddy viscosity can be regarded as isotropic.These flow conditions are rarely fulfilled in real flows.Because of the mentioned assumption of isotropic eddy viscosity, k   model achieved good results in homogeneous turbulence region.In the wall region, where eddy viscosity anisotropy is increased, wall functions are introduced in the k   model [5].In this manner, the influence of the near wall anisotropy is decreased by introducing the assumption of the near wall velocity distributions.Standard k   model in combination with wall functions was used for numerical simulations in [6], [7] and [8].In order to avoid the use of wall functions in the wall region and to allow integration of equations in complete flow domain, Wilcox introduced k   turbulence model in 1998.[9].Numerical simulation of flows using this model is performed in [10].Introducing variables  and  , which determine scales of turbulence, allows twoequation models to account the history effects like convection and diffusion of turbulent energy [10].In order to obtain more accurate calculations in regions with positive pressure gradient, Menter proposed SST turbulence model [11].SST turbulence model is a combination of good properties of k   model in the free-stream region and k   model in homogenous turbulence regions and in the wall region.Numerical simulations of flows using SST turbulence model were performed in [12], [13] and [14].
The aim of this paper is to present the difference in results of these turbulent models, regarding the velocity profiles in vortex core region and the performance of the small bulb non-regulated turbine.

COMPUTATIONAL DOMAIN AND NUMERICAL MESH
All numerical simulations were performed for small bulb turbine displayed in Figure 1.Geometry of this turbine was defined in [15].Runner of the turbine has three, while the guide vanes have six fixed blades.Profiles of all blades have constant thickness with semicircular rounded leading and trailing edges.Optimal operating regime of the turbine is defined with the following quantities: discharge .There are four parts in the flow domain of the turbine: inlet pipe, guide vanes, runner and draft tube.Discretisation of the draft tube domain was performed using ICEM CFD and block structured hexahedral mesh was applied [16].
Runner and guide vanes were discretized using Turbo Grid and block structured mesh was generated [17].CFX-mesh was used for generation of unstructured tetrahedral mesh of the inlet pipe.In the wall region parameters of mesh were adjusted for the use of wall functions.Meshes of flow domains were connected between each other by the use of general grid interface (GGI).A stage model was applied for runner and guide vanes connection modeling.The connections inlet pipeguide vanes and runner-draft tube were modelled by the use of the frozen rotor model.Mesh of the turbine is displayed in Figure 2. and Figure 3.

TURBULENCE MODELS GOVERNING EQUATIONS
Numerical simulations were performed with the use of RANS equations and Boussinesq hypothesis.By analogy with Newton hypothesis, Boussinesq hypothesis is given by the expression: In fact, equation ( 1) represents the decomposition of turbulent stress tensor into anisotropic and isotropic part.Newton hypotesis also can be expressed by equation: as the sum of anisotropic and isotropic tensor of viscous stresses.For homogeneous shear flows where 0, ij S i j   equation (1) leads to equality of normal turbulent stresses: that is physically incorect in most real flows [18].According to equations ( 1) and ( 2) tangential turbulent stresses and tangential viscous stresses are given by: Molecular viscosity is isotropic for most incompresible flows and it can be expressed as second order isotropic tensor: Contrary to it, eddy viscosity is anisotropic in almost all real flows, which represents the main disadvantage of Boussinesq hypothesis as it was shown in [19], [20].Expression (1) can be transformed as folows: Eddy viscosity in standard k   model is represented by [4]: Transport equations for turbulent kinetic energy and dissipation of the standard k   model are [8]: .
Obviously, equations for k   model contain five constants whose values are as follows: The k   model gives good results in free-stream and wake regions [3].For high Reynolds numbers turbulence in wake region is homogenous and eddy viscosity is regarded as isotropic.In the wall region, for lower values of Reynolds numbers, wall functions are introduced in the k   model.Logarithmic near-wall velocity distribution is given by: 1 ln , , .
Scalable wall functions are introduced in the software Ansys CFX with the objective to avoid the use of logarithmic velocity distribution in subdomains with positive pressure gradients, where zero values of friction velocity can be achieved.Scalable wall functions are given by: In a formulation of scalable wall functions velocity * u is an alternative to friction velocity u  .Other variables in expressions (8)  .However, with the introduction of scalable wall function, viscous sublayer calculations are excluded.For fully resolving boundary layer flow, automatic near-wall treatment is introduced in CFX.Automatic near-wall treatment ensures automatic smooth shift from wall functions to a low Reynolds number near wall formulation as the mesh is refined.As already stated, k   model was introduced in order to avoid the use of wall functions in the wall region.The eddy viscosity of k   model is given by [3]: In a boundary layer at low Reynolds numbers, k   model gives better results than k   model [21].According to expression ( 10) k   model use is problematic in free-stream regions where 0 Transport equations of turbulence kinetic energy and turbulence frequency of k   model are given by [3]: Constants of k   model are: 9 100  model, Menter [11] introduced Here the   , , The function 1 F is the blending function which is equal to unity in viscous sublayer and logarithmic region.In a fully developed turbulence region this function is null.Calculation of the constant 1  is enabled by the use of the following quantities:  is calculated by the use of: Eddy viscosity is given by the following expression: where certain quantities in the former equation are given with: 0.31, tanh(arg ), 500 arg max 2 ; , 0.09 and  is the magnitude of vorticity tensor.

BOUNDARY CONDITIONS AND CONVERGENCE
Mass flow boundary condition was appointed on the inlet section of the inlet pipe (Figure 1).According to the mass flow inlet boundary condition, a uniform distribution of velocity in the inlet section of the inlet pipe was assumed.Average static pressure boundary condition was appointed on the outlet section of draft tube.In this manner, good robustness in working of solver was achieved. and SST models.For all performed simulations first layer thickness of mesh in the wall region was defined according to 1 y   .First order upwind scheme was used for discretization of convective term in momentum equation and convective terms in transport equations for turbulence properties [23].
Steady state simulations are performed for seven different discharge rates and three selected turbulence models, so total of twenty one simulations were performed.Good convergence of solutions was achieved for all models and all operating regimes.Convergence was improved by adjusting the time scale factor value.For all performed simulations RMS was about 4 10  , so enough convergence of solutions was achieved [9], [16].Figure 4 displays RMS residuals of equations for SST model in the optimal operating regime.

NUMERRICALLY ESTIMATED PERFORMANCES OF TUBINE
Brutto head of the turbine was calculated according to difference between the average total pressure on the inlet and outlet sections of the turbine.Brutto power was calculated according to brutto head.Netto power of the turbine was calculated according to torque and speed of the runner.Torque on the runner was calculated by numerical integration of static pressure on the runner blades.Calculated dependence of head, power and efficiency on the discharge rate is displayed on Figure 5. Performance curves on Figure 5 are obtained by interpolation of values for seven different operating points, which are defined by seven different discharge rates.Figure 5 shows that the head is variable.That is because the turbine is unregulated.In order to maintain the constant head it is necessary to introduce system of regulation of the turbine.It is estimated, from performed numerical simulations, that in the optimal operating regime, the efficiency of turbine runner and the turbine casing is about 0.86 and 0.93 consequently.These two efficiencies give the value of turbine efficiency of 0.8.The main reason for low turbine efficiency is high speed of turbine runner, giving relatively high hydraulic friction losses.
According to Figure 5 it can be seen that optimal operating regime of turbine is defined by the head about 3 m and the power about 8.4 kW.Differences in head, power and efficiency for different turbulence models, according to Figure 5, are very small.That is because averaging of variables decreases visual difference of velocity and pressure field for different turbulent models.For industrial flow simulations, where it is mainly needed an estimation of turbomachine performance, making a choice between two equation turbulence model is not so significant.Selecting an appropriate turbulence model is significant when it is necessary to approximate velocity and pressure fields more exactly.As it is impossible to estimate differences in velocity and pressure fields for two different turbulence models by analyzing the turbine performances, further comparisons of selected two equation turbulence models were performed by analyzing velocity and pressure fields in the turbine.

VELOCITY FIELD IN DRAFT TUBE
It is difficult to visualize and compare complete velocity and pressure fields of a turbine.For this reason the comparison of flow fields can be performed in certain regions of the turbine's flow domain.In a draft tube behind the runner a swirl flow is formed that is characterized by a very complex flow structure.In turbines, the goal is to have swirl flow with the smallest possible swirl intensity, and the weakest swirl is achieved in the optimal operating regime.In performed analysis a comparison of velocity field is performed in a cross section behind the runner of the turbine, along two straight orthogonal lines, axis X and Y, as it is shown in Figure 6.For velocities circumferentially averaging, ten concentric circles with centers at intersection of X and Y axes, were appointed behind the runner, as it is shown in Figure 6.It is known that two equation turbulence model gives poor results in predicting turbulent swirl flow [19].Consequently, it is interesting to examine the differences in velocity profiles which are obtained for two different equation turbulence models in selected region shown in Figure 6.Profiles of averaged circumferential velocity which are obtained by averaging of circumferential velocity along selected concentric circles, in optimal operating regime of the turbine are displayed in Figure 7.In velocity distributions that are shown in Figure 7 one can see small negative values of circumferential velocity close to the wall of draft tube.This is a consequence of a weak swirl behind the turbine's runner in the optimal operating regime, and bad projected geometry close to external diameter of turbine's runner.It is evident that runner's geometry in optimal operating regime does not fulfill projected flow conditions defined by the constant swirl ( rW const  ).However, these small deviations of circumferential velocity in optimal operating regime have no significant influence on turbine's performance.Outside of the optimal operating regime, for smaller discharge values, these deviations of circumferential velocity can be significantly greater [24].Profiles of averaged radial velocity along selected concentric circles in turbine's optimal operating regime are displayed in Figure 8. Figure 9 shows profiles of averaged axial velocity along selected concentric circles in optimal operating regime of the turbine.It is evident that locally averaging gives clearer differences between turbulent models than calculating performances of turbine.Greater differences in turbulent models are evident in the vortex core region.The reason is the presence of swirl flow [6], as well as the influence of flows with positive pressure gradient around the trailing edge of runner hat that is present downstream of the runner.Downstream, this influence produces inhomogeneous eddy structure.
Consequently, the biggest differences between k   and k   models are observed in this region.However, volume of vortex core region is relatively small compared to the volume of region of runner and guide vanes.Further downstream, this flow region is positioned around the central axis downstream, so there is no direct influence of the flow in this region on the flow around runner blades.These are some of the reasons why differences of flow fields in vortex core region have small influence on the resulting differences in turbine's performance simulations.A good mutually agreement between the results of turbulence models in homogeneous turbulence regions is observable on previous figures.The near-wall region wasn't analyzed by circumferential averaging, so in this region there are not differences between models in above figures.This analysis was performed by displaying the velocity profiles along X and Y axes, previously shown in Figure 6.Circumferential velocity profiles are shown in Figure 10 and Figure 11.According to the displayed profiles, it a good agreement of circumferential velocity can be seen in homogenous turbulence regions for all of the used turbulence models.Further, circumferential velocity changes sign in the wall region of draft tube.That is because swirl on the runner's outlet is weak in optimal operating regime, and because geometry of the runner was bad projected close to the outside diameter.Also, asymmetry of casing geometry of this turbine contributes to the asymmetry of the circumferential velocity profile.This influence is most evidently shown in Figure 12, which displays axial velocity distribution in the vertical section of the whole turbine.Velocity distribution on Figure 12 is obtained in optimal operating regime using SST turbulence model.Distribution of axial velocities along X and Y axes are displayed in Figure 13 and Figure 14.Difference in calculated value of axial velocity is the biggest in the vortex core region.Although all analyzed models are based on the assumption that turbulent viscosity is isotropic, there are differences in obtained results in the vortex core region where turbulent viscosity is anisotropic.Distribution of radial velocities along X and Y axes are displayed in Figure 15 and Figure 16.

PRESSURE AND VELOCITY DISTRIBUTION IN THE TURBINE RUNNER
Static pressure distribution on outlet contours of turbine runner is displayed in Figure 18.There are very small differences in pressure field between selected models.Flow regions in the area of the leading and trailing edges of runner blades have the lowest static pressure values.Wider low pressure regions are positioned from leading edges close to the outer radius of the runner.Along these regions the possibility for cavitation development is the greatest.Blade loading diagrams are good indicators of pressures differences on runner blades.Blade loading diagram along stream-wise direction from leading to trailing edge on half span from hub to shroud of the runner is displayed in Figure 19.Here one can see significant differences in pressure values only around leading and trailing edge of runner blades.It follows from energy equation, because there are very small differences in mean velocity values in the area of the runner blades for different turbulence models.Mean values of velocity in passages between runner blades primarily depend on runner geometry and continuity equation, so dependence of momentum and turbulence model is not significant.Consequently, according to Bernoulli's equation, the pressure distribution in passages between runner blades, also primarily depends on runner geometry and continuity equation.That is because in optimal operating regime flow in passages between runner blades approximately correspond to homogenous turbulence.In homogenous turbulence regions selected turbulence models give very similar results.Velocity distribution in the wall region is calculated using wall functions, so there are no significant differences in values of velocities in the wall region, because selected turbulence models use similar wall functions.Absolute velocity fields for different turbulence models along runner passages are displayed in Figure 20.Differences between velocity values for selected turbulent models are the biggest in wake regions in the area of leading and trailing edges.It can be seen that velocity values along the passages between blades are very similar values for all selected turbulence models There is are big differences in area-averaged absolute velocities from inlet to outlet of runner for different turbulence models.Numerically calculated reaction ratio for different turbulence models in the turbine's optimal operating regime is given in Table 1.Similarity in pressure distribution in the area of the runner blade profile for different turbulence models is one of the reasons for similarity in calculated values of turbine's performance.

CONCLUSIONS
This paper deals with numerical simulations of a flow in a bulb turbine using several different turbulence models.Calculated values of head, power as well as the efficiency of the bulb turbine, using selected turbulence models, did not have great differences among themselves, as a consequence of variables averaging.As it was impossible to estimate the differences in the use of turbulence models by analyzing the turbine's performances, a comparison of selected turbulence models were obtained by the analysis of turbine's pressure and velocity fields.
Performed analysis showed that in turbines there are flow regions with completely different turbulence structure, so one cannot use only one turbulence model to describe the whole flow field.For example, the difference in Reynolds number in turbine's inlet pipe and the draft tube is of the order of  .There is also the influence of turbines geometry and the runner rotation.The latter is the cause of the swirl flow formation characterized by the anisotropic turbulence, which cannot be taken into account in the adequate manner by the use of two equation turbulence models.This is the region of greatest differences in the velocity profiles regarding model selection, as it was shown in this paper.Although vortex core region has no major influence on the turbine performance, it occupies the imagination of many researchers.Flow modeling in this region is a big challenge for researchers.
In the regions where turbulence can be regarded as isotropic, small differences in velocity profiles for considered turbulence models were obtained.Small differences in mean velocity profiles in homogenous turbulence regions of runner passages have as a consequence small differences in pressure distribution on blades.This may be one of the reasons why prediction of lift coefficient for runner blades based on potential flow theory gives good results for optimal operating regime of turbines.In other words it can be said that there are small differences between homogenous turbulence and potential flow theory regarded the lift coefficient calculations for small angle of attack of the blade.
There is no universal turbulent model, as turbulence modeling is an inexact art.Most commercial codes offer a variety of turbulence models.It is the phantasy and the knowledge of the researcher who makes their use possible, for them to work correctly in order to resolve the imposed problem.

Figure 1 .Figure 2 .
Figure 1.Overall dimensions of the small bulb turbine.

.
In the low Reynolds number zones k   model gives better results than k   model.However, the k   model does not adequately include differences of eddy viscosity in the wall region and in free-stream regions.In free-stream regions k   model gives better results than the k   model.In order to consolidate good properties of k   and k 

Figure 4 .
Figure 4. RMS of SST model in the turbine's optimal operating regime.All walls were regarded as smooth, no slip walls.Scalable wall functions were used for k   turbulence model.Automatic near wall treatment was used for

Figure 6 .
Figure 6.Location of calculated velocity profiles, at a distance of 300 mm from the runner outlet.

Figure 7 .
Figure 7. Distributions of circumferentially averaged circumferential velocity in the turbine's optimal operating regime.

Figure 8 .
Figure 8. Distributions of circumferentially averaged radial velocity in the turbine's optimal operating regime.

Figure 9 .
Figure 9. Distribution of circumferentially averaged axial velocity in the turbine's optimal operating regime.

Figure 10 .
Figure 10.Circumferential velocity distribution along X axis in the turbine's optimal operating regime.

Figure 11 .
Figure 11.Circumferential velocity distribution along Y axis in the turbine's optimal operating regime.

Figure 12 .
Figure 12.SST model axial velocity distribution in the turbine's optimal operating regime.

Figure 13 .
Figure 13.Axial velocity distribution along X axis in the turbine's optimal operating regime.

Figure 14 .
Figure 14.Axial velocity distribution along Y axis in the turbine's optimal operating regime.

Figure 15 .
Figure 15.Radial velocity distribution along X axis in the turbine's optimal operating regime.

Figure 16 .
Figure 16.Y axis distribution of radial velocity in the turbine's optimal operating regime.Axial velocity distribution along Y axis close to wall of draft tube is displayed in Figure 17.It can be seen the difference in velocity profiles very close to wall between k   and k   models.The difference is the

Figure 17 .
Figure 17.Near wall axial velocity distribution along Y axis in the turbine's optimal operating regime.Scalable wall function ignores all points of mesh for which it states 11.06 y   , while automatic near wall treatment also includes these points.Consequently, first point of mesh with calculated values of variables is farther from the wall for k   model, than k   based models.Therefore, very close to walls velocity values obtained from k   model in Figure 17 are greater than

Figure 18 .
Figure 18.Static pressure distribution on the outlet runner contours in the turbine's optimal operating regime.

Figure 19 .
Figure 19.Blade loading chart of the runner on 50% span in turbine's optimal operating regime.

Figure 20 .
Figure 20.Absolute velocity distribution in runner passages on 50% span in the turbine's optimal operating regime.
The near-wall mesh resolution depends on the choice of wall functions.Standard wall function cannot be applied on arbitrarily fine meshes.Scalable wall functions overcome this drawback, because they exclude parts of mesh for which