- Email: [email protected]

Contents lists available at ScienceDirect

Energy Conversion and Management journal homepage: www.elsevier.com/locate/enconman

Three-dimensional Improved Delayed Detached Eddy Simulation of a two-bladed vertical axis wind turbine Hang Lei a, Dai Zhou a,b,c,⇑, Yan Bao a, Ye Li a, Zhaolong Han a,d a

School of Naval Architecture, Ocean and Civil Engineering, Shanghai Jiao Tong University, Shanghai, People’s Republic of China Collaborative Innovation Center for Advanced Ship and Deep-Sea Exploration (CISSE), Shanghai, People’s Republic of China c State Key Laboratory of Ocean Engineering, Shanghai Jiao Tong University, Shanghai, People’s Republic of China d Cullen College of Engineering, University of Houston, Houston, USA b

a r t i c l e

i n f o

Article history: Received 11 August 2016 Received in revised form 26 November 2016 Accepted 30 November 2016 Available online 18 December 2016 Keywords: Computational fluid dynamics H-VAWTs Polyhedral mesh 3D IDDES Aerodynamics

a b s t r a c t The aerodynamic performance of a two-bladed vertical axis wind turbine is investigated using the turbulence model of the Improved Delayed Detached Eddy Simulation and the polyhedral mesh. The sliding mesh technique is used to simulate the rotation of the rotor. Meanwhile, the results obtained by the shear stress transport k-x model are presented as contrast. Then, the simulated power coefficients at different tip speed ratios and the wake velocity are validated by comparison with the experimental data from available literature. It is shown that the power coefficients and wake velocity predicted by the Improved Delayed Detached Eddy Simulation are closer to the experimental data than those by the shear stress transport k-x model. The pressure distributions predicted by the two turbulence models show different degrees of discrepancies in different scales of flow separation. By comparing the vorticity magnitude graphs, the Improved Delayed Detached Eddy Simulation is found to be able to capture more exquisite vortices after the flow separations. Limited by its inherent ability, the shear stress transport k-x model predicts vortices that are less realistic than those of Improved Delayed Detached Eddy Simulation. Hence, it may cause some errors in predicting the pressure distributions, especially when the blades suffer dynamic stall. It is demonstrated that the Improved Delayed Detached Eddy Simulation is regarded as a reliable model to analyze the aerodynamic performance of vertical axis wine turbines. Ó 2016 Elsevier Ltd. All rights reserved.

1. Introduction Wind power is regarded as one of the most promising energy resources due to the cleanliness without air pollution and easy accessibility [1]. As a typical device that converts the wind energy into electricity, wind turbines have attracted a lot of interest during the past years [2]. Wind turbines can be classified into two main types based on the axis of rotation: horizontal axis wind turbines (HAWTs) and vertical axis wind turbines (VAWTs) [3]. Because of a large number of investments, the technology for HAWTs has rapidly developed in the past decades [4]. VAWTs have been disregarded for a long time due to limitations on the previous technology, however, recently VAWTs have received a growing amount of attention due to several unique advantages [5].

⇑ Corresponding author at: School of Naval Architecture, Ocean and Civil Engineering, Shanghai Jiao Tong University, Shanghai, People’s Republic of China. E-mail addresses: [email protected] (H. Lei), [email protected] (D. Zhou), [email protected] (Z. Han). http://dx.doi.org/10.1016/j.enconman.2016.11.067 0196-8904/Ó 2016 Elsevier Ltd. All rights reserved.

Modern VAWTs have been developed into a variety of types, which mainly include the Darrieus type, the Savonius type, and the H-Rotor type [6]. The H-rotor vertical wind turbines (HVAWTs), which can be regarded as a variation of the Darrieus type wind turbine were developed in the United Kingdom during the 1970s, and belong to the family of lift force driven wind turbines [7]. In addition, the aerodynamic characteristics of the H-VAWTs are more complex than other types, which will be a meaningful work to explore. With the popularity the Computational Fluid Dynamics (CFD) method, many complex aerodynamic problems have been resolved by this favorable tool [10]. Therein, Bianchini et al. [11] studied the virtual incidence effect on the airfoils of Darrieus wind turbines. Butbul et al. [12] studied the impact of inertial forces on the aerodynamics of VAWTs. Shahizare et al. [13] investigated the output power of VAWTs in urban area by considering the effects of the Omni-direction-guide-vane angle. Wang et al. [14] simulated the aerodynamic performance of a novel vertical axis wind turbine with adaptive blades. All the publications mentioned above

236

H. Lei et al. / Energy Conversion and Management 133 (2017) 235–248

Nomenclature c Cp D H k LIDDES LLES LRANS Q Sij t uj V1 Vt W xj y+

chord length power coefficient rotor diameter blade span turbulent kinetic energy IDDES length scale LES length scale RANS length scale rotor torque mean strain rate time velocity components free stream wind speed line velocity of blades velocity of the wind relative to the rotating blade coordinate variables dimensionless wall distance

Greek symbols a angle of attack with pitch angle a0 angle of attack without pitch angle b pitch angle b⁄ a constant in SST k-x D grid scale h azimuth angle k tip speed ratio

employed the CFD method. In addition, two-dimensional (2D) aerodynamic models were widely used in pneumatic analysis of VAWTs. Danao at al. [15] used the 2D Reynolds-averaged NavierStokes equations (RANS)-based CFD method to study the aerodynamic performance of a VAWT in unsteady wind. Ismail and Vijayaraghavan [3] utilized a 2D numerical method to investigate the NACA-0015 airfoil profile with optimized shape which was installed with an inward dimple and a Gurney flap. Buchner et al. [9] performed a 2D unsteady RANS model to analyze the dynamic stall of a two-bladed H-VAWT, on which a wind tunnel test was conducted so as to verify the simulation results. Gharali and Johnson [8] and Almohammadi et al. [16] studied the 2D dynamic stall of a single blade and a full straight blade respectively, and Gharali declared that a three-dimensional (3D) model was a better choice to investigate the dynamic stall phenomena. A 3-D unsteady RANS method was utilized by Chowdhury et al. [17] to compare the performances of a VAWT in the upright and tilted configurations by using the SST k-x model. It was shown that the torque produced by the tilted condition was higher than the torque produced by the straight one. Siddiqui et al. [18] engaged five different simplified geometric models, and also found that the 3D CFD simulation results were closer to the experimental data. Lee and Lim [19] studied a variety of parameters that affect the performance of VAWTs, including rotor diameter, chord length, pitch angle, and helical angle, and those results provided some guidance for the design of VAWTs. On the other hand, the Large Eddy Simulation (LES) is also a powerful turbulence model that has attracted more attention in recent VAWTs’ studies. The 3D LES model was used by Li et al. [20] to investigate and validate the aerodynamics of H-VAWT with a short segment of airfoil blades. Elkhoury et al. [21] attempted to understand the performance of a VAWT under a variable-pitch by using the LES, showing a good agreement with the experiments. Although LES is considered as a high-fidelity method, it requires a huge amount of computational cost. In order to obtain accurate

l lt q sij

Oij

x

molecular viscosity turbulent viscosity air density tensor of stress vorticity rotor angular velocity

Abbreviations 2D Two-dimensional 3D Three-dimensional CFD Computational Fluid Dynamics DDES Delayed Detached Eddy Simulation DES Detached Eddy Simulation HAWT Horizontal Axis Wind Turbine H-VAWT H-Rotor Vertical Axis Wind Turbine IDDES Improved Delayed Detached Eddy Simulation LES Large Eddy Simulation NREL National Renewable Energy Laboratory N-S Navier-Stokes RANS Reynolds-averaged Navier-Stokes SST Shear Stress Transport TSR Tip Speed Ratio TKE Turbulent Kinetic Energy

results and save the computational cost, a method combining the RANS model with the LES model was suggested, with one of them called the Detached-eddy simulation (DES) [22]. In the DES method, the small eddies near the boundary were modeled using a RANS formula and the large eddies in the far field were simulated by a LES type model [23]. The DES method has already been used to simulate the aerodynamics of HWATs successfully. For instance, Li et al. [24] studied the National Renewable Laboratory (NREL) phase VI wind turbine using the DES turbulence model. The results from their CFD simulations showed that the DES model had the ability to predict accurate fluctuations of the frequencies presented in the experiment. Johansen et al. [25] considered only a single blade of the NREL phase VI, and the parked blade case and pitching blade case were studied. They found the DES model still had almost no obvious improvement in predicting the global blade characteristics. Therefore, there exists some controversy in using the DES model to predict the aerodynamics of HAWTs for higher efficiency and better accuracy. In light of the controversy in HAWTs and the lack of use in VAWTs [26], whether the DES model can predict the performance of VAWTs accurately or not is an interesting topic that promotes the current study. Further, the Improved Delayed Detached Eddy Simulation (IDDES), established by Shur et al. [27], is a development of the classical DES model. It combines the SST k-x based model and the Delayed Detached Eddy simulation (DDES) [28], modifying the inaccurate treatment in the grid-induced separation for the DES. Therefore, the IDDES aims at wall modeling by the LES, and can resolve turbulence activity close to the wall more effectively when compared to the pure DES model [27]. In other fields, Hu et al. [29] studied the flow structures behind a backward-facing step using the IDDES model, and Xiang and Luo [30] used the IDDES model to simulate the massive separation around triple cylinders with good accuracy. For all the reasons mentioned above, the main purpose of this paper is to study the aerodynamics of an H-VAWT with the IDDES

H. Lei et al. / Energy Conversion and Management 133 (2017) 235–248

237

model by using the commercial CFD solver as it has been widely used in simulating the aerodynamics of wind turbines [26]. The polyhedral mesh and sliding mesh are adopted to generate the mesh topology. In order to demonstrate the advantages of the IDDES model, the results simulated by the RANS model (SST k-x model) are also presented. This paper focuses on the power coefficients and torques of the H-VAWT, the pressure distributions on the blades, and the vortex structures of the rotor. It also tests the ability of the IDDES model to simulate the aerodynamics of the H-VAWTs. 2. IDDES method The IDDES method is adopted as the turbulence model to simulate the aerodynamics and performance of the H-VAWT. In the IDDES model, the modified Menters SST two-equation eddyviscosity [31] is used as the RANS model, which modifies the dissipation-rate term of the turbulent-kinetic energy (TKE) transport equation. Hence, the TKE equation for the IDDES model can be written as:

@ðqkÞ @ðquj kÞ @ þ ¼ @t @xj @xj

1 @k qk3=2 þ sij Sij l þ lt ; @xj r LIDDES

ð1Þ

where t, k, q, uj, l, lt, sij and Sij are the time, turbulent kinetic energy, density, velocity, molecular viscosity, turbulent viscosity, tensor of stress, and mean strain rate, respectively. The LIDDES defined as the IDDES length scale [32], is written as:

LIDDES ¼ ~f d ð1 þ f e ÞLRANS þ ð1 ~f d ÞLLES ;

ð2Þ

where the LIDDES and the LLES are the length scales of the RANS and the LES, respectively. They are defined as: 1=2

LLES ¼ C DES D;

LRANS ¼

k ; b x0

ð3Þ

where b⁄ is a constant in the SST k-x equal to 0.09, D = min{max {CwDmax, Cwd, Dmin}, Dmax} is the grid scale, Cw is an empirical constant, d is the distant to the nearest wall, Dmin = min{Dx, Dy, Dz} and D = max{Dx, Dy, Dz}, ~f ¼ maxfð1 f Þ; f g is interpreted max

d

dt

Fig. 1. Geometric model of an H-rotor VAWT: (a) 3D model; (b) top view of the rotor; (c) front view of the rotor.

B

with more details in Ref. [27]. 3. Numerical simulation Generally, the flow can be considered incompressible when the Mach number less than 0.3 [33]. As the Mach number of the present H-VAWT is 0.0832, the 3D unsteady incompressible NavierStokes equations are used as the governing equations. The IDDES model acts as the turbulence model. STAR-CCM+ (CD-adapco Group, USA) [34], a computational fluid dynamics solver based on the finite volume method, is employed to solve the N-S equations.

Fig. 2. The definition of blade pitch angle (b).

3.2. Computational domain 3.1. Geometric model The 3D geometrical model of the H-VAWTs is shown in Fig. 1(a). The airfoil section of the blades is NACA-0021 airfoil with a chord length of c = 0.265 m. The diameter of the rotor is D = 2 m (Fig. 1 (b) and (c)), and the blade span is H = 1.2 m (Fig. 1(c)). The blade pitch angle b (as shown in Fig. 2) is set as 4° in accordance with the experimental settings [35]. The solidity is defined as Nc/ D = 0.265 for the present H-VAWT, where N, c, and D are the number of blades, chord length, and rotor diameter respectively. All the other information of the wind turbine can be found in the experiments [35]. The shaft and the support arms are neglected in the simulations due to their ignorable influences [19].

In order to simulate the performance of the H-VAWTs, an appropriate 3D computational domain is constructed. A cross section of the computational domain perpendicular to the axis of rotation is shown in Fig. 3(a). The computational domain contains two parts: the rotating zone and the stationary zone, which is 6D in width and 15.5D in length [36]. The height of the computational domain is 4 times the blade span, as shown in Fig. 3(b). The whole 3D computational domain is shown in Fig. 3(c). A sliding mesh technique is used to simulate the rotating characteristic of the rotor. The contact surfaces of the rotating zone and the stationary zone are set as the interface, on which the two zones can exchange

238

H. Lei et al. / Energy Conversion and Management 133 (2017) 235–248

Fig. 4. Topology structure of computational grid around the blade: (a) boundary layer grids; (b) grids near blade; (c) grids on blade surface.

Fig. 3. Computational domain and boundary conditions specification: (a) top view; (b) front view; (c) 3D plot.

data after each time step. The mesh size of the two zones in the interface should be generated consistently as much as possible, so that the effectiveness and accuracy of the data exchange can be guaranteed. The boundary condition is shown in Fig. 3(a) and (b). In order to simulate the experimental condition and compare the numerical results to the experimental data, an interface is set on the contact surface of the rotating zone and the stationary zone. The blade’s surface, the inlet and the outlet are set as a no slip wall, velocity inlet (8 m/s), and pressure outlet with zero gauge respectively. Slip walls are set on the four sides of the stationary zone.

layer as it is able to effectively reduce the number of grids, reduce interpolation errors, and improve the solution efficiency [38]. The accuracy of the numerical simulation results is affected by the number of meshes significantly. Hence, four different types of meshes are constructed for comparison. As shown in Table 1, the power coefficients based on the IDDES model are listed at TSR = 2.24, which is corresponding to the maximum power coefficient [35]. The TSR can be defined as k = xD/2V1, where x is the angular velocity of the rotor, and V1 is the inflow velocity. With the increase of cell number, the power coefficients become higher and thus closer to the experimental data. As can be observed, the error between the predicted power coefficient and the experimental data reaches the smallest when the cell number is 2,089,128. At the same time, by further increasing the cell number to 2,516,288 (index 4), the predicted power coefficient only increases by 1% when compared with the index 3. Balancing the computational time and accuracy, mesh index 3 (2,089,128 cells) is adopted as the standard for the subsequent simulations.

3.3. Mesh dependency test 3.4. Solver settings The purpose of the current study is to investigate the performance of the wind turbine and the aerodynamic characteristics of the blades. Hence, the mesh around the blades should be refined enough. Meanwhile, prismatic boundary layer cells are needed to be made on the surface of the blades. The total thickness of the boundary layer is 0.025 m with a growth ratio of 1.25, as depicted in Fig. 4(a). The mesh topology around the blades is shown in Fig. 4, the inlet flow velocity is 8 m/s, and the Reynolds number base on the chord length (0.265 m) is about Re = 1.62 105–2.89 105. Hence, the height of the first prismatic cells is set to be 2 105 m which ensures y+ < 1 to meet the requirements of the DES model for grid resolution [37]. In addition, as shown in Fig. 4(b) and (c), an unstructured polyhedral cell is used to mesh the outside boundary

Implicit unsteady segregated flow method is chosen to solve the discretized continuity and momentum equations, and the secondorder difference formula is used for the discrete differential equations. The SIMPLE scheme is used to solve the pressure-velocity coupling equation and the IDDES turbulence model. The convergence criterion is set to 104, which means that the converged solutions are regarded when the residuals are equal to or less than 104. To improve the computational efficiency, 8 parallel CPU cores are selected to solve the equations. All the simulations are operated on a Small-Scale Server with two Intel(R) Xeon(R) CPUs (E52630 v3), 128 GB memory, and 2T hard drive. It takes 208 h to achieve six revolutions in calculation.

239

H. Lei et al. / Energy Conversion and Management 133 (2017) 235–248 Table 1 Mesh dependency test based on the errors of the power coefficient. Mesh index

Number of mesh

Predicted CP by IDDES

Experimental value

Relative error (%)

1 2 3 4

1,080,927 1,440,603 2,089,128 2,516,288

0.141 0.152 0.178 0.180

0.173 0.173 0.173 0.173

18.50 12.14 2.89 4.05

In the unsteady state calculation, the setting of the time step is a very critical problem [22]. There were many studies investigating the effect of the time step on the calculation results. Elkhoury et al. [21] studied a 3-bladed VAWT, the time step was set to the values corresponding to the time that the rotor rotated 2p/300 and 2p/600 rad, respectively, and they found that the two time steps had almost no effect on the instantaneous torque. Therefore, the time step is set to 0.001 s (which corresponds to time required for the rotor to rotate 2p/300 rad) in the present study so as to reduce the computational time. For obtaining accurate aerodynamic force and stable flow field information, the maximum physical time should be sufficiently long, such that the stable flow field information can be obtained after several rotations of the rotor. Fig. 5 shows the instantaneous torque values changing over time, which demonstrates that the adjacent peaks of the torque have almost no change (error less than 1%) after 4–5 rotations of the rotor. Therefore, the performance of the H-VAWT and the flow field are considered to be stabilized after 5 rotations, indicating that the maximum physical time is equal to 6 rotation periods of the rotor. 4. Validation of the CFD model One of the important parameters for a wind turbine is the Power coefficient (Cp), which is a measure of turbine efficiency and defined as [10]:

Qx ; q V 31 DH 2

Cp ¼ 1

ð4Þ

where Q, x, V1, D, and H are the rotor torque, angular velocity of the rotor, inflow velocity, rotor diameter, and blade span respectively. Here, the experimental results [35] are also utilized to compare the predicted Cp with the IDDES model (Cp1), and the SST k-x model (Cp2). The power coefficients are compared in Fig. 6(a). It is found that the variation trend of the power coefficients predicted by the two models both fit well with the experiments, and the values all approximate to the experimental data when the TSRs are less than 1.250. However, the values of Cp2 are smaller than the experimen-

Fig. 5. Periodic torque curve for one blade over several rotations of rotor.

Fig. 6. Validation of the CFD model between the simulations (IDDES model and SST k–x model) and the experiments [35]: (a) Comparison of power coefficients at different TSRs; (b) Comparison of average wake velocity at 1D.

tal data and show a large error when the TSRs are between 1.450 and 2.478. The values of Cp1 are larger than the values of Cp2 when the TSRs exceed 1.250, and the values of Cp1 are closer to the experimental data in the same grid resolution. In light of this comparison, the IDDES model is able to predict the power coefficients more accurately than the SST k-x model, and is regarded as a reliable method to investigate the performance of the VAWTs. It should be noted that the values of Cp1 are slightly higher than the experimental data at some points. Compared with the experiment, the maximum error, which occurs at the point of TSR = 2.478, is 9.03%. This may be caused by the neglect of the shaft and support arms. The same phenomena was also found by Siddiqui et al. [18], who reported the approximated error of 10–12% when the support arms was excluded and only the blades were left. Furthermore, the mechanical energy loss in the experiments was not considered in the simulations, causing the numerical results of Cp1 to be slightly higher. The wake velocity profiles are also compared between the experiments [35] and the simulations, as shown in Fig. 6(b). The

240

H. Lei et al. / Energy Conversion and Management 133 (2017) 235–248

IDDES model can more realistically reflect the average wake velocity, both at the lower TSR and the higher TSR. The SST k-x model shows a smaller wake velocity than the experiment at the lower TSR due to too much viscosity, which affects the velocity distribution. Meanwhile, the turbulence is regarded as isotropic, and the transport of momentum in the far field has not been correctly considered in the SST k-x model. Whereas, the large eddies in the far field are resolved using the LES method in the IDDES. Hence, the IDDES can more accurately predict the wake velocity than the SST k-x model can.

5. Results and discussions Torque is a major parameter in evaluating the generating capacity of a wind turbine [26]. The relationship between the torque and the power coefficient in Eq. (4) is positively related about the angular velocity. The torque ripples of one-bladed and two-bladed at TSR = 1.032 and TSR = 2.245 respectively, are plotted in Fig. 7, where the results of the SST k-x model were plotted as well. Fig. 7(a) and (c) shows the instantaneous torque ripples of one-bladed during one revolution. The maximum instantaneous torques occur at the azimuth angle of h = 74° (h is the azimuthal angle) at TSR = 1.032 in Fig. 7 (a), and azimuth angle of h = 90° at TSR = 2.245 in Fig. 7(c). The relative errors of the maximum instantaneous torques predicted by the two models are 11.0% (TSR = 1.032, Fig. 7(a)) and 17.1%

(TSR = 2.245, Fig. 7(c)). The regions of positive torque are all located in the upwind spacing. Most of the negative torque occurs in the downwind spacing due to the increase of the drag force. Fig. 7(b) and (d) shows the torque ripples for two-bladed torque ripples at TSR = 1.032 and TSR = 2.245, respectively. Two peak torques are observed during one revolution, which correspond to the effects of the two blades. Since the torque of the rotor is the vector sum of the torques of each blade, it can be seen that when the torque of one blade reaches the peak value, the torque of the other one is close to zero. Hence, the maximum values of the torque for one-bladed and two-bladed have almost no difference during one revolution (see Fig. 7(a)–(d)). It can also be found that the regions of negative extremal torque of the rotor appear twice during one revolution, and the rotation of the rotor is affected by the inertia in the two regions. Meanwhile, the locations of the negative torque regions and the peak torque change along with the TSRs. From Fig. 7, it is revealed that the average rotor torques obtained by the SST k- x model are still smaller than those of obtained by the IDDES model. However, the torque ripples simulated by the two models all show the same phenomenon. The torque ripples at a larger TSR (TSR = 2.245) seem more smooth than at a lower TSR (TSR = 1.032) in both one-bladed and two-bladed cases. Different from the HAWTs, the aerodynamic forces of VAWTs change periodically with the rotor’s rotation due to the unceasing variation of the Angle of attacks (AOAs). The AOA is defined as the angle between the inflow velocity and the chord line, which is

Fig. 7. Torque variation ripples of One-bladed and Two-bladed during one revolution: (a) One-bladed, TSR = 1.032; (b) Two-bladed, TSR = 1.032; (c) One-bladed, TSR = 2.245; (c) Two-bladed, TSR = 2.245.

H. Lei et al. / Energy Conversion and Management 133 (2017) 235–248

241

depicted in Fig. 8 and denoted as a0 with b = 0°, and a with b = 4°. The equations for the velocity of the wind relative to the rotating blade can be defined as: !

!

!

W ¼ V1 þ Vt; Vt ¼

xD 2

ð5Þ

;

ð6Þ

!

!

where V 1 and V t are the free stream wind speed and the blades rotating speed, respectively. The vector relation for them is also shown in Fig. 8. Hence, the W can be expressed as:

W ¼ V1

qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ 1 þ 2k cos h þ k2 ;

a0 ¼ tan1 a ¼ tan1

sin h ; k þ cos h

ð7Þ

sin h b; k þ cos h

ð8Þ

Fig. 9. Schematic of AOA (a) changes with azimuth angle (h) at different TSRs, (solid line: b = 4°, dotted line: b = 0°).

ð9Þ

maximum negative pressure are relatively large at almost 11.32%, and increase with the increase of the TSRs. Further, the pressure distributions obtained by the two models are more consistent when the TSR is low (TSR = 1.032), but become different when the TSR increases to 2.245, and especially show a larger difference on the suction side. The reason could be that some laminar separation bubbles emerge on the trailing edge and shed to the suction side. As shown in Fig. 10(a) and (b), the larger positive pressure occurs at the two ends of leading edge and a smaller positive pressure occurs at the middle of leading edge. Fig. 11 shows the 3D streamline of the two-bladed VAWT as simulated by the IDDES model at TSR = 1.032 and h = 74°. Since the pressure difference is very large between the pressure side and the suction side, an increasing pressure gradient exists between the two sides. It leads to the streamline of a blade moving from the top and bottom to the middle of the blade. Fig. 12 shows the velocity vector of one plane parallel to the inflow. The velocity decreases gradually when the air flows through the blades. The velocity in the downstream spacing is smaller than that in the upstream spacing. All of these flow regimes lead to the cross-flows appearing at the top and bottom of a blade, and a flow area with an apparent larger velocity at the two ends of a blade. Hence, the pressure at the top and bottom ends is relatively larger than that of the middle of the blade. Fig. 13 presents the pressure coefficients on the middle section of one blade at 30° intervals during one revolution, as predicted by both the IDDES model and the SST k-x model. It is found that the pressure coefficient curves predicted by the two models are almost coincident at the azimuth angles of h = 26° and h = 56°, as the vortices around the inner and outer surfaces of the blade predicted by the two models are basically identical at the azimuth angle of 26° (see Fig. 14(a) and (b) h = 26°). However, the pressure coefficient curves show differently at other azimuth angles. In Fig. 9, it can be found that the AOAs will be greater than 15° when the azimuthal angle exceeds 60° at TSR = 1.032, which means that the blade is experiencing the dynamic stall. At the moment, the flow separation has emerged and different scales of vortices are formed. As can be seen in Fig. 13, the pressure coefficients in trailing edge and inner face show a large difference at h = 86° and h = 116°. In the two azimuth angles, the trailing edge vortex forms and sheds to the inner face (see Fig. 14(c) and (d) h = 86°), and the laminar separation vortices are attributed to the laminar-turbulent transition on the surface of the blades. Since the laminar separation vortices have different effects on the distributions of velocity and pressure, the predicted pressure coefficients are affected by the ability of the IDDES and the SST k-x models to distinguish those vortices and as a result the different pressure coefficients form. When the blade

where k = xD/2V1 is the TSR, and h is the azimuth angle. According to Eq. (9), the relation between the AOA and the azimuth angle can be plotted in Fig. 9, which shows the increase of the absolute value of the AOA with the decrease of the TSRs. The setting of the pitch angle could reduce the values of the AOA for the airfoil in the upstream spacing, shrinking the range of dynamic stall. In fact, the AOAs of the airfoil in the downstream spacing are smaller than those plotted in Fig. 9, since the induced velocity in the downstream spacing is lower than the free stream velocity. It also leads to the fact that dynamic stall in the downstream spacing is alleviated. For the symmetrical airfoil blades, Sheldahl and Klimas [39] pointed out that the dynamic stall often occurred at the values of AOAs ranging from 10° to 15°, and the dynamic stall always took place at k < 4. Hence, the large AOA is a sufficient condition of dynamic stall, and it is necessary to study the aerodynamic forces of the blades when they operate at a lower TSR. From Fig. 9, the maximum AOA of a blade airfoil can reach about 80° and the dynamic stall always takes place at TSR = 1.032 for the present blades. In this section, the pressure distributions of one blade at the high AOA will be discussed. Fig. 10 shows the pressure distributions on one blade at the points of maximum torque errors predicted by the two different models in Fig. 7(a) and (c). The points occur at h = 74° with TSR = 1.032 and h = 90° with TSR = 2.245. It can be seen that the errors of maximum positive pressure predicted by the two turbulence models are small in the two cases. However, the errors of

Fig. 8. The definition and the schematic of AOA (a).

242

H. Lei et al. / Energy Conversion and Management 133 (2017) 235–248

Fig. 10. Contours of the instantaneous pressure distribution on the surface of one blade: (a) IDDES, TSR = 1.032, h = 74°; (b) SST k–x, TSR = 1.032, h = 74°; (c) IDDES, TSR = 2.245, h = 90°; (d) SST k–x, TSR = 2.245, h = 90°.

Fig. 11. Diagrams of streamline predicted by the IDDES at TSR = 1.032, h = 74°.

rotates to the downstream zone, the pressure coefficient curves (see Fig. 13, h = 266°, h = 296° and h = 326°) predicted by the two models still show some differences. The velocity and turbulence intensity of the inflow in the downstream zone are different from the upstream due to the disturbance of the blade in the upstream zone. It makes the aerodynamics of the blade in the downstream zone more complex than those in the upstream zone. In the downstream zone, the blade pitches downward, and the absolute values of the AOAs are diminished gradually. As a result of those reasons above, the vortices in the downstream zone obtained by the two models are more different relative to those in the upstream zone. As shown in Fig. 14(c) and (d), the vortex magnitude and the posi-

tion of the separation point predicted by the IDDES are different from those of the SST k-x model. It can also be seen that the pressure coefficients obtained by the two turbulent models still show a slight difference at h = 356°. In fact, the AOA corresponding to the azimuth angle of 356° is less than the stall AOA. However, the flow field requires some time to recover after a dynamic stall. The two models still show different vortices (see Fig. 14(g) and (h), h = 356°), which causes the difference in the predicted pressure coefficient. On the other hand, the SST k-x model is effective for the fully developed turbulence due to the limitation of the inherent characteristics. Although the Wall Function is used, it is still not desired.

H. Lei et al. / Energy Conversion and Management 133 (2017) 235–248

243

Fig. 12. Velocity vector of a plane predicted by the IDDES at TSR = 1.032, h = 74°.

In present study, the flow regimes of the H-VAWTs during one rotation are firstly characterized into three stages: attach flow, mildly flow separation and massively flow separation. In contrast to the SST k-x model, the IDDES could resolve all the regions of attached flows, mildly separated flows and massively separated flows, formation [27]. Hence, it shows a large difference in predicting the pressure coefficient curves by the two models (see Fig. 13, h = 86°, h = 266°) when the mildly separated flows occur (see Fig. 14, h = 86°, h = 266°), especially the discrepancy in the trailing edge. However, the pressure coefficient curves at h = 176° and h = 206° obtained by the two models demonstrate the slight differences relative to the case of mildly separated flows. This is because the massively separated flows take place (see Fig. 14, h = 176°, h = 206°), and the fully developed turbulence can be resolved by the two models. Therefore, the slight discrepancies in predicting the pressure coefficients are caused by the natural characteristics of the two models treatment of the N-S equations which result in the discrepant abilities to capture the vortices. Based on this, the IDDES model is regarded as a more accurate approach to deal with the aerodynamic problems of the H-VAWTs after flow separations formation relative to the SST k-x model. Thus it can be seen that the predicted pressure coefficients are affected by the ability of the turbulence models to cope with the flow separations. Almohammadi et al. [16] studied the pressure coefficients of a three-bladed H-VAWT using the transition SST model and the SST k-x model. It was found that the effect of the laminar separation bubbles could be suppressed if only the SST k-x model was used without the translation model, and the leading edge bubbles were only to be distinguished by the transition model. In their studies, the pressure coefficients of the airfoil obtained by the transition SST model and the SST k-x model were not in accordance when the azimuth angles ranged from 72° to 264°. The interval of the azimuth angle showing the accordance and discordance of the pressure coefficients in the present study is essentially consistent with the results of Ref. [16]. Fig. 14 shows the instantaneous vorticity magnitude of the rotor’s middle section at a lower TSR = 1.032 (larger AOA) predicted by the IDDES model (the left side) and the SST k-x model (the right side), respectively. As shown in Fig. 14(a) and (b), the results predicted by the two models both show that the large vortex is not formed around the blades at the position of 26°, the airflow attaches to the surface of the blades and the flow separation does not appear. It is clearly seen that the vortex structures around the blade predicted by the IDDES model are different from those predicted by the SST k-x model when the blade rotates to the azimuth angle of 206° (see Fig. 14(a) and (b), h = 206°). In the results of IDDES, there are two bubbles appearing on the leading edge and

trailing edge respectively. A series of different magnitudes of vortices around the airfoil are observed. This rich vortex structures are not found in the visualization of flow field obtained by the SST k-x model, where only two banded vortex structures are formed on the leading edge and trailing edge of the airfoil. This phenomenon strongly supports the difference of the pressure coefficients predicted by the two models as shown in Fig. 13 (h = 206°). When the blade rotates to the azimuth angle of 86° (see Fig. 14 (c) and (d), h = 86°), the dynamic stall develops further and the detached vortex derived from the trailing edge moves to the middle of the airfoil. Meanwhile, a reversed vortex at the trailing edge is induced again. As one can see, the results predicted by the two models all emerge as a large detached vortex on the suction surfaces of the blades, but the details of the detached vortex structures and magnitudes are not quite the same. There is a consecutive anticlockwise detached vortex on the suction surface of the IDDES model, while the detached vortex is inconsecutive in the SST k-x model. The similar phenomena also can be found when the blade rotates to the azimuth angle of 316° (see Fig. 14 (e) and (f), h = 316°). The AOA increases further when the blade rotates to the position of 136° as shown in Fig. 14(e) and (f). The energy of the induced reversed vortex at the trailing edge is enhanced gradually, therefore, two large bulges of vortices appear at the leading edge and the trailing edge, and the vortex at the leading edge is shedding and deviating from the leading edge to the inner side of the airfoil. Although the similar phenomenon can be observed in the SST k-x model, the strength of the vortex core is much smaller than that of the IDDES model (see Fig. 14(f), h = 136°). Subsequently, as shown in Fig. 14(g) and (h), the large scale flow separation appears in the surface of the airfoil at the position of 176° due to the AOA reaching 70°. The induced reversed vortex on the trailing edge and the leading edge detached vortex interfere with each other, which leads to the two vortices cyclical shedding from the two edges of the airfoil alternately. As the visualization obtained by the IDDES model shows (see Fig. 14(g), h = 176°), three different magnitudes of vortices separate from the leading edge, and a large vortex attaches to the leading edge and flows to the inner of the airfoil. This large vortex in the leading edge was not found in the results of the SST k-x model (see Fig. 14 (h), h = 176°). Meanwhile, the vortices shedding from the trailing edge obtained by the two models show a large difference in their position. As can be seen from Fig. 14, IDDES could show richer vortex structures and larger magnitudes of vorticity around the middle section of the rotor. This is because the LES, which engaged in the IDDE model out of the boundary zone, is able to capture the large scale vortex more effectively and obtain more sophisticated flow field information when compared with the RANS model.

244

H. Lei et al. / Energy Conversion and Management 133 (2017) 235–248

Fig. 13. Pressure coefficient curves of the middle section of one blade during one revolution.

Since the large scale vortex plays a main role in the momentum transport and energy transfer of the turbulence, the IDDES model could reflect the information flow field with more detail than that of the RANS model. Fig. 15 shows the Q-criterion iso-surface colored by the magnitude of the resultant velocity from the IDDES model and the SST k– x model for TSR = 1.032 and TSR = 2.245. Q-criterion is a method to describe the instantaneous vortex structures, which could

reflect the 3D structures of the flow field. The Q is a second invariant with respect to the velocity gradient tensor, defined as [37]:

Q¼

1 Xij Xij Sij Sij ; 2

ð10Þ

where

Xij ¼

1 @ui @uj ; 2 @xj @xi

Sij ¼

1 @ui @uj ; þ 2 @xj @xi

ð11Þ

H. Lei et al. / Energy Conversion and Management 133 (2017) 235–248

Fig. 14. Instantaneous contours of vorticity magnitude at different azimuth angles for TSR = 1.032, IDDES (the left side), SST k–x (the right side).

245

246

H. Lei et al. / Energy Conversion and Management 133 (2017) 235–248

(a)

(b)

(c)

(d) Fig. 15. Iso-surface contours of Q-criterion colorized by norm of resultant velocity (U) for Q = 250: (a) IDDES, TSR = 1.032; (b) SST k–x, TSR = 1.032; (c) IDDES, TSR = 2.245; (d) SST k–x, TSR = 2.245.

where Oij, Sij, and u are the vorticity, mean strain rate, and velocity respectively. As shown in Fig. 15(a) and (b), the 3D vortex structures around the blades predicted by the IDDES model are more plentiful than those predicted by the SST k-x model when the TSR is equal to 1.032. While the blade rotates to the azimuth angle of 26°, the wing tip vortices obtained by the two models are almost the same since the flow is attached to the surface of the blades, which all

show two trailing edge vortices on the two ends of the blades diffusing to the rear. The columnar shedding vortex can be seen from the results of the IDDES, but it does not appear in the results of the k-x model. A great difference appears in blade 2, which is in a state of dynamic stall. A lot of smaller vortices around blade 2 can be observed by the IDDES model, especially the smaller scale leading edge vortices, trailing edge shedding vortices and wing tip vortices. Whereas, those fine vortex structures are not observed obviously

H. Lei et al. / Energy Conversion and Management 133 (2017) 235–248

in the visualization field of the SST k-x model. Although wake away from the blades are all captured by the two models, the IDDES model could show more detailed vortex structures than the SST k-x model. Li et al. [20] pointed out that the VAWT energy extraction was affected by the vortex shedding from the blades in the upwind region severely when the dynamic stall occurs. The vorticity magnitude can reflect the energy transportation and dissipation of the flow field, which also affects the turbulent kinetic energy and the energy dissipation rate. Hence, the ability of the turbulence model to capture the vortices affects the prediction of aerodynamic characteristics directly. Fig. 15(c) and (d) demonstrate the vortex structures when the TSR is 2.245. Since the AOAs of the two blades are lower, the flow separation does not happen. Two strip vortices form in the two ends of the trailing edge, and a series of shedding vortices are developed moving toward the downstream. Hence, it can be found that the vortices obtained be the two models have almost no difference when the flow separation does not happen. From the 3D vortices structures, the IDDES model is able to identify more complex vortices shedding from the leading edge and trailing edge out of the boundary region. Those small scale vortices just have a crucial impact on the dissipation of energy, and thereby affect the pressure and velocity distribution in the flow field. Therefore, it can be explained again that the pressure prediction using the IDDES model is relatively more realistic than that of the SST k-x model when the blades are suffering the dynamic stall. As can be seen, the pressure coefficients in Fig. 13 are different between the two models. Hence, the IDDES is a promising model to analyze the performance and aerodynamics of VAWTs in the state of dynamic stall.

6. Conclusions The present study firstly employs the 3-D Navier-Stocks equations with the IDDES model and the polyhedral mesh to resolve the incompressible aerodynamic problems for an H-VAWT. In order to analyze the aerodynamics of VAWTs more carefully, both the results by the IDDES and the SST k-x model are shown. The main conclusions are listed as: (1) The power coefficients predicted by the IDDES model are closer to the experimental results than those of the SST kx model at the higher TSRs (between 1.450 and 2.478). Meanwhile, the IDDES can more realistically predict the wake velocity than that of the SST k-x at a lower TSR. (2) At a lower TSR (TSR = 1.032), the cross-flow appearing in the two ends of a blade leads to the higher pressure at the two ends of the leading edge and lower pressure at the middle of the leading edge. Since the inherent treatments of N-S equations of the IDDES and SST k-x models, the pressure distributions on the middle section of one blade present larger differences in mildly separated flow than those in massively separated flow. (3) By the analysis of the results of the flow field visualization, it is found that the vortices around the blades predicted by the IDDES model are richer and more realistic than those obtained by the SST k-x model, especially when the blades suffer the dynamic stall. Combined with the vortices and the inherent characteristics of the two models, this can explain the differences of the predicted pressure coefficients. Meanwhile, it can also support the idea that the IDDES model is a promising method to predict the performance and aerodynamics of the H-VAWTs.

247

Acknowledgements The authors are grateful for the financial support from the National Natural Science Foundation of China (Nos. 51278297 and 51679139) and the Major Program of the National Natural Science Foundation of China (No. 51490674), Research Program of Shanghai Leader Talent (No. 20) and Doctoral Disciplinary Special Research Project of Chinese Ministry of Education (No. 20130073110096).

References [1] Howell R, Qin N, Edwards J, Durrani N. Wind tunnel and numerical study of a small vertical axis wind turbine. Renew Energy 2010;35(2):412–22. [2] Cai X, Gu R, Pan P, Zhu J. Unsteady aerodynamics simulation of a full-scale horizontal axis wind turbine using CFD methodology. Energy Convers Manage 2016;112:146–56. [3] Ismail MF, Vijayaraghavan K. The effects of aerofoil profile modification on a vertical axis wind turbine performance. Energy 2015;80:20–31. [4] Borg M, Shires A, Collu M. Offshore floating vertical axis wind turbines, dynamics modelling state of the art. Part I: Aerodynamics. Renew Sustain Energy Rev 2014;39:1214–25. [5] Jin X, Zhao G, Gao K, Ju W. Darrieus vertical axis wind turbine: basic research methods. Renew Sustain Energy Rev 2015;42:212–25. [6] Islam M, Ting DSK, Fartaj A. Aerodynamic models for Darrieus-type straightbladed vertical axis wind turbines. Renew Sustain Energy Rev 2008;12 (4):1087–109. [7] Mohamed MH. Performance investigation of H-rotor Darrieus turbine with new airfoil shapes. Energy 2012;47(1):522–30. [8] Gharali K, Johnson DA. Dynamic stall simulation of a pitching airfoil under unsteady freestream velocity. J Fluids Struct 2013;42(4):228–44. [9] Buchner AJ, Lohry MW, Martinelli L, Soria J, Smits AJ. Dynamic stall in vertical axis wind turbines: comparing experiments and computations. J Wind Eng Ind Aerodyn 2015;146:163–71. [10] Posa A, Parker CM, Leftwich MC, Balaras E. Wake structure of a single vertical axis wind turbine. Int J Heat Fluid Flow. http://dx.doi.org/10.1016/j. ijheatfluidflow.2016.02.002. [11] Bianchini A, Balduzzi F, Ferrara G, Ferrari L. Virtual incidence effect on rotating airfoils in Darrieus wind turbines. Energy Convers Manage 2016;111:329–38. [12] Butbul J, MacPhee D, Beyene A. The impact of inertial forces on morphing wind turbine blade in vertical axis configuration. Energy Convers Manage 2015;91:54–62. [13] Shahizare B, Nik-Ghazali N, Chong WT, Tabatabaeikia S, Izadyar N, Esmaeilzadeh A. Novel investigation of the different Omni-direction-guidevane angles effects on the urban vertical axis wind turbine output power via three-dimensional numerical simulation. Energy Convers Manage 2016;117:206–17. [14] Wang Y, Sun X, Dong X, Zhu B, Huang D, Zheng Z. Numerical investigation on aerodynamic performance of a novel vertical axis wind turbine with adaptive blades. Energy Convers Manage 2016;108:275–86. [15] Danao LA, Edwards J, Eboibi O, Howell R. A numerical investigation into the influence of unsteady wind on the performance and aerodynamics of a vertical axis wind turbine. Appl Energy 2014;116:111–24. [16] Almohammadi KM, Ingham DB, Ma L, Pourkashanian M. Modeling dynamic stall of a straight blade vertical axis wind turbine. J Fluids Struct 2015;57:144–58. [17] Chowdhury AM, Akimoto H, Hara Y. Comparative CFD analysis of vertical axis wind turbine in upright and tilted configuration. Renew Energy 2016;85:327–37. [18] Siddiqui MS, Durrani N, Akhtar I. Quantification of the effects of geometric approximations on the performance of a vertical axis wind turbine. Renew Energy 2015;74:661–70. [19] Lee YT, Lim HC. Numerical study of the aerodynamic performance of a 500 W Darrieus-type vertical-axis wind turbine. Renew Energy 2015;83:407–15. [20] Li C, Zhu S, Xu YL, Xiao Y. 2.5 D large eddy simulation of vertical axis wind turbine in consideration of high angle of attack flow. Renew Energy 2013;51:317–30. [21] Elkhoury M, Kiwata T, Aoun E. Experimental and numerical investigation of a three-dimensional vertical-axis wind turbine with variable-pitch. J Wind Eng Ind Aerodyn 2015;139:111–23. [22] Spalart PR, Jou WH, Strelets M, Allmaras SR. Comments on the feasibility of LES for wings, and on a hybrid RANS/LES approach. Adv DNS/LES 1997;1:4–8. [23] Sagaut P. Large eddy simulation for incompressible flows: an introduction. Springer Science and Business Media; 2006. [24] Li Y, Paik KJ, Xing T, Carrica PM. Dynamic overset CFD simulations of wind turbine aerodynamics. Renew Energy 2012;37(1):285–98. [25] Johansen J, Sørensen NN, Michelsen JA, Schreck S. Detached-eddy simulation of 447 flow around the NREL phase-VI blade. ASME 2002 wind energy symposium. American Society of Mechanical Engineers; 2002. p. 106–14. [26] Balduzzi F, Bianchini A, Maleci R, Ferrara G, Ferrari L. Critical issues in the CFD simulation of Darrieus wind turbines. Renew Energy 2016;85:419–35.

248

H. Lei et al. / Energy Conversion and Management 133 (2017) 235–248

[27] Shur ML, Spalart PR, Strelets MK, Travin AK. A hybrid RANS-LES approach with delayed-DES and wall-modelled LES capabilities. Int J Heat Fluid Flow 2008;29 (6):1638–49. [28] Spalart PR, Deck S, Shur ML, Squires KD, Strelets MK, Travin A. A new version of detached-eddy simulation, resistant to ambiguous grid densities. Theoret Comput Fluid Dyn 2006;20(3):181–95. [29] Hu R, Wang L, Fu S. Improved Delayed Detached Eddy Simulation of flow structures behind a backward-facing step. In: 54th AIAA aerospace sciences meeting. p. 1104. [30] Xiao ZX, Luo KY. Improved delayed detached-eddy simulation of massive separation around triple cylinders. Acta Mech Sin 2015;31(6):799–816. [31] Menter FR. Two-equation eddy-viscosity turbulence models for engineering applications. AIAA J 1994;32(8):1598–605. [32] Xiao L, Xiao Z, Duan Z, Fu S. Improved-Delayed-Detached-Eddy simulation of cavity-induced transition in hypersonic boundary layer. Int J Heat Fluid Flow 2015;51:138–50. [33] Larose GL, D’Auteuil A. Effect of local air compressibility on the aerodynamics of rectangular prisms at Mach number below 0.3. 6th international colloquium on bluff body aerodynamics from 7/20/2008 to 7/24/2008.

[34] USER CD-adapco. Guide: Star-ccm+; 2009. [35] Li QA, Maeda T, Kamada Y, Murata J, Yamamoto M, Ogasawara T, et al. Study on power performance for straight-bladed vertical axis wind turbine by field and wind tunnel test. Renew Energy 2016;90:291–300. [36] Zhang LX, Liang YB, Liu XH, Jiao QF, Guo J. Aerodynamic performance prediction of straight-bladed vertical axis wind turbine based on CFD. Adv Mech Eng 2013;5:905379. [37] Ghasemian M, Nejat A. Aerodynamic noise prediction of a horizontal axis wind turbine using Improved Delayed Detached Eddy Simulation and acoustic analogy. Energy Convers Manage 2015;99:210–20. [38] Lanzafame R, Mauro S, Messina M. Wind turbine CFD modeling using a correlation-based transitional model. Renew Energy 2013;52:31–9. [39] Sheldahl RE, Klimas PC. Aerodynamic characteristics of seven symmetrical airfoil sections through 180-degree angle of attack for use in aerodynamic analysis of vertical axis wind turbines. Albuquerque, NM (USA): Sandia National Labs; 1981.

Copyright © 2021 COEK.INFO. All rights reserved.