FLUID FLOW MODELLING THROUGH AN AXIAL-FLOW MICROHYDRO TURBINE

Original scientific paper In this study, a three-dimensional fluid field of axial-flow type microhydro named Agnew has been investigated. All numerical simulations were performed using OpenFOAM, an open source computational fluid dynamic (CFD) code, to investigate the rotor-stator interaction and losses occurring in the turbine. A sliding interface plane was used to pass the disturbance of stationary domain to rotary domain and vice versa. Several turbulence models were also examined together with wall function in order to take into account the turbulence in the flow. The obtained results show that a high resolution advection scheme, mixing plane to model the rotor-stator interaction together with a turbulence intensity of I = 6 % at the inlet, best matches with the experimental results.


Introduction
Through the last decades, application of the microhydro turbines has acquired a greater attention.Microhydro is a useful way to provide power of houses, workshops or villages that need an independent supply.An investigation on available potentials of the microhydro in Iran revealed that most of them are suitable to install the axial flow type turbines [1].Therefore, on the basis of a mutual corporation between University of Glasgow and the Iranian Research Organization for Science and Technology (IROST), a laboratory for testing the water turbines has been designed and equipped to examine performance of an axial-flow type microhydro turbine named Agnew [2].This turbine has the mechanism to adjust stagger angle of runner blades, which provides the maximum efficiency in a wide range of operating conditions [3].
Numerous studies have been carried out both experimentally and numerically to investigate the flow field in hydro turbines and their components.Jain et al. [4] modelled a Francis turbine and compared two sets of boundary conditions: (i) pressure inlet and pressure outlet, and (ii) mass flow inlet and pressure outlet.The obtained results prove that the second sets of boundary conditions are more appropriate for CFD analysis of the hydraulic turbine.
Vu et al. [5] predicted efficiency for a propeller turbine by means of the k−ε turbulence model at the steady state.A comparative CFD simulation between a full turbine simulation including six different blade geometries and another simulation including periodic conditions was also carried out.Both computations give relatively similar results, but some differences were found in prediction of the draft tube flow.One possible reason is the difference in the turbulence intensity level which is developed at the runner outlet, and leads to different results in the velocity profile below the hub.
Prasad et al. [6] applied the k−ω SST turbulence model for numerical simulation of an axial-flow hydro turbine.They computed efficiency of the turbine in some critical points and compared their results with the experimental values.They realized that the value of the efficiency which was computed from a full 3D analysis is in a close agreement with the experimental value at this regime.
Guénette et al. [7] predicted the hillchart of a bulb turbine through the k-ε turbulence modelling simulations.Their published results using RANS models contain useful predictions for runners opening around the best efficiency points.
Keck et al. [8] used a full stage simulation method to calculate the hill chart of a Francis turbine and compared it with the experimental results.They used TASK Flow code to solve three-dimensional equation of Reynolds Averaged Navier-Stokes (RANS), while the turbulence effects are modelled using standard k−ε model.
Choi et al. [9] validated CFD modelling for performance improvement of a 500 kW Francis turbine by the means of ANSYS CFX in order to analyse 3D turbulent flow at constant rotational speed.The average values of flow parameters like velocities and flow angles at the inlet and outlet of runner, guide vane and stay vane of turbine are computed to derive the flow characteristics.Their obtained results were found in a good agreement with the in site experiments, especially for the characteristic curve.
Shojaeefard et al. [10] showed that the swirl components of the inlet velocity vector are the most important performance parameters in the shape optimization of draft tubes.They found that the pressure recovery factor increases with the height and angle values over the design ranges.
In the present study, the axial micro hydro turbine named Agnew is simulated using the open source CFD package of OpenFOAM.The flow computation concentrates on a steady state simulation and comparison on the consequence of using different types of advection discritization as well as different turbulence models such as k−ε, SST, Reynolds Stress Model and Explicit Algebraic Reynolds Stress Model.In addition, the effects of surface roughness, rotor-stator interactions and turbulence intensity of an entering fluid flow on performance of the turbine are also investigated.The main results are integrated into a set of performance curves of the complete Agnew turbine.The flow characteristics of the entire turbine are assessed and compared with the experimental data at different measured operating points.

Numerical modelling 2.1 Geometry description
In the case of three-dimensional viscous flow analysis, a geometric model of runner and draft tube was developed separately and then assembled through proper interfaces to form a complete flow domain.An Agnew turbine is comprised of a 45° inclined axial microhydro turbine the blades of which can be adjusted at different angles.The 45° bend has stood in the inlet part of the turbine that leads to decrease diameter of the inlet pipe to diameter of turbine shroud from 254 to 155 mm.This reduction in the diameter causes increasing of kinetic energy and turbulence intensity at the inlet.Since the flow enters the runner axially, the 45° bending can be neglected and the inlet part should be considered far enough from the runner.

Grid generation
The computational domains are discretized by a combination of the clustering O-grids around the blades and H-grids topology in upstream and downstream domains [11].Using the O-grid around the blade yields a well-established boundary layer resolution and nearorthogonal elements on the blade.For the draft tube, the hexahedral meshes were also made.The criteria used to generate the mesh involved a minimal angle above 25 degrees, having a determinant greater than 0,6 and an aspect ratio below 200 [12].The four different grid densities were studied and are compared in Tab. 1.In these simulations, the rotational speed and the mass flow rate were calculated as 1000 rpm and 47,2 l/s, respectively.In order to evaluate sensibility of the results versus the grid size, it seems necessary to carry out a study on the grid independence.The reference values are the mechanical torque generated on all the runner blades and the pressure generated at the inlet.A Grid Convergence Index (GCI) method is used to estimate the discretization error [13].The GCI that is based on the RE method combines the often reported relative difference between two discrete solutions with a (r p −1) factor required in the denominator.The GCI takes a further step for converting the error estimation into an error or uncertainty band, which is again appropriate when one does not have a high degree of confidence in the error estimate.The GCI for a numerical solution of the fine grid and extrapolated values are defined as below: where, e  21 and  21  are the approximate relative error and grid refinement factor which are respectively calculated from: where,  should be greater than 1,3 [14].This value is found based on experience, rather than formal derivation.
is a critical variable to the conclusions being reported and ℎ is a representative cell, mesh or grid size.For example, regarding the three-dimensional calculations defined as below: In Eq. ( 2), p represents the accuracy order being defined as: where, which  32 =  3 −  2 ,  21 =  2 −  1 and   denote the solution on the k th grid.In Tab. 2, this calculation procedure is recorded for four resolution grids.The local order of accuracy p for the efficiency ranges from 2,43 to 4,9, with an average value of  ave = 3,2.The avearage values of p for the output power and the inlet pressure are proved to be 0,85 and 4,7, respectively.According to the GCI method, it is found that the third grid can be acceptable.All the values of  ext 32 are presented in Tab. 2 [13].Fig. 3 shows that the coarse grid allows obtaining the power values for the turbine components of each simulation, indicating the number of cells in each part.The effect of grid size on efficiency is shown in Fig .4 [15].The Y + of N 3 grid density varies in the range of 1 to 30, which is the acceptable range for an automatic wall function treatment in the boundary layer of SST turbulence model [16].

Grid generation
The interaction between the stationary and rotating domains was modeled by the means of three mathematical interface models: mixing plane, frozen rotor and sliding plane.All these models are common in this fact that they lower the time-dependent interaction between the stator and the rotor domains to a steady-state problem, but they are different in the way how the upstream disturbances are transported across the interface [8].
As shown in Fig. 5, at the mixing plane, the disturbances are averaged circumferentially over the interface, but the frozen rotor ignores this interaction between the rotor and the stator.The most realistic model is the sliding interface in which this transient interaction is correctly predicted.
Figure 5 Three models for rotor-stator interface [8] To simulate the previously described stage, a mixing plane model was utilized for the steady-state problems, and the result of the frozen rotor was taken as an initial estimation of the unsteady solution, where the latter, i.e. sliding plane, was implemented for the unsteady simulations.

Turbulence modelling
The best turbulence model is simply the one that matches with the experimental data in the best way [17].Because of producing the virtual diffusion, the standard kε model is not suggested for the problems which include intensive non-isotropic and non-equilibrium effects in the flow.In other words, the values of turbulent viscosity were predicted by this model a bit greater.Because of the swirling flow and the curving streamline, the k-ε becomes problematic near the wall regions and cannot model the boundary layer well.Therefore, this model has a greater difference by the experimental data compared to the other models as illustrated in Figs. 6 and 7. On the other hand, the results of SST are acceptable although not quite accurate.In the SST, separation, swirling flow and reverse gradient near the leading edge and trailing edge are modelling better than k−ω and k−ε.Therefore, the results are nearer to the experimental data in comparison with the k−ε.The EARSM allows an extension of the current k−ε turbulence models to capture the effects of physical phenomena such as secondary flows, flows with streamline curvature and system rotation [13].Therefore, the obtained results are the same by the two equation models.

Boundary conditions
The turbine runner was defined in a moving reference frame with given rotational speed, while the draft tube was considered as a stationary reference frame.Since the simulations are based on the concept of Multiple Frame Reference (MFR), the method of transferring the disturbances between rotating and non-rotating components of the machine should be determined.The interaction between the stationary and rotating domains was modelled by a mixing plane.The outlet of the runner domain and the inlet of the draft tube are taken into account as an interface of the mixing plane.
The measured mass flow rate has been specified at an inlet which is normal to the boundary surface.The averaged static pressure has been set equal to the atmospheric pressure at the outlet boundary, located at the end of the draft tube.
It is also necessary to estimate the turbulence intensity and the length scale at the inlets.Estimating the turbulence intensity in the inlet is often difficult.Because of the high-speed flows inside complex geometries such as those in the turbomachines, the value of intensity normally ranges from 5 % to 20 % [17].The turbulence intensity, often referred to as the turbulence level, is defined for a fully developed pipe flow as [11]: where  ′ is the root-mean-square of the turbulent velocity fluctuations and U is the mean velocity.  ℎ is the Reynolds number based on the pipe hydraulic diameter,  ℎ .The turbulence length scale, L, is a physical quantity describing the size of large energy containing eddies in a turbulent flow.In the pipe flows, this can be estimated from the hydraulic diameter.In a fully developed pipe flow, the turbulence length scale is considered 7 % of the hydraulic diameter.Hence [11]: From Eqs. ( 9) and ( 10), the intensity and the length scale were approximated about 6 % and 0,0154 m, respectively.As the bend is located before the turbine, the turbulent intensity was assumed 6 %.In order to investigate the effect of turbulence intensity, 3 values of that (i.e. 6 %, 10 % and 15 %) were studied.
Being coarse in the physical geometries, the walls are modelled with two smooth and coarse types.The surface roughness may be considered as an aspect of geometry, which leads to an increase in the turbulence near the wall, with a large influence on the boundary layer growth and loss at low and medium Reynolds numbers.If the roughness is known and the boundary layers are fully turbulent, its effect on skin friction can be reasonably predicted.The conventional wall functions are made sensitive to the effects of fine-grain surface roughness through introduction of a dimensionless roughness height.If the average roughness height is denoted by k, then the dimensionless roughness height k + can be defined in terms of the friction velocity u* [18], by the following equation [11]:

Solution method
In order to determine the efficiency and power curves and in order to analyse the losses in a draft tube, the steady state simulation is believed to provide acceptable results [15].Accuracy of the prediction is sensitive to the advection scheme.Upwind, high resolution and some blend factors to blend between the first and second order advection schemes are specified to calculate the advection terms in the discrete finite volume equations.Implementation of the advection schemes can be cast in the following form: where  up is the value at the upwind node, and  is the vector from the upwind node to the integration point.Particular choices for  and ∇ yield different schemes as described below.By the upwind setting, the advection terms are found of the first-order accuracy.This is equivalent to specifying a blend factor of 0. This scheme is very robust, but it will introduce diffusive discretization errors that tend to steep spatial gradients.A value of 1,0 uses the second order differencing for the advection terms, whereas values between 0,0 and 1,0 blend the first and second order differencing, with an increased accuracy and reduced robustness as β approaches to 1,0.At higher values, some overshoots and undershoots may appear, where excessive diffusivity can occur at lower values.In high resolution setting, the blend factor values vary from 0 to 1.In the flow regions with low variable gradients, the blend factor will be close to 1,0 for more accuracy.In areas where the gradients change sharply, the blend factor will be much closer to 0,0 in order to prevent the possible overshoots and undershoots and also to maintain its robustness [11].
All the numerical solutions for the previously described stage have been attained by OpenFOAM, which is an unstructured multiple element open source code of Technical Gazette 22, 6(2015), 1517-1526 Finite Volume Method (FVM).The code employs a fully implicit solution strategy and uses a multigrid accelerated factorization technique for solving the discrete system of the linearized equations [19].
Convergence is judged by examining the level of the residuals of the equations, as well as monitoring the values of kinetic energy and turbulence quantities.The normalized RMS and maximum residuals of the momentum, mass and energy equations are also monitored.The adopted convergence criteria for the RMS residuals of pressure and velocity have been computed 10-5 in all the simulations.In addition, the convergence is verified when the pressure at the inlet and the mass flow at the outlet become steady [17,21].

Results and discussions
As discussed before, the intensity of turbulence in the inlet fluid flow was determined.The greater values of intensity were also examined to investigate the sensitivity of solution per different values of it.Figs. 9 and 10 indicate that variations of the inlet turbulence intensity have small effects on performance characteristics of the turbine.Therefore, the calculated value of I = 6 % was taken for the next simulations.Figs. 12 and 13 depict the results for different advection schemes.The upwind scheme gives the most robust performance of the solver but still suffers from the numerical diffusion.So using this advection scheme is not recommended to obtain the final results.When β=1, the solution is more accurate but less robust.The high resolution scheme has a good flexibility and it can adapt itself in the flow regions with low and high variable gradients.According to Figs. 14 and 15 and based on the aforementioned discussion, the high resolution is suitable to calculate the advection terms in the discrete finite volume equations.It can be seen that as the discharge flow increases, the output power of the turbine also increases.The overall efficiency of the turbine increases by increasing the discharge and reaches the maximum value at the design discharge and then starts decreasing.The output power predicted by the CFD reveals a good consistence with the model testing results reported by the manufacturer.Some deviations in the overall efficiency were also noticed; however the trend of both curves was the same.
Poor performance of the draft tube causes drop of the flow and reduces the turbine efficiency.Efficiency of the draft tube which is obtained from the distribution of velocity and pressure at the inlet and outlet is given below [9].Fig. 14 shows that the maximum efficiency of the draft tube is 75 %, while the minimum efficiency of the draft tube is 70 % in the worst design of Kaplan turbine [16].Hence, design of the draft tube has some problem.Maximum cone angle of the draft tube has been mentioned 12° in the literature [22], while cone angle of the available draft tube is 13°.Area ratio of the inletoutlet and length of the draft tube are also the main parameters for designing the draft tube.
The pressure distribution on the runner is shown in Fig. 15.It can be seen that the flow is accelerated in the runner resulting in a pressure drop.On the suction side of the blades, the velocity maximum can be observed in the same position, where in the minimum pressure could be detected [12].indicates that a low velocity zone is formed near outlet of the draft tube.The velocity decreases and the pressure increases from the inlet to the outlet of the draft tube due to the increased cross-area of passing flow.This may contribute to increase the power generation in the turbine [17].
Angular velocity of the runner cone plays an important role in creating the separation and thereupon efficiency of the micro turbine.The effect of angular rotation is investigated in the present study.Early separation at end of the hub leads to reduce the effect of the draft tube upon pressure recovery.The separation of recirculation zone increases the average mean velocity, and hence the drops will become greater.
As depicted in figures 18 to 21, the streamlines leave the runner corn for ω=500 to 900 rpm earlier than ω=1200.As angular velocity of the runner cone increases, the swirl decreases.At ω=1200 rpm, the streamlines are straight and attached to the runner cone.The runner cone does not entrain the flow anymore, it decrease intensity of the swirl issued from the runner.The contours of the axial wall shear stress on the runner cone are also illustrated in Figs.18 to 21.For ω=500, 700 and 900 rpm, the wall shear stress has larger positive values near end of the cone indicating separation.Separation is obtained earlier for ω=500 rpm.This corroborates the fact that the streamlines leave the runner cone earlier for this angular velocity.The separation occurs later for a ω=1200 rpm.This means a better utilization of the cone, i.e. a better pressure recovery.The flow is attached to the runner cone up to end of it [24].Further experimental study on the draft tube of Agnew microhydro turbine has been presented in [25].Technical Gazette 22, 6(2015), 1517-1526

Conclusion
The experimental approach for evaluating the performance of a hydro turbine is costly as well as time consuming.Conversely, the CFD approach is faster and a large amount of results can be produced at almost no added cost.The CFD approach for prediction of the efficiency of a micro hydro turbine was presented and compared with the model results obtained from a manufacturer.A comparison between the computed and experimental efficiencies indicates that the best efficiency regime indicated by computation from both the numerical and experimental approaches is the same and the value of computed efficiency from a full-3D analysis is in a close agreement with the experimental value at this regime.The investigation shows that from accuracy and cost and time points of view, the SST model is one of the best alternatives for modelling the turbulent flows.The results show that the effects of surface roughness are negligible in this case.However, the current research shows that the scheme of high resolution for advection, mixing plane for a rotor-stator interface and I = 6 % for the turbulence intensity at the inlet leads to better results.Accuracy of the computed results can also be improved by making finer grids.The difference between the efficiencies computed from both numerical approaches and experimental values may be ascribed to a numerical error (sensitive to mesh size, discretization, convergence level and round-off), a model error (sensitivity to turbulence model and interface models) and a systematic error (sensitivity to geometry, boundary conditions and steady/transient assumptions).

Figure 1
Figure 1 Compatibility of actual blade with modelled geometry The Agnew turbine has four runner blades of a tip-totip diameter of 150 mm and a 2,5 mm tip clearance over the runner blades.Stagger angle of the runner varies from 0 to 30°.Two stagger angles of 30° and 15° are considered as quite open and semi open operating conditions, respectively.The draft tube is the pipe whose area increases from the inlet toward the outlet to reduce the velocity of discharge water and minimize the kinetic energy loss at the turbine outlet.Figs. 1 and 2 compare

Figure 2
Figure 2 Overview of draft tube geometry

Figure 3
Figure 3 Dependence of runner blade torque on grid size

Figure 6
Figure 6 Sensitivity of power for different turbulence models

Figure 7
Figure 7 Sensitivity of efficiency for different turbulence modelsThe SSG model is a new six equation Reynolds stress model and predicts much more accurately than the k−ε and SST (nevertheless more costly).This model is based on the transport equations for all components of the Reynolds stress tensor and the dissipation rate[4,18].Turbulent viscosity term is not in the RSM equations.This makes the solution tend to instability.Existence of the turbulent viscosity in the RANS equations acts as a damper of disturbance and helps to converge the solution.Therefore, with eliminating this term in the RSM model, the numerical instability grows unintentionally and the convergency becomes rather difficult.Most of the turbulent flows demonstrate a non-isotropic nature.In this simulation, the values of principal Reynolds stress are close together.Therefore, the results of the two equation models and RSM are much close.Fig.8shows the values of principal stress.According to the results of Figs. 6, 7 and 8, the SST model can lead to the best accuracy, cost and time.

Figure 8
Figure 8 Value of principal stress

Figure 9
Figure 9 Sensitivity of efficiency to intensity turbulence at inlet

Figure 10
Figure 10 Sensitivity of torque to intensity turbulence at inlet Fig. 11 illustrates the effect of wall roughness on power generation and efficiency.It can be concluded that the effect of variations of surface roughness on skin friction is negligible.Hence, the walls were assumed smooth.Figs.12 and 13 depict the results for different advection schemes.The upwind scheme gives the most

Figure 11
Figure 11 Sensitivity of power to wall surface roughness

Figure 12
Figure 12 Sensitivity of torque to advection scheme

Figure 13
Figure 13 Sensitivity of efficiency to advection scheme

Figure 14
Figure 14 Efficiency of draft tube

Figure 15
Figure 15 Pressure distribution on runner

Figure 16 Figure 17
Figure 16 Pressure distribution on runner

Figure 18 Figure 19
Figure 18 Velocity streamlines and wall shear contours for N=500 rpm of open blades working with 3 pumps

Figure 20 Figure 21
Figure 20 Velocity streamlines and wall shear contours for N=900 rpm of open blades working with 3 pumps

Table 1
Four different grid densities