DESIGN AND HYDRODYNAMIC ANALYSIS OF HORIZONTAL-AXIS HYDROKINETIC TURBINES WITH THREE DIFFERENT HYDROFOILS BY CFD

The conversion of kinetic energy that comes from low-head water currents to electrical energy has gained importance in recent years due to its low environmental and social impact. Horizontal axis hydrokinetic turbines are one of the most used devices for the conversion of this type of energy [1], being an emerging technology more studies are required to improve the understanding and functioning of these devices. In this context, the hydrodynamic study to obtain the characteristic curves of the turbines are fundamental. This article presents the design and hydrodynamic analysis for three horizontal axis tri-blade hydrokinetic turbine rotors with commercial profiles (NACA 4412, EPPLER E817 and NRELS802). The Blade Element Momentum (BEM) was used to design three rotors. The DesignModeler, Meshing and CFX modules from the ANSYS® commercial package were used to discretize the control volumes and configure the numerical study. In addition, Grid Convergence Index (GCI) analysis was performed to evaluate the precision of the results. The computational fluid dynamics (CFD) was used to observe the behavior of the fluid by varying the speed of rotation of the turbines from 0.1 rad s-1 to 40 rad s-1, obtaining power coefficient of 0.390 to 0.435. For a maximum shaft power of 105W. In addition, it is evident that for the same conditions the rotor designed with the EPPLER E817 profile presents better performance than built with the NACA4412 and NREL S802.


INTRODUCTION
Environmental and energy problems worldwide have forced researchers to seek sustainable energy solutions. They have found renewable alternatives such as wind, sun and water to generate clean and environmental-friendly energy [2]- [5]. The energy available from water resources is the world's largest and cheapest [6]. Large hydroelectric have been criticized because they require dams and infrastructure. Civil works for the construction of these projects caused damage to the environment and ecosystems [7]. Accordingly, small hydroelectric plants have become one of the most economical and environmentally-friendly sources of energy generation in rural areas [8]- [9]. Hydrokinetic turbines are an emerging technology and they have gained strength and recognition because they provide a low-cost energy solution with low environmental impact [1]. This technology is still under development and it is in the research stage [10]. It is necessary to carry out more hydrodynamic studies of the operation of these devices to understand the phenomena that allow the conversion of hydraulic energy into mechanical [6]. The use of algorithms and mathematical models has facilitated the design of the hydrokinetic's rotors. Some authors have used semi analytical models to improve performance characteristics [11] or the Blade Element Momentum (BEM) theory to calculate the local forces on a turbine blade from being generated [12] [13]. The optimal use of water is very important for this type of turbines due to the low density of electricity generation. There are studies focusing on the improvement of the efficiency and operation of hydrokinetic turbines using different aerodynamic profiles and changes in curvature, the thickness [10] [11] [14][15] and yaw angle [16]. Several numerical and experimental investigations have been carried out to evaluate and improve the performance of hydrokinetic turbines [17]- [19]. Several authors have studied different hydrodynamic profiles [18], [20]- [22] and multi-elements hydrofoil [20] [21] . However, these studies evaluate flows, sizes and hydrodynamic profiles different from those reported in this research. For this reason, the power coefficients are different. The present study was focus on the hydrodynamic analysis of three rotors by CFD. Each rotor was previously designed with different profile (NACA 4412, NREL S802 and EPPLER E817). The objective was to analyze the fluid speed and pressure near the turbine and to construct the Cp vs TSR (Tip Speed Ratio) curve of each rotor. There are no reports available of power coefficient curves for the rotors whit the parameters presented in this study. The paper is organized as follows. Design of three rotors are presented. Description of the simulation method. Finally, the computational results are exposed and discussed.

Choice of aerodynamic profile
The design begins with the determination of the rotor radius taking into account the following considerations: test bench was the swimming pool of ITM which has a depth of 1.5 m, the gap between the tip of the blade and the surface or the bottom must have a distance of not less than 2R [25]. In this regard, if the depth is 1.5 m and a distance of 2R must be kept upwards, 2R downwards and the turbine measures 2R, the radius of the turbine must be 0.25 m.
The density of the fluid (ρ) and the kinematic viscosity (ν) at 20 °C is used. The flow velocity (U) of the test will have values between 0.8 ms -1 and 1.4 ms -1 . Reynolds number (Re) was calculated with the initial assumption of a chord (c) length for this type of rotors is approximately R/10. The Fluid velocity is considered as the relative velocity (w) between the chord of the turbine and the velocity of the flow (U), this relative speed varies between 5 ms -1 and 6 ms -1 with this data Reynolds are obtained between 200,000 and 300,000.
Previously, three families of profiles were found that the authors stood out for their performance in aquatic applications [11], [17]. Normally wind profiles have been used for the design of hydrokinetic turbines but these are aquatic applications. For this reason, there is an interest in the study of the families NACA 44XX, Eppler E8XX and NREL S8XX. The Reynolds values obtained previously are entered into Matlab using the XFOIL [26] program of the Massachusetts Institute of Technology (MIT). This program searched the NACA family for the tilt and profile that had the maximum lift-drag ratio, from the NACA 4410 profile to the NACA 4424 profile, varying the angle of attack (α) of each profile from -10 ° to 30 ° with steps of 0.001 °. The same process was carried out with the families of hydroprofiles E8XX and S8XX. The profiles with the best lift-drag ratio per family for a Reynolds of 200,000 were NACA 4412. For the E8XX family, a search was made between profiles E817 to E837, the profile with the best performance was E817. Finally, for the S8XX family, we searched from profile S801 to S835, finding that the profile with the best performance was S802. From each of these profiles the coefficient of lift was obtained (C L ) and the drag coefficient (C D ).
The parameters of the chosen profiles are reported and compared in the Table 1. This profile NACA 4412 has been widely used for wind applications. The Eppler E3XX family has been highlighted by some authors, this group of profiles has a good performance in aquatic applications and some have been widely used as the E387 [11], [27], Others like the E817 have been used for aquatic applications as propellers [28] but not for generation turbines. The National Renewable Energy Laboratory (NREL) created the NREL S802 profile. The S8XX family has been used for air and water applications [29], [30].
The Airfoiltools tool was used to obtain the "X" and "Y" coordinates [30], from the official page the data is downloaded into an Excel® file where they are processed so that the coordinates of the intrados are on one line, the coordinates of the extrados are left in other. These values are saved in a pad and then exported to the Design-Modeler as lines.

Rotors modelling
The coordinates in "X" and "Y" of an aerodynamic profile are saved in a notebook and then this file is exported as a curve to the DesingModeler module of the ANSYS program. By means of the standard tool, 10 exactly equal profiles were replicated on the Z axis with a gap between profiles given by the division between the radius of the rotor and the number of sections n of 0.025m. Using the value of the chord and twist of the BEM methodology, each profile is scaled and rotated. Once the ten scaled and rotated profiles are covered by performing the skin operation, in this way the solid 3D blade was obtained. Using the standard tool, two more blades were created and placed at 120 ° and 240 °. Then the hub was designed with the following dimensions: a 0.48R diameter cylinder with a total length of 0.4R and a semi-sphere on the top of the same diameter. The vanes and the bucket join together to finally create the rotor The rotor built with the NACA 4412 profile (from now on NACA 4412) is similar to the rotor built with the NREL S802 profile (from now on NREL S802) but the rotor built with the EPPLER E817 profile (from now on EPPLER E817) is wider. These models were used to obtain the control volume that was simulated. A third part of the complete rotor was simulated to save computational resources, as shown in

Simulation Setup
Fig shows the control volume created for the simulation.
The largest control volume contains a rotating flow (small control volume). The turbine is located into the rotating flow. By means of a Boolean subtraction operation, the solid body of the blade and the hub were eliminated. In this way, the fluid remains around the blade for simulation. Both control volumes are meshed. Proximity and curvature were used to densify the region surrounding the blade and the rotor because this zone is where the phenomenology generated by the torque is produced.
Boundary condition of the fluids are shown in Fig . The fluid inlet is in the green zone and the outlet in the red zone. Interfaces were used on the flat faces of the control volume and defined as symmetric. The inlet was configured as a speed inlet, while the outlet was made as a pressure outlet and was set as a free discharge with outlet at relative pressure 0 Pa. The symmetric face configuration was selected because the geometry presents third order rotational symmetry. In addition, this configuration helps to reduce computational costs. Other values necessary to configuration are shown in 2. Turbulence model was chosen because it is combination of a k-ω model and k-ε model. It used the k-ω model in the inner boundary layer and k-ε model in the outer region. This turbulence model allows to accurately capturing the phenomena close to the boundary layer that generate lift and drag forces. Furthermore, this model also adequately captures phenomena that occur far from geometry. In addition, this model is widely used in the literature for this kind of simulations [31]- [34] The value of inlet speed is very common in rivers. The rotation speed range is where the maximum efficiency point is normally found. Simulation was configured as transient with a time of 4s and adaptive time step that varies between 0.01s and 0.00001s. The external domain was considered station-ary and the fluid surrounding the rotor was configured as rotating, with a speed of rotation ranging from 0.1rad s -1 to 40rad s -1 . The frame change was configured as "Frozen rotor" and the interface model was of general connection. The periodic interfaces were configured as rotational periodic.
After concluding the configuration of the simulation, torque on the surface of the blade was selected as the output variable. 12 cores were used in parallel for the double precision solution. When the first result was obtained, the torque was parameterized to continue with the mesh independence study.

Mesh independence studies
Mesh independence studies guarantee the independence of results with respect to the number of elements and the characteristics of the mesh. Fig 8 NACA 4412 rotor full pressure profile.The mesh is composed of tetrads from 300.000 to 30´000.000 elements. These studies consisted in observing the behavior of torque when number of elements of the mesh is increased. Fig 4 shows the number of elements vs relative error. The mesh is independent when torque variation is less than 3% compared to the densest mesh (relative error). Thus, the chosen mesh does not increase the error of the results and guarantee economy of computing resources. A mesh of 3.79 million elements was chosen for NACA 4412 rotor since it meets the mesh independence parameters with the least amount of elements to save computational resources. For NREL S802 a mesh of 3.79 million elements was chosen and for EPPLER E817 that has a more complex shape, a mesh of 11.1 million elements was used. It should be noted that all meshes  With the selected meshes, the study of power coefficient against TSR was parameterized. The angular velocity was varied from 0.1 rad s -1 to 40 rad s -1 to obtain power characteristic curves of each rotor. In addition, planes were created in different sections of the blade to obtain velocity and pressure profiles at those points for each of the rotors. In order to ensure the precision and measure the uncertainty of the numerical results, a mesh uncertainty study was performed. Table 3 shows the details and main metrics of three meshes used for the E817 uncertainty study. It can be seen that the metrics of each one meet the acceptability ranges. Fig. 5 shows three points for each mesh described in Table 3. For the coarse mesh, the torque is 3.23084 Nm, for the medium mesh 3.18678 Nm and for the fine mesh 3.15496 Nm. This shows as the number of elements increases the torque decreases.
With the previous results, the GCI is calculated, using torque as the main variable. Since the value of the convergence condition (R) is 0.722, the convergence condition for torque is said to be monotonic. In addition, the GCI was calculated for mesh 1 (GCI12) and for mesh 2 (GCI23) with a safety factor of 1.25 obtaining values of 0.203 and 1.449 respectively. The results show that as the mesh is densified, the error decreases, which is why the finest mesh was chosen to continue with the next part of the investigation, which consists of changing the TSR.

Hydrodynamic analysis
As shown in Fig 6, NACA 4412 rotor located in the center of the image. It can be observed, around the turbine the behavior of the speed represented in contours. As expected by the physical principles that model the problem, the speed of the water above the turbine is greater than found downstream. This is caused by loss of kinetic energy of the water when it hits the blades. The velocity in posterior part of the hub is lower, in the region speeds lower than 0.4 m s -1 due to part of the fluid transfer kinetic energy to the blades. As seen in Fig , a low velocity region is also seen downstream of the turbine just as the theory indicates. There are four sections at 30% of the blade, 50%, 70% and 90%. In addition, each row belongs to a different TSR. It can be seen for TSR = 1 the behavior of the fluid in the lower part of the profile is irregular, this performance remains throughout the entire blade. In addition, it can be seen, in the images presented for this TSR, fluid direction lines go in through the upper part of the profile. This phenomenon contributes to the low coefficient of lift. Consequently, torque will be significantly affected.
As the TSR increases, the direction in which the fluid hits the profile improves. This make increasing the lift coefficient. Although in the sections near the hub irregularities continue to be generated in the direction of the fluid, upon reaching the middle of the vane this behavior is smoothed, thus reducing losses. The same analysis was performed for the NREL rotor S802 and E817. The differences were not obvious enough to be shown in images but later it will be shown how these small differences make the operation of the rotors different.   which is located in the center of the image. The pressure behavior represented in contours can be observed around the turbine. The reddish colors represent overpressures, while the green and yellow colors represent under pressures in relation to ambient pressure. It can see that the fluid has a positive pressure. As the fluid approaches the turbine, the pressure begins to increase. Once the fluid passes through the blades, the pressure drops suddenly reaching negative values. A moment after the fluid completely crosses the turbine, the pressure returns to normal. This phenomenon can be corroborated with the physics that models the problem since they have the same behavior.
shows the behavior of the pressure around different sections of the blade by means of contours for a TSR = 1. It can be seen that the highest pressures are concentrated in the part of the extrados in all sections of the blade. The above is because in this area the fluid hits. In addition, in the intrados of all sections low pressures are present. This is due to the angle with the fluid hits the blade. Due to this, the lift coefficient is low in almost the entire blade, this generates low torque performance in the turbine. This phenomenon is repeated with less intensity for TSR = 2 and also shows the same behavior for the other two rotors. On the surface, pressure contours are generated when the turbine has TSR = 4. It can be seen that the highest pressures occur at the leading edge and that the pressure differences between the intrados (bottom view) and the extrados (top view) are not significant. This phenomenon indicates that the TSR that maximizes the power coefficient is close to this value, for the reasons mentioned in the previous figure. Fig 11 shows the power coefficient results from the simulation of the NACA 4412, NREL S802 and EPPLER E817 rotors for different TSRs. In the image, it can see three curves in the form of a parabola corresponding to the numerical power coefficient of each rotors and a straight line corresponds to the maximum power coefficient obtained by means of the BEM methodology without taking into account the losses. It is observed the maximum numerical power coefficient is close to the maximum Cp without losses. The numerical results present a similar behavior to the results obtained through the BEM methodology, which is a validated methodology with experimental results and which offers high reliability. On the other hand, the maximum of the three functions is between TSR = 3 and TSR = 4, this coincides with the design TSR which is 3.5 for the NACA 4412, 3.6 for the NREL S802 and EPPLER 817. To be more precise, a trend line of type 5 polynomial was created. The equation was obtained to calculate the local maximum. For NACA 4412, EPPLER 817, and NREL S802 the TSR = 3.9 is the optimal value at which this rotor must operate to extract the maximum energy from the water. For TSRs close to 6 the NREL and EPPLER rotors rotate in the opposite direction due to the high relative speeds. This makes the Cp take negative values. Table 4 shows the power coefficient values and the power for each profile from the simulation. The power values obtained with the BEM theory are also presented. It can be seen the EPPLER profile is the one with the highest power coefficient and, in turn, the highest power.  The power values obtained with the BEM theory do not take into account the losses, which means that the results are above the numerical ones since the latter does consider some losses.

CONCLUSIONS
Three rotors were designed using the BEM methodology with different TSRs that maximize the power coefficient of each design, with Cp results greater than 0.442. The differences of chord and twists between rotor EP-PLER E817 and the other two are notorious, while the differences between the NACA 4412 and NREL S802 cannot be seen in the design, but in numerical form.
The rotor with better performance is the NREL S802 since it has a higher power and power coefficient. The high speeds must be produced in the lower part of the blade and the lower ones in the upper part, in order to generate the pressures to drive the rotor. In addition, a laminar flow will have higher power coefficients.
To obtain the proper performance of the rotor, the relative angle between the fluid and the blade must be optimal, thus high pressure regions appear in the intrados and low pressures in the extrados, since this phenomenon will generate the torque. Although the three rotors presented similar phenomenology, the differences can be seen in the power coefficient. The rotor with the highest power coefficient was the EP-PLER, followed by the NACA and finally the NREL. The torque of each one of the rotors according to the coefficient of power found is 117.0W for the EPPLER, 106.6W for the NACA and 105.6W for the NREL.