Numerical investigation of 3D unsteady flow around a rotor of vertical axis wind turbine darrieus ty

Abstract

Keywords: Wind Turbine Darrieus; Aerodynamic; CFD; DDES

Introduction

Renewable energy sources are becoming increasingly interesting. Wind energy is attractive due to its inexhaustible potentials. In the idea of decentralizing part of the electricity produced from wind energy, one way is to use directly the wind energy that is in the immediate vicinity of the places of consumption, such as the urban areas. Small vertical axis wind turbines (VAWT) provide solutions where the installation of a horizontal axis wind turbine (HAWT) of several hundred kilowatts would not be feasible in urban areas Malcolm. [1]; Jamieson. [2]. The integration of wind turbines into urban environments can take advantage of the wind deviation that is produced by buildings to concentrate the flow of air to the machines and potentially improve production Beller, C et al. [3]; Sunderland K. et al. [4]. In 1931, Darrieus revealed his conception of vertical axis wind turbine (VAWT), which is the area of research of today. These wind turbines, due to their design can operate with winds coming from all directions and are less subject to these disturbances than wind turbines with horizontal axis. They have been designed to respond to the constraints caused by the turbulence of the urban environment. They relatively produce lesser sound and can easily be integrated into the architecture of buildings Martinat G. et al. [5]. However, these turbines have a low power coefficient compared to the horizontal axis wind turbine. Therefore, many studies are done by many authors to study the parameters which are responsible for the low power coefficient of these turbines. Numerical simulation makes it possible to predict the flow around the rotors of these wind turbines to explore aerodynamic parameters. Thus, many researchers have investigated numerically aerodynamic performance of the wind turbine. Peng Bo et al. [6] presented a numerical simulation on the wind turbine with k-ω turbulence model at different attack angles. Ping Ma et al. [7] applied the RANS with the SST K-ω models k –ω. Mertens et al. (2003) [8] studied, by using CFD, the performance of H-rotors in skewed flow. They concluded that there is an increase in the power coefficient when H-rotor is in the skewed flow. Simão Ferreira C. et al. [9] applied the DDED with the turbulence models SA to simulate the complex flow field around the wind turbine. They concluded that the power generated by the wind turbine was in good agreement with the experimental data for the fixed flow conditions, but discrepancies in the predictions were observed for the blocked flow regime. To understand the physics of the unstable flow around the rotor of the wind turbine subjected to the effects of wall rotation and to calculate its aerodynamic properties of the Darrieus wind turbine, numerical simulation in 3D constitutes a good alternative. The aim of this article is to describe and analyze the unsteady flow predicted numerically considering the effects of arms like blade-arms interference, blade-wake interactions around the rotor type H and the effect of tip vortices. In order to explore the unsteady and

turbulent flow around the rotors of the Darrieus turbine and to analyze the aerodynamic parameters responsible for the turbine performance, a 3D numerical modeling using the DDES method is carried out. The Detached Eddy Simulation (DES) model has been proposed by (Shur et al. [10,11]. The DDES is a hybrid solution that merges the Reynolds averaged Navier–Stokes and Large Eddy Simulation methods (Spalart P. et al. [12]. The accuracy of the DDES approach has generally been higher than that of Reynolds mean Navier - Stokes methods, a number of works have demonstrated this, see (Mertens et al. [8]; M.H. Mohamed et al. [13]. The rotor model chosen for this study is similar to that used in experimental works of Qing’on li et al [14]. In section 3 of this study, a two-dimensional numerical model is used to justify choosing the modeling of the turbulence, of the time step, and the verification of the stability of the results. Therefore, the results presented concern the influence of the mesh and the influence of the turbulence model on the torque coefficients. In section 4, a configuration of a three-dimensional Darrieus wind turbine and a calculation domain are described, then an analysis of the unsteady flow and global coefficients of aerodynamic performance are presented. A comparison of these numerical results with experimental data from the bibliography of our wind turbine model is presented for validation.

Model Geometry

The H-rotor is constituted by three straight blades with airfoil. The H-rotor blades generate lift and drag forces by the rotation motion in the wind, and the torque is produced by the force in the tangential direction. The power is generated by the blades of the rotor. This prototype was used in the experimental work of Qing’on li et al (2015) at the laboratory of Mie University in Japan. In order to validate the numerical results of the investigation with experimental data from the bibliography of Qing’on li et al. [14], we have digitized the geometrical parameters of the rotor used in this investigation experimental. These geometric parameters are presented in the Table 1. Geometric parameters of the rotor Designation

table 01, the solidity of the rotor is equal to 0.40. The solidity is defined according to Ion Paraschivoiu. [15] by the relation: σ = N.C / D

Mathematical Formulation

The delayed detached eddy simulation (DDES) model is a modified version of the Spalart-Allmaras model and can be considered a more practical alternative to LES for predicting the flowaround airfoils. The DES approach combines an unsteady RANS version of the Spalart-Allmaras model with a filtered version of the same model to create two separate regions inside the flow domain: one that is LES-based and another that is close to the wall where the modeling is dominated by the RANS based approach. The LES region is normally associated with the high-Re core turbulent region where large turbulence scales play a dominant role. In this region, the DES model recovers the pure LES model based on a one-equation sub-grid model. Close to the wall, where viscous effects prevail, the standard RANS model is recovered. In this numerical simulation, two turbulence models Spalart-Allmaras and SST k-ω were used for the turbulence modeling. DDES with the Spalart-Allmaras Model The standard Spalart-Allmaras model uses the distanceto the closest wall as the definition for the length scale. This Length scale, d, is generally taken as the shortest distance at any point to the closest wall in RANS model. The DDES replaces everywhere with a new length scale dmax defined by Fluent theory. [16]: d max = d − f d max( 0.d − CDES ∆ max )

Where d is the distance from … to the wall, dmax is the maximum computational cell dimension. The empirical constant CDES has a value of 0.65. ∆max, is the maximum local grid spacing. fd is given by f d = 1 − tanh(( 20rd )3 )

For turbulence model SST k-ω in the freestream the constant is: β*= 0.0828 And where FDES is expressed as   Lt FDES = max  (1 − FSST ),1 ∆ C  DES max 

With FSST =0, F1,F2 where F1 and F2 are the blending functions of the SST K-ω model of Shur M. L. et Al. [11]

Two-Dimensional Model

To start the numerical investigation, two-dimensional simulations are used in a preliminary configuration of the numerical procedure in order to study the influence of the effect of the turbulence models, the effects of the mesh and the time step. A calculation of the approximate performance parameters is carried out, namely the torque and power coefficients. This calculation is performed for different values of rotor speed and an infinite wind speed of 8 m/s. Computational Domain and Mesh in 2D In this 2D numerical study, obviously, the effects of arms and the shaft were neglected. Figure 1 shows the computational domain with the boundary conditions of the zones. To take into account the rotor blade rotation, the computational domain is divided into two sub-domains for the moving mesh formulation; a rotating sub-domain (turbine) and a fixed sub-domain. The geometric parameters of the calculation domain are 60 R by 30 R. These dimensions are chosen to avoid the effects of far numeric fields as indicated in the bibliography, see see Marco Racitiet al. [18].

With k = 0.41 DDES with the SST K-ω model The dissipation term of the turbulent kinetic energy is modified for the DES turbulence model as described in Menter’s work. [17]

The turbulent length scale is the parameter that defines this RANS model:

Figure 3. Independence solution. Figure 2. Computational Domain Mesh.

In order to capture the flow field around the rotating blade, a very refined mesh was built in the rotating subdomain. The position of this sub-domain is respectively at

38. R to the outlet boundary. Theses parameters of position

allows complete development of the rotor wake. In the rotating domain, the geometrical parameters of the rotor digitized are identical to the model, where the 2D model includes the rotor with three rotating blades. The boundary conditions used for the computational domain are: on inlet, velocity inlet with a constant velocity equals 8 m/s, this choice is justified by the fact that this type of wind turbine operates at very low wind speeds; on outlet pressure outlet with atmospheric pressure has been set because the wind turbine operates under this pressure. For the boundaries of the other side’s (top and bottom), the boundary condition of symmetry was used. The outer boundary of the rotating sub-domain and the inner boundary of the flow field are defined as interfaces. Different values of the angular speed of the sub-domain have been imposed. Figure 2-a shows the grid of the computational domain. A structured mesh is used in the sub-domain fixed,this subdomain is subdivided into nine zones. The structured mesh is used to minimizes the file size of the mesh, it allows us to reduce the computation time and increasing the accuracy. Figure 2-b shows the grid of the rotating sub-domain. The mesh on both sides of the interfaces has approximately the same characteristics of the cell size. For the rotating sub-domain, great attention has been given where the computational grids have been constructed, the unstructured mesh has been chosen for the rotating sub-domain domain. Figure 2-c shows the grid near-wall of airfoil, a boundary layer mesh is applied around the wall of

the blades. The local meshes around the blades are refined for accurate and efficient resolution. To guarantee more flexibility of the mesh of the 3D geometry, the parallelepiped is divided into several blocks of the same parallelepiped shape. The fixed and rotating (rotor) regions are independent of each other by the two interfaces, allowing a dynamic simulation. The mesh of the two areas in the interface should generate as much as possible coherently, so that the efficiency and accuracy of the data exchange can be guaranteed between the two interfaces. To ensure sufficient modeling of the boundary layer, the height of the first mesh was less than 0.001 times chord of the profile in order to ensure that a y+ value is less than 1, which is the limit of the turbulence model which was chosen for the simulations. A growth rate of 1.1 was applied to the inflation of the mesh around the blade, which gives about 20 layers inside the boundary layer. For the remaining domain, a growth rate of 1.2 was used. A smoothing algorithm in the mesh software was used to reduce the angular skewness of the cells so that the maximum was less than 0.7. A series of tests of the time step independence of mesh quality are carried. The results obtained for different sizes of time steps and the refinement of mesh are shown in Figure 3. It can be seen that the difference in torque does not exceed 4% along with a turbine revolution. This value is very satisfactory to make calculations. Solver Setup In the numerical investigation to resolve the equations, the solver Fluent is used. The method of a finite volume that is implemented in the solver is used. The Simple algorithm was used to handle the pressure-velocity coupling that exists. A 2nd-order special interpolation scheme for pressure was used, along with a 2nd-order upwind discretization scheme for the momentum equation and modified

turbulent viscosity. Residuals convergence criterion for each physical time step has been set to 10-5. A time step was chosen small enough. The value of the chosen time step corresponds to the time for which the turbine changes its angular position with Δθ = 10. The number of inner iterations for each time step is set to 20. The values of 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. In order to analyze the effect of turbulence modeling on the numerical results, two models were used, the SST k-ω and Spalart-Allmaras models. According to the literature review (Menter F. R. et al. 2003), these two models are very suitable for capturing the dynamic stall and to modeling the rotating areas.

Results And Discussion

Figure 4-a and Figure 5-a show the absolute velocity fields around the rotor, for the model of turbulence SST k-ω and Spalart-Allmaras respectively. During the rotation of the blades, the angles of attack change continuously, this change induces a constantly evolving circulation on the blades. This circulation gives a wake, which is characterized by the formation of a strongly vortices zone, in the downstream part. The viscous friction losses in this zone are significant because of this recirculation phenomenon. The turbulent flow generated is particularly complex, with the separation of boundary layers and the presence of unsteady turbulent structures that can lead to aerodynamic interactions and aerodynamic noise.

Figure 6. Variation of the torque coefficient versus the azimuth angle.

The prediction of the two models of turbulence of these dynamic phenomena is different. The same flow conditions were taken to make the comparison. Figure 4-b and Figure 5-b show the velocity fields around the blade turning at an angular position of 45° for SST k-ω turbulence model and Spalart-Allmaras respectively. The model SST k-ω predicts a larger recirculation zone at the trailing and upper edges. This proves a resulting decrease in aerodynamic performance. Dynamic stall occurs at a smaller angle of attack. With the interference of the blades, the dynamic stall is more important. This means that the SST K-ω turbulence model is better than Spalart-Allmaras turbulence model, especially near the walls. Besides that, in order to represent the influence of turbulence modeling on rotor aerodynamic performance prediction, we present in the following sections the numerical results of the torque coefficient and the power coefficient. The equations of these coefficients are given by: The torque coefficient expression: CQ = Q/2 ρRAU∞2

A is the area calculated, R is the radius and U∞ is the infinite velocity. The power coefficient is defined by: CP = Q.ω/0.5ρAU3

Figure 6 shows the predicted torque coefficient for the wind turbine blade versus azimuthally angle for an entire revolution of the rotor. It is observed clearly that prediction of the torque coefficient is depended on the turbulence

Figure 7. Variation of the Power Coefficient versus the tip speed ratio.

models. When the first blade of rotor moving in the upstream region, the torque values increase until reaching its maximum at azimuthally angle equal 95°. The torque values start ​​ to decrease and become low in the downstream region. The reason for this tendency is due to the influence of the attack angle and the interference wake effect of the preceding blade. The model Spalart-Allmaras predicts a higher torque at the end of rotor rotation than that predicted by the model SST k-ω, this shows that this model predicted a higher power than that which is predicted by SST k-ω model. Figure 7 shows the power coefficient comparison of the rotor predicted by the turbulence models SST k-ω and Spalart-Allmaras. These results are compared between them and then with the experimental data from the bibliography for validation. Six tip speed ratio (TSR) from 0.5 to 3 are chosen. The tip speed ration is defined by: λ= ω.R/U∞

were ω is the angular velocity of the rotor, R is the rotor radius and U∞ is the infinite velocity of wind. The values of TSR choosen make it possible to analyze the characteristic bell shape curve that describes the behavior of the power coefficient of a Darrieus turbine. The results give an acceptable agreement obtained with the experimental data. The numerical results show that the model SP predicts a higher power coefficient compared to the model SST k-ω at the same tip speed ratio. The maximum power coefficient predicts by the model SST k-ω is closer to the experimental data, therefore the SST K-ω

model is better suited to numerical computation of aerodynamic rotor performance than a wide tip speed ration range. It is observed that the numerical results are overestimated compared to the experimental data for different TSRs. These differences 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 phenomenon of numerical diffusion is also a source of error in the numerical simulation of complex flows.

Three-Dimensional Model

The 3D numerical modeling considers the effects that are caused by the presence of arms and the effect of tip vortices on the flow field and on the performance of the Darriues turbine. The aim of this section is to describe and analyze the unsteady flow predicted numerically considering the effects of arms like blade-arms interference, bladewake interactions, and the tip vortices. The same wind turbine Darrieus configuration is considered, namely three straight VAWT blades. The profile NACA0021 and the rotor diameter D = 2m are the same for the two-dimensional

analyzes but with the rotor height H = 1.2 m for the third-dimensional. Computational Domain and Mesh The discretization of the domain of computation led to two distinct sub-domains for the moving mesh as the 2D numerical configuration. A rotating sub-domain that contains the rotor and a fixed sub-domain which determines the global computational domain. Figure 8 shows the domain of the 3D solution. The boundary conditions selected for this simulation are the up-stream speed and the downstream static pressure. The boundary condition without sliding is applied on the blades and arms of the rotor. The symmetry boundary condition was used for the two sides (top and bottom). The boundary conditions around the external surface of the rotating domain and the internal surface of the large domain have been defined as interface. A structured multi-block mesh was chosen for the large, fixed domain and an unstructured mesh for the rotating domain which is a moving mesh. Figure 9-a shows the mesh of the two sub-domains The detail of mesh rotor is presented in Figure.9-b. Solutions independent of the mesh have been found using this unstructured and structured mesh topology with approximately 11,000,000 meshes. The numerical configuration is the same configuration as that of 2D modeling since it was evaluated and analyzed with the time steps and the turbulence modeling. The computing power used for the simulation of the 3D modeling has been increased since the size of the model requires a more powerful calculator. 20 sub-iterations were necessary for the solution to converge at each time step. A computer with 12 processors and 2.10 GHz clock frequency is used. The simulations required a calculation time around 45 days.

Results And Discussion

3D numerical calculations made it possible to explore the unsteady flow around the rotor and the wake over a large field downstream. thus, the results give a visualization of the vortices emitted by the tip of the blade and the vortex of the arms. Figure 10 shows the velocity field around the turbine rotor for a TSR equal 1.75 presente the development of the downstream wake, generating vortices from the blade structure and arms. the wake develops downstream of the rotor over a distance of more times the diameter of the rotor.Two pairs of separated vortices appear. The first tourbillon is formed at the leading edge of the airfoil. The second vortex is formed from the trailing edge rotates in the opposite direction. Together, they form the doublet of two vortices rotating in opposite directions, which go downstream to meet the second blade. These vortices form a wake that disturbs the leeward part. Under these conditions, the downstream blades move in the turbulent wake of the downwind part. The velocity field values are affected by the rotor operation. The wake is irregular and the interference between the blades is observed. The wake has regions in which the flow velocity is very low compared to the infinite upstream speed. The downstream flow is divided into two parts: an inner part where the axial velocity is slowed by the vortices, and an outer part in which the velocity is accelerated by the vortices. The wake of the blades is very bulky therefore; the blades are a source of disturbances and high aero-dynamic instability as shown in the figure. The wake moves away from the rotor a great distance. In the Darrieus rotor, the arms are necessary to connect the blades to the tower. These structures have an influence

on the performance of the turbine. The effects of the arms like blade-arms interference and blade-wake interactions are often described using correction factors but are not yet fully understood. To understands and identify the characteristic of the vortex shape behind the rotating arms, the contours of velocity around the arms surfaces were visualized and are shown in Figure 11, which identifies the characteristics of the vortex shape behind the rotating arms. In addition, the dynamic of arms vortices is correctly captured with their inboard motion in the horizontal plane of the wake. On a rotor blade of finite length, the distribution of lift is not uniformand tends towards a zero value at the blade tips. So the lift gradient is very strong at the tip of the blade and this is manifested by the emission of a vortex called the tip vortices or marginal vortices. This vortex emission is in the direction of flow, therefore orthogonal to the blade. Figure 12 shows the marginal vortices that develop in the tips of the blades. Decomposition of the tip vortex occurs near the blade and in the distance the wake becomes very turbulent with irregular vortex structures. The investigation shows that this tip vortex is the vortex core has a high velocity gradient. It should be noted that the review of successive velocity fields over time shows a fluctuation in the position of the vortex nuclei. To analyze the variation of the torque during a rotor revolution, Figure 13 shows the variation of the instantaneous torque obtained for the three rotor blades as a function of the azimuth angle. It is noted that the maximum torque of the blades 1, 2 and 3 is reached respectively at the azimuth angle 95°, 215° and 335°. The maximum value of the

predicted coefficient of torque obtained is around 38 Nm at the TSR value which is equal to 1.75. This result is closer to the experimental data of Qing’on li et al. [14]. The variation of the instantaneous torque is identical for the three blades; it increases up to a maximum value which corresponds to the dynamic stall angle of the blade. After that the value of the torque decreases because the flow around the blade is completely separated. This torque even becomes negative at the end of the rotation cycle.

Figure 14. Power coefficient validation with the experimental data of Qing’on li et al. [14].

To analyze the aerodynamic performance behavior of the rotor, it is necessary to establish the power coefficient variation curve which describes the behavior of the Darriues wind turbine as a function of the tip speed ratio. Simulations are performed with a variation in the rotation speed of the rotor. Six rotational speeds corresponding to a variation of the tip speed ratio from 1 to 3 have been chosen. Figure 14 shows the power coefficient for the wind turbine blade versus tip speed ratio. The variation of the power coefficient has very good compatibility with the experimental data of Qing’on li et al. [14] for all the tip speed ratios. We note that the maximum value of the CFD coefficient is of the order 0.22 at TSR = 1.8 which is a value very closes to that the experiment data estimated at 0.185. This consistency is due to the taking into accounted the three-dimensional effects, such as the loss of the tip of the blade and the effects of the arms, like blade-arms interference and blade-wake interactions. There is always this acceptable difference which does not exceed 5% between the two curves. This difference is due to the numerical modeling and the numerical diffusion.

Conclusion

In this research a 3D numerical investigation of unsteady flow around a rotor of Darriues turbine type H has been performed. we first made simulations of the wind turbine in 2D to analyze the effect of different parameters of numerical modeling namely, the mesh, the time step and the turbulence modeling to discover their effects on the quality of the results, thus, a quantitative analysis of the performance and the flow around the rotor has been

presented. By comparing the velocity field, it is proven that the Spalart-Allmaras model can capture more vortexes after flow separations. The SST k-ω model predicts vortices that are less realistic than those of Spalart-Allmaras Model. From these results, the approach of the Detached Detached Eddy Simulation with the SST k-ω model can be considered as good alternative in the numerical prediction. Simulations with model 3D were performed to quantify the effects that is caused by the presence of arms and marginal tip vortices on the flow field and on the performance of the turbine. Numerical results of aerodynamic performance were compared to experimental data to valid the numerical model. The comparison of the power coefficients shows a good agreement for this 3D model. The knowledge acquired here can be used to continue the study of the wake of the rotor which contains vorticity, which is called a vortex emission. The dynamics of letting go of vortex play a crucial role in the functioning of a Darrieus machine. The vortex illustrated the level of circulation around the blades therefore gives an indication of the level of effort. The energy extracted by the rotor is based on this variation in circulation.

Nomenclature

The cross-section area of the wind turbine Chord length, m Diameter Rotor, m Power, Nm/s Torque, Nm turbine radius, m Torque coefficient, Power coefficient Renolds Number, infinite velocity of wind, m/s

Greek symbols σ Solidity, ω Angular velocity, rad⁄s θ Azimuthal angle, ° α angle of attack, ° λ Tip speed ratio ρ Air density, kg/m3

Abbreviations

AOF Angle Of Attack NACA National Advisory Committee For Aeronautics TSR Tip Speed Ratio URANS Unsteady Renolds Average Navier Stocks VAWT Vertical Axis Wind Turbine

Data Availability Statement

The authors confirm that the data that supports the findings of this study are available within the article. Raw data that support the finding of this study are available from the corresponding author, upon reasonable request.

Conflict Of Interest

The author declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Ethics

There are no ethical issues with the publication of this manuscript. [11]

References

1. The article references can be accessed from the .pdf file.

Share and Cite