NUMERICAL INVESTIGATION OF THE BLADE PROFILE EFFECT ON THE AERODYNAMIC PERFORMANCE OF A VERTICAL-AXIS WIND TURBINE

A thorough understanding of the various parameters that affect the vertical axis wind turbine performance and lifespan will make design this model of turbine not difficult task. This article presents a study of the blade profile effect on the aerodynamic performance of a vertical axis wind turbineDarrieus H-rotor. A complete series of simulations based on complete unstable URANS calculations are performed for a three-blade Darrieus wind turbine. Fluent’s software is used for the numerical solution. The aerofoils of blade chosen for testing were NACA 0012, 0015, 0018, 0021 and -6 pitch angle is taken as the reference case for comparison. In order to quantify the influence of the profile of the blade on the performance, the characteristics of the flow field around the rotor of the different configurations are studied. For different values of the tip speed ratio, dynamic quantities, such as torque and rotor power, are presented and analyzed. Also in this investigation of the flow flied, two parameters are carefully investigated: mesh resolution and time step size. In the analysis, it appears that these parameters affect result accuracy. Finally, the numerical result shows that the variation of blade profile directly affected directly the power production.


INTRODUCTION
Electricity can be generated from the wind by wind turbines. Which are broadly classified as horizontal axis wind turbines (HAWT) and vertical axis wind turbines (VAWT). HAWTs are currently dominant for largescale energy generation. Recently, special efforts have been made to develop the vertical wind turbine (VAWTs) in the built environment [1].This focus is motivated by future demand for sustainable and decentralized energy supply in urban areas. VAWTs have several advantages over horizontal axis wind turbines (HAWTs) in this environment [2]. These turbines are extremely quiet and this is due to their special design where the blades do not create the usual conical noise that occurs with conventional horizontal axis wind turbines. These wind turbines also have the advantage, their independence from the wind direction. Wind conditions at low altitude are favorable for VAWT implementation. Recently, the research is focused on the Darrieus vertical-axis machine using an H-rotor in which the blades are straight and parallel to the axis of rotation [3-4-5].
Several parameters can affect VAWT performance and lifespan. It is therefore necessary to understand the effects of these different parameters on VAWT design. Both numerical and experimental investigations regarding the influence of different geometrics parameters on the performance of these VAWTs are needed to achieve optimal operating shapes and sizes, suitable for certain locations. Among these parameters there is the profile of the rotor blade. The aerodynamic designer seeks an airfoil shape that achieves the best aerodynamic performance while taking into account the tradeoffs with other disciplines and control. Traditionally, the wind turbine performance is definedin terms of the power extraction and expressed as the power coefficient (C P ).

389
The vertical axis wind turbine has an axis of rotation perpendicular to the oncoming airflow and the aerodynamics of VAWT has an inherent unsteady aerodynamic behavior due to the phenomenon of dynamic stall and the wake coming from the blades in the windward part and the axis. The complexity of the problem and the need to find new design approaches led the authors to focus on research using VAWT's Computational Fluid Dynamics (CFD) models where an approximation of the continuity equations and of momentum of the Navier-Stokes equation is solved. A CFD problem modeled with an appropriate turbulence model will yield simulation results that approximate the physical realities and results of the experiment.
Among these research works found in the literature, we can quote, Healy et al. [6] who is investigated the influence of aerofoil camber on VAWT performance and concluded that the closer the cambered aerofoil is the better its power output. Beri, Het al. [7] in 2011 also, conducted an investigation on the performance of VAWT with cambered aerofoils. He concluded that cambered aerofoil's can extract the maximum extraction of power from the wind. Qin et al. [8] used the URANS method in their numerical modeling of the flow around the rotor of an H-type straight-blade VAWT. They presented the stall phenomenon and the strong blade wake interactions. Okeoghene Eboibiet al. [9] studied the influence of the blade chord on the aerodynamics and performance of vertical axis wind turbines with experimental and computational fluid dynamics methods. He indicated the performance increases with increase in Reynolds numbers. Castelli M. et al. [10] performed numerical simulations using RANS model for VAWT with a three, four and five-bladed rotor. He has concluded a decrease of the torque during the revolution of a vertical-axis wind turbine (VAWT) by increasing the blade number. For the turbulence modeling, Wang and Tao [11] propose a campaign of simulations based on complete unstable RANS calculations. They compared the numerical results to the experimental data. They concluded that both turbulence models SST and K predicted with reasonable accuracy for low angles of attack were the flow is not detached. Ferreira et al. [12] notes in his work that VAWT presents an unstable aerodynamics due to the variation of the incidence angle with azimuthal position and relative velocity. He concluded that the dynamic stall phenomenon has a significant impact in turbine power. Nobile et al. [13] demonstrated that SST model predicted well as compared to other models in their CFD analysis to predict the dynamics stall at low TSR. A 2D study was done on a three bladed VAWT with NACA0018 profile using different RANS turbulence model.
In this study, a numerical modeling of the flow around a three-bladed wind turbine rotor Darrieus type H (VAWT) is performed. The approach URANS Unsteady Reynolds-Averaged Navier-Stokes with turbulence models k-ω SST and Spalart-Allmaras is used. The Fluent's software is used for the numerical simulation. This work concentrates in particular to study the parameter which affects the performance of these turbines, the airfoil of blade for to find the best configuration of these turbines. For this purpose, two dimensional simulations are carried out to estimate the aerodynamic performance of H-rotors with different airfoil shapes NACA 0012, 0015, 0018, 0021, 0024 because they are best for the vertical axis wind turbine due to their aerodynamic symmetry [2]. The simulations are performed for different tip speed ratio (TSR) values from 1.5 to 3. The ability of the moving mesh model of a computational fluid dynamics (CFD) resolver is used. Torque and power variations as a function of tip speed ratio are calculated and compared with experimental data from bibliography. To understand the physics of the interaction between viscous regions, non-viscous regions and vortices, this work has also focused on the analysis of flow in these regions.

AERODYNAMICS OF H-ROTOR
The Figure1 shows the Darrieus wind turbine. The H-rotor consists of straight blades with a bearing surface cross section profile. The strut and support arms serve as structural support.

-wR
When the turbine is subject to a velocity field perpendicular to its rotation axis, consequentially, the blade is under the influence of a force parallel to the rotation which generates a torque that rotates the turbine.

Figure 1.Three-bladed H-rotor VAWT
The direction and amplitude of the relative velocity change with the rotation of the blade. The azimuthally angle directly affects the values of the angle of attack, as shown Figure 3. The maximum and minimum effective velocities are reached when the blade is at the azimuthally angle of 0˚ and 180˚ according to Marco R. Castelli et al. [10]. Figure 2 shows the directions of the forces acting on the blade during its rotation. We have two main forces, the lift and the drag. The projection of these two forces gives the normal and tangential components.

Figure 2 .Forces and velocity distribution on Darrieus rotor
These forces act on the blade, the tangential force is parallel to the direction of rotation is the force and the normal force is perpendicular to the direction of rotation.
The angle of attack is varied depending on the azimuthally angle and the speed ratio (TSR), which is defined as the ratio between the speed of rotation of the blade (R) and the speed (U∞) of far stream not disturb.
where α= tan -1 [sin θ/(cosθ+ λ)] The incidence angle is related to the tip speed ratio, according to the figure 3. it is noted that, if the TSR increases the incidence angle decrease.
These values of incidence angle help to increase the efficiency of the turbine delaying dynamic separation because the dynamic stall angle increases. FT= L sin α− D cosα (4) where L is the lift force and D is the drag force FN= L cosα+D sin α where ρ is the density and U∞ is infinite velocity Journal of Thermal Engineering, Research Article, Vol. 6, No. 6, Special Issue 12, pp. 388-402, December, 2020

392
The averaged tangential force (F Ta ) The moment is calculated by this expression: The moment coefficient is calculated by this expression: The area (A) in 2D simulations is calculated by: The overall power (P) is defined by: P= Q · ω (12) where ω is the rotational velocity The power coefficient is defined by: The solidity is given by the equation follow:

MODEL GEOMETRY
The rotor studied here is similar to the rotor of Qing'on li et al. [8], this three-bladed wind turbine similar to the wind turbines used in the firm wind turbine. Figure 4 shows the Rotor with the azimuthally position, that was identified by the angular coordinate of the pressure centre of blade No. 1 at 0.25 chord for airfoil [10] [14].  Table-1 shows the geometrical basis parameters of the VAWT rotor. These parameters were taken accurately and digitized for a rotor. The blades are untwisted and are fitted with a constant pitch angle equal 6 degrees. This geometry of the digitized rotor was used in order to generate fluid volumes which constitute the various fields of calculation.

COMPUTATIONAL DOMAIN AND MESHING
In order to minimize the effects of the domain limits on the simulation results, several numerical tests have been done on the dimensions to determine the most suitable as indicated in the bibliography [5]. The selected geometrical parameter values of the domain of the solution 2D, is 60 R by 30 R. R is rotor radius of the VAWT. In order to allow a full development of the wake the outlet boundary conditions respectively 20R upwind and 38 R downwind. The figure 5 shows this 2D solution domain.

Figure 5. 2D computation domain
This 2D model includes the rotor with three rotating blades without support arm or shaft. The radius of the rotor is 1 m and the chord of NACA profile is 0.265m. Therefore the solidity of wind turbine is equal to 0.4. To take into account the rotor rotation, the flow field is divided into two distinct zones of mesh: an external fluid field which is the wind flow tunnel and a field that represents the internal fluid the turbine rotor. Thus, we construct the computational domain by grouping the mesh of the two areas.
A moving mesh technique was applied to the rotating circular area in order to predict the overall torque generated by the blades. For this purpose the correct choice of distance between the rotor and the input speed limit condition is very important, Hence the left side boundary of the domain is defined as velocity inlet [14].The boundary on the right of the domain was specified as outlet boundary condition , the static atmospheric pressure was imposed. On the sides of the domain the high and low , the condition of symmetry is used. at the inlet, the constant speed condition is used. the wind profile is constant and equal to 8 m / s with a turbulence of 5% of the free flow velocity. A no-slip boundary condition is specified at the blade , Both interface boundaries of the rotor domainand the outer domain were defined as interface boundary designed for sliding mesh application.

394
The figure 6 shows the mesh of the computational domain. The outer domain that represents the far field has a rectangular outer boundary. The computational domain is divided into nine parts (split with edges) for used the structured mesh at the far field of the rotor in the stationary domain and unstructured mesh around the rotor blade in rotated domain. In order to resolve the unsteady flow around the rotor blade, a boundary layer mesh is used around the blades. The first grid spacing on the airfoil was determined to make y+ less than 1, the growth ratio was limited to less than 1.12 to ensure numerical stability.

Figure 6. 2D Domain mesh
So that the results of simulation are not affected by the change of step of time and mesh size, a successive study was carried out on the choice of mesh. we have carries out a series of tests of independence of the time step, Three cases different were examined for a time step Δt equal 0.001, 0.002 and 0.003 which corresponds to 1 °, 1.5 ° and 2 ° azimuthal angle. The results obtained for these different sizes of time steps are shown in figure 7 where it can be seen that the difference in torque does not exceed 4% along a turbine revollution.

FLOW SOLVER AND TURBULENCE MODEL
To simulate the flow, the unsteady incompressible Reynolds Averaged Navier-Stokes equations are solved. In the URANS resolution method the time-averaged equations are solved. The K-ω shear stress transport (SST) model developed by Menter [15] is used to model the turbulence. For the pressure-velocity coupling a centered Simple algorithm is used. The unsteady terms are implicit second-order discredited. The computations have been obtained using Fluent a Computational Fluid Dynamic Package. This software uses the finite volume method. For the convection and diffusion terms, a third-order MUSCL scheme discretization is used. It should be noted that Rezaeiha et al [16] in his work, selected these same numeric parameters To calculate the time step, the following relation is used: (π/ ((180 * ω))). Therefore The azimuthal increment for the unsteady simulations of the flow auround the wind turbine rotor is equal to 1 °. To ensure that scaled residuals drop below,The number of iterations per time step used is 10. The values for the unsteady simulations are obtained after a computation time for 4 rotations of the rotor, these values were considered sufficient to obtain a convergence of the computations.

RESULTS AND DISCUSSION
The velocity field analysis allows to understand the performance of the turbine. The figure 8 shows the velocity field around the turbine rotor for a TSR equal 1.75. So, the visualization of the flow makes it possible to identify the wake downstream. this wake develops downstream of the rotor over a distance of more than 4 times the diameter of the rotor. the wake has regions in which the flow velocity is very low compared to the infinite upstream speed. the velocity field values are affected by the rotor operation . the wake is irregular and the interference between the blades is observed. These phenomena are less important for lower values of TSR. When the speed of the free flow is high, the wake propagates rapidly downstream. For to analyze the unsteady flow field around the profile blade, one presente in the figure 9 the velocity field at six different angles azimuth.The continual change of the angle of attack induces a circulation of the flow which changes continuously on the blades which generates the vortex formation.Two pairs of stall vortice appear. The first vortice is formed at the leading edge of the airfoil. A second vortice from the trailing edge rotates in the opposite direction. Together, they form the doublet of two vortices rotating in opposite directions, which goes downstream to meet the second blade. These vortex form a wake that disrupts in the downwind part . In these conditions, the downstream blades move in the turbulent wake on the downwind part. Although dynamic stalling increases the performance of the turbine, this also includes major drawbacks. It will cause an increase in the vibrations of the aeroelasticity and the fatigue of the blade. As can be seen in Figure 11, the dynamic stall will start earlier and the vortices are larger. It can be observed that in the same position of the blade profile during its rotation, the incidence angle of the flow separation changes gradually as a function of the TSRs.
The figure on the left shows the flow around the blade profile before the boundary layer separation. The figure on the right shows the flow separation and the formation of the vortex. Dynamic stall angle change with increasing TSR. This is the reason why the torque produced by the rotor increases upstream until reaching the maximum value, then decreases to minimum values. This decrease in torque is due to the flow separation on the profile wall, due to the reversal of flow in the boundary layer that produces vortices. The Darrieus wind turbine is particularly sensitive to dynamic stall because the change of the angle of incidence is great, especially at low speed reports.  Table 2 shows the different stall angle in function of the TSR. Dynamic stall angle depends on TSR. For lower TSR equal 0.5 the dynamic stall occurs at smaller angle of azimuth (38 0 ). The low tip speed ratios can also result in large separation of the flow around the airfoil of wind turbine. As the TSR increases, dynamic stall occurs at larger azimuthally angles.

Averaged Torque Coefficient
The average torque coefficient of the Darrieus wind turbine model was calculated for different TSRs. Figure 11-a shows the variation of the average torque coefficient as a function of the azimuthally angle for a TSR is equal to 1.75. The maximum value of the average torque coefficient is of the order of 0.25 for a value of TSR equal to 1.75. Figure 11-b shows also the maximum value of the average torque coefficient as a function of the azimuthally angles. This value is maximum for azimuth angles of 100 0 , 230 0 and 350 °. The averaged torque coefficient has nearly an identical Periodic variation for all the blades.

Torque Coefficient
Figure12 shows the variation of the instantaneous torque coefficient of a blade with the NACA 0021 and NACA 0012 profile as a function of The TSR. The maximum values of torque are generated in the upstream part of rotation of the turbine. The product torque decreases as the value of TSR increases. The angular position of the maximum instantaneous torque is almost identical for all TSR values in the range of 100 ° to 120 °. Note that the blade with the profile NACA 0012 delivers a stable instantaneous torque and lower than the blade with the profile NACA 0021, in conclusion, the blade with the profile NACA 0021 is recommended for the turbines operating in zones where the speed the wind is strong.

Power Coefficient
The figure 13 shows the CP variation versus TSR. The results show that the NACA 0021 profile blades produce higher power compared to the others profile NACA at the same tip speed ratio except for NACA 0024.The profileNACA 0021 delivers a maximum power coefficient estimated at 0.45 at tipspeed ratio equal 2 bat the NACA 0012 provides the minimum coefficient in the order of 0.22 at tip speed ratio equal 1.75. The conclusion is that at low TSR, The thicker profile blades perform aerodynamic performance better than thinner profile blade. The reason is that thicker profile has greater resistance to dynamic stall. But at hightTSR the performance decrease of the thicker profile and that resistance to the dynamic stall disappears.These results prove also, that our turbine only works better in the range of TSR between 0.5 and 3.

Figure 13.Power coefficient variation CP
In these conditions, the angles of attack reach critical values with values of low TSR as 0.5, which promotes aerodynamic stall too early. This dynamic stall is the main cause of reduction of lift [14]. On the other hand higher speed ratios as in our case gives incidence angle values too low which results in a loss of lift and a drag increase. We can so estimate that the NACA 0021 is more efficient for the turbine VAWT

Numerical Model Verification and Validation
The numerical results are compared with the experimental data of Qin et al. [4] that are available in the scientific literature.The figure 14 shows the comparison of the power coefficient (CP) variation versus TSR. Similar trends in the performance curve are obtained from both the experimental data and the numerical modeling 400 compared to the experimental results for different TSRs. These deviations in the results are due to several reasons, notably the simplifications adopted in the geometry of the simulation model to the 2D, the neglect of the effect of the arms and the shaft, the roughness of the surface of the blades, the phenomenon of numerical diffusion is also a source of error in the numerical simulation of complex flows. Indeed, a numerical diffusion is assigned to all the numerical resolution schemes of the principal equations, because this last one results from the truncation errors which are a consequence of the writing of the principal equations in a discrete form. Thus, in our case, the flux around the rotor is complex (in rotation), numerical diffusion is important because the fluid flows are not aligned on the mesh cells and this digital diffusion is inversely related to the resolution of the mesh.

CONCLUSION
Numerical analysis was performed in order to understand the effect of airfoil on the behavior of a straightbladed vertical-axis wind turbine. This work has been focused on the evaluation of the aerodynamic performance of a vertical axis wind turbine according to six types of blade profile. The method Unsteady Reynolds Averaged Navier Stokes is used for this work. The URANS model is is able to predict the generation of the wake and the rotor aerodynamic performance. It also shows an acceptable sensitivity to the grid refinement and the time step, which makes it suitable for simulation. The transient turbulence model SST used as the turbulence closure gives a better results and closest to experiments data for the bibliography.
The moving mesh method is used to handle the rotation motion of the blades. A very fine mesh has been adopted around the profile and in the wake zone to capture the phenomenon of flow separation and vortexes. In the simulations with varying grid refinement, these simulations presented different behaviors of the flow field this proves that the quality of the mesh influences the numerical results.
The numerical results are compared to the experimental data of the bibliography for validation. The power coefficient predicted by the numerical simulation presents the same variation as that of the experimental data. It is observed that the power which is absorbed from wind by wind turbine mainly depends on upstream region of azimuth angle of 0-180 0 the power coefficient produced by the wind turbine is sensitive to the dynamic stall phenomenon and the turbulence level around the blades. It has also been observed that the power coefficient results of the simulations are overestimated with respect to the experimental results. These differences in the results are due to the simplifications adopted in the geometry of the 2D simulation model, the roughness of the blade surface and the phenomenon of numerical diffusion. The analysis of the numerical results confirms the torques generated on the blades of rotor is highly dependent on the airfoil shape selected and the TSR. The thicker profile blade performed aerodynamic performance better than thinner profile blade. This paper shows the aerodynamics involved at different profiles for a 2-D rotor. A better improvement can be achieved in the future investigation of a 3-D case. The URANS method will take into consideration the 3-D nature of the vortices developed during the rotation of the rotor with arms and mast. This is necessary because the aerodynamic performance of the blade profile of the rotor in VAWTs can have a substantial impact on both the design and power generation of a wind turbine.