The Darrieus wind turbine: Proposal for a new performance prediction model based on CFD (original) (raw)

The Darrieus wind turbine: Proposal for a new performance prediction model based on CFD

Marco Raciti Castelli, Alessandro Englaro, Ernesto Benini*
Department of Mechanical Engineering, University of Padova, Via Venezia, I, 35131 Padova, Italy

A R T I C L E I N F O

Article history:
Received 18 December 2010
Received in revised form
18 May 2011
Accepted 20 May 2011
Available online 29 June 2011

Keywords:

Darrieus wind turbine Vertical-Axis-Wind turbine Performance predition

A B S T R A C T

This paper presents a CFD model for the evaluation of energy performance and aerodynamic forces acting on a straight-bladed vertical-axis Darrieus wind turbine. The basic principles which are currently applied to BE-M theory for rotor performance prediction are transferred to the CFD code, allowing the correlation between flow geometric characteristics (such as blade angles of attack) and dynamic quantities (such as rotor torque and blade tangential and normal forces). The model is proposed as a powerful design and optimization tool for the development of new rotor architectures for which test data is not available.

After describing and validating the computational model against experimental data, a full campaign of simulation is proposed for a classical NACA 0021 three-bladed rotor.

Flow field characteristics are investigated for several values of tip speed ratio, allowing a comparison among rotor operation at optimum and lower CpC_{p} values, so that a better understanding of vertical-axis wind turbines basic physics is obtained.
© 2011 Elsevier Ltd. All rights reserved.

1. Introduction and background

The use of BE-M (Blade Element - Momentum) models for the design and analysis of vertical-axis wind turbines has aroused a large credit, not only in research and academic communities but also in industrial appliances, thanks to general acceptable accuracy of the results, as well as widely available literature and code simplicity. BE-M theory, first introduced by Glauert [1] in order to predict the structural dynamics and performance of airplane propellers, was adjusted for vertical-axis wind turbine aerodynamics by Templin [2], obtaining the most elementary approach based on momentum theory, named single-streamtube model. Strickland [3] extended Templin’s approach by considering the single streamtube as composed of a number of adjacent and aerodynamically independent smaller streamtubes, thus allowing the use of airfoil characteristics based on local (rather than on average) Reynolds number. This approach, named multiplestreamtube model, was further developed by several authors [4-9] by placing two actuator discs behind each other, thus combining the multiple streamtube model prediction tool with the

[1]double actuator disc theory. Through continuous refinements [10,11] (including also the most recent optimization techniques, such as genetic algorithms [12]), BE-M based design tools proved to be able of accurately predicting Darrieus wind turbine performance and are to be considered as the basis of most of the performance prediction tools currently used in wind turbine industry.

Aerodynamic performance prediction tools based on momentum theory, while presenting the advantage of a low computation effort, are nevertheless limited by the availability of extended airfoil databases (up to 180∘180^{\circ} ) which should also refer to rather low local Reynolds numbers (close to Re=105R e=10^{5} ), which are typical of a vertical axis wind turbine operation [13,14].

Since most airfoil database available in literature is derived from aeronautical applications, hardly extending far beyond stall and also referring to relatively high Reynolds numbers (above Re=106R e=10^{6} ), it has been recognized that there is a paucity of reliable airfoil data needed to describe the aerodynamics of wind turbine blades [15]. The airfoil section data requirements for application in vertical axis wind turbines are in fact broader in scope than those the aircraft industry usually concerns itself with. This limits the analysis on vertical axis wind turbines to a narrow number of blade profiles, substantially some NACA series of symmetrical airfoils, for which extended database at low Reynolds numbers is available in literature [16]. Moreover, as pointed out by Burton et al. [17], dwelling upon considerations on the energy extraction process rather than

1. • Corresponding author. Tel.: +39 049 8276767; fax: +39 049 8276785.
     E-mail addresses: marco.raciticastelli@unipd.it (M. Raciti Castelli), alessandro. englaro@studenti.unipd.it (A. Englaro), ernesto.benini@unipd.it (E. Benini). ↩︎

Nomenclature

Aij[359×90]A_{i j}[359 \times 90] matrix containing the local values of blade relative angle of attack as a function of both angular position and time
Ai[ mm2]A_{i}\left[\mathrm{~mm}^{2}\right] rotor swept area
Aii[1×90]A_{i i}[1 \times 90] vector containing the local values of blade relative angle of attack at time step ti as a function of azimuthal position j
A.O.A. [rad] average value of blade relative angle of attack
c[ mm]c[\mathrm{~mm}] blade chord
cP(θ)[−]c_{P}(\theta)[-] differential pressure coefficient
CP(θ)[−]C_{P}(\theta)[-] rotor instantaneous power coefficient
Cp, ave [−]C_{p, \text { ave }}[-] average power coefficient during a full rotor revolution
Ci(θ)C_{i}(\theta) [-] rotor instantaneous thrust coefficient
CT(θ)[−]C_{T}(\theta)[-] rotor instantaneous torque coefficient
Drotor [mm]D_{\text {rotor }}[\mathrm{mm}] rotor diameter
Fa(θ)[N]F_{a}(\theta)[\mathrm{N}] instantaneous rotor force in xx direction
Hdomain [mm]H_{\text {domain }}[\mathrm{mm}] computational domain height
Hrotor [mm]H_{\text {rotor }}[\mathrm{mm}] rotor height
Hwindtunnel [mm]H_{\text {windtunnel }}[\mathrm{mm}] computational domain height, validation model
j[−]j[-] azimuthal position along blades trajectory
n[−]n[-] number of azimuthal positions which are not considered for the calculation of the relative value of local angle of attack
N[−]N[-] number of rotor blades
p(x,y)p(x, y) total pressure at point of coordinates (x,y)(x, y)
P(θ)[W]P(\theta)[\mathrm{W}] rotor instantaneous power
Pave [W]P_{\text {ave }}[\mathrm{W}] average power during a full rotor revolution
Rrotor [m]R_{\text {rotor }}[\mathrm{m}] rotor radius
Re⁡[−]\operatorname{Re}[-] blade Reynolds numbers
Sdomain [mm2]S_{\text {domain }}\left[\mathrm{mm}^{2}\right] computational domain section
t[ s]t[\mathrm{~s}] \quad simulation time
ti[−]t_{i}[-] \quad simulation time step
T[ s]T[\mathrm{~s}] \quad rotor period of revolution
Trotor (θ)[Nm]T_{\text {rotor }}(\theta)[\mathrm{Nm}] rotor instantaneous torque
TSR [-] tip speed ratio
U[m/s]\mathrm{U}[\mathrm{m} / \mathrm{s}] absolute wind velocity at blade position
Ue[ m/s]U_{e}[\mathrm{~m} / \mathrm{s}] absolute wind velocity at blade position, component along xx axis
Uf[ m/s]U_{f}[\mathrm{~m} / \mathrm{s}] absolute wind velocity at blade position, component along yy axis
Ut[ m/s]U_{t}[\mathrm{~m} / \mathrm{s}] blade tangential speed at blade position
Vtest section [m/s]V_{\text {test section }}[\mathrm{m} / \mathrm{s}] mean wind velocity at rotor test section
V∞[m/s]V_{\infty}[\mathrm{m} / \mathrm{s}] unperturbed wind velocity at computational domain entrance
w[ m/s]w[\mathrm{~m} / \mathrm{s}] relative wind velocity at blade position
Wdomain [mm]W_{\text {domain }}[\mathrm{mm}] computational domain width
Wwind tunnel [mm]W_{\text {wind tunnel }}[\mathrm{mm}] computational domain width, validation model
α\alpha [rad] blade relative angle of attack in the plane of the airfoil cross section
αpen [rad]\alpha_{\text {pen }}[\mathrm{rad}] positive stall angle for NACA 0021 at Re⁡=300.000\operatorname{Re}=300.000 (assumed 0.21 rad)
αneg [rad]\alpha_{\text {neg }}[\mathrm{rad}] negative stall angle for NACA 0021 at Re⁡=300.000\operatorname{Re}=300.000 (assumed -0.21 rad)
αtij[rad]\alpha_{t i j}[\mathrm{rad}] local values of blade relative angle of attack as a function of azimuthal position j along blades trajectory for time step i
εsb[−]\varepsilon_{s b}[-] solid blockage correction factor
δ[[⋅]\delta\left[{ }^{[\cdot}\right] angle between absolute wind velocity at blade position and unperturbed wind direction
θ[[⋅]\theta\left[{ }^{[\cdot}\right] blade azimuthal coordinate
Ψ[rad]\Psi[\mathrm{rad}] angular sector covered during a period of rotor revolution
γ[∘]\gamma\left[{ }^{\circ}\right] angle between absolute wind velocity at blade position and blade translational speed at blade position
ρ[kg/m3]\rho\left[\mathrm{kg} / \mathrm{m}^{3}\right] air density (assumed 1.225 kg/m31.225 \mathrm{~kg} / \mathrm{m}^{3} )
σ[−]\sigma[-] \quad rotor solidity
ω[rad/s]\omega[\mathrm{rad} / \mathrm{s}] rotor angular velocity
on the specific turbine design, all classical aerodynamic tools are unable to visualize the basic structure of the flow field inside the rotor volume.

The limitations in low Reynolds number accurate aerodynamic databases can be overcome by using advanced CFD (Computational Fluid Dynamic ) codes, which can outflank the lack of airfoil data thanks to their inherent ability to determine the aerodynamic components of actions through the integration of the Navier-Stokes equations in the neighborhood of the wind turbine blade profile.

A first attempt of improving the efficiency of a vertical-axis wind turbine through the application of the emerging CFD capabilities was performed by Vassberg et al. [18] through the simulation of the dynamic motion of a turbine blade spinning about a vertical axis and subjected to a far-field uniform free-stream velocity flow field.

Also Ferreira et al. [19,20] presented a systematic CFD analysis of a two-dimensional straight-bladed rotor configuration. The effect of dynamic stall in a 2D single-bladed VAWT (Vertical-Axis Wind Turbine) was investigated, reporting the influence of the turbulence model in the simulation of the vortical structures spread from the blade itself.

Recently, Kumar et al. [21] proposed a low Reynolds number VAWT design and optimization procedure based on both CFD simulations and BE-M calculations, while Raciti Castelli and Benini [22] presented a model for the evaluation of energy performance and aerodynamic forces acting on a helical single-bladed VAWT
depending on blade inclination angle: the analysis was based on five machine architectures, which were characterized by an inclination of the blades with respect to the horizontal plane in order to generate a phase shift angle of 0∘,30∘,60∘,90∘0^{\circ}, 30^{\circ}, 60^{\circ}, 90^{\circ} and 120∘120^{\circ} between lower and upper blade sections, for a rotor with an aspect ratio of 1.5 .

D’Alessandro et al. [23] developed a new computational approach, based on both CFD and wind tunnel measurements, for the simulation of a drag-driven VAWT energy performance, in order to gain an insight into the complex flow field structures developing around a Savonius rotor.

Raciti Castelli et al. [24] performed a numerical analysis validation campaign for a Darrieus micro-VAWT through a systematic comparison with wind tunnel experimental data. This work proved that it is possible to determine the best near-blade grid element dimension through statistical analysis of some indicators, such as the y∘y^{\circ} parameter, in order to maximize the accuracy of the numerical prediction of rotor performance while maintaining a reasonable computational effort.

The main goal of the present work is to propose a new straightbladed VAWT performance prediction model, based on numerical simulations, in order to determine the rotor power curve. The basic principles which are currently applied to the BE-M theory for rotor performance prediction are transferred to the CFD code, allowing the correlation between flow geometric characteristics (such as blade angles of attack) and dynamic quantities (such as rotor torque and power). In particular:

• the kinematic and dynamic quantities which are classically determined from the “Momentum” section of a BE-M code have been calculated using CFD, thus being able to better investigate the flow field close to the rotor by considering the effects of dynamic stall, streamtube expansion and velocity components normal to unperturbed wind direction;
• the kinematic and dynamic quantities which are classically determined from the “Blade Element” section of a BE-M code have been calculated using CFD, thus being able to overcome the lack of low Reynolds number accurate aerodynamic databases thanks to the inherent ability of the CFD code in determining the aerodynamic components of actions through the integration of the Navier-Stokes equations in the neighborhood of the blade profile [25].

After describing the computational model and the relative validation procedure, a full campaign of simulation is proposed for a classical NACA 0021 three-bladed rotor. Flow field characteristics are investigated for several values of tip speed ratio, allowing a comparison among rotor operation both at optimum and lower values of angular velocity, so that a better understanding of vertical-axis wind turbines basic physics is obtained.

2. Model geometry

Fig. 1 shows a schematic description of the adopted survey methodology, consisting in the coupling of an analytical code to a solid modeling software (capable of generating the desired rotor geometry depending on the desired design geometric parameters) which is linked to a finite volume CFD code for the calculation of rotor performance. CFD results are post-processed using a second analytical code for calculation of rotor main kinematic and dynamic characteristics.

The aim of the present work was to numerically analyze the aerodynamic behavior of a three-bladed Darrieus VAWT operating at different angular velocities for a constant wind speed of 9 m/s9 \mathrm{~m} / \mathrm{s}. The main geometrical features of the tested rotor are summarized in Table 1. The solidity parameter σ\sigma is defined as Nc/Rrotor \mathrm{Nc} / \mathrm{R}_{\text {rotor }}, as suggested by Strickland [3].

Fig. 1. Scheme of the survey methodology.

Table 1
Main geometrical features of the tested model.

Rotor azimuthal position was identified by the angular coordinate of the pressure centre of blade No. 1 midsection (set at 0.25 - c for NACA 0021 airfoil), starting between the 2nd and 3rd Cartesian plane octants, as can be seen in Fig. 2.

3. Spatial domain discretization

As the aim of the present work was to reproduce the operation of a rotating machine, the use of moving sub-grids was necessary. In particular, the discretization of the computational domain into macro-areas has led to two distinct sub-grids:

• a rectangular outer zone, determining the overall calculation domain, with a circular opening centered on the turbine rotational axis, which was identified as Wind Tunnel sub-grid, fixed;
• a circular inner zone, which was identified as Rotor sub-grid, rotating with rotor angular velocity ω\omega.

3.1. Wind tunnel sub-grid

Fig. 3 shows the main dimensions and the boundary conditions of the Wind Tunnel sub-grid area.

The computational domain width was set almost 80 rotor diameters so that the value of solid blockage correction factor [26] was less than 0.32%0.32 \% and would cause an increase in the cube of the velocity at rotor test section of less than 1%1 \% with respect to the value of unperturbed wind velocity at computational domain entrance, in formulas:
εzb=1/4AcSdomin =1/4Drotor ⋅Hrotor Wdomain ⋅Hdomain \varepsilon_{\mathrm{zb}}=1 / 4 \frac{A_{\mathrm{c}}}{S_{\text {domin }}}=1 / 4 \frac{D_{\text {rotor }} \cdot H_{\text {rotor }}}{W_{\text {domain }} \cdot H_{\text {domain }}}
Vtest section =(1+εzb)V∞V_{\text {test section }}=\left(1+\varepsilon_{\mathrm{zb}}\right) V_{\infty}

Fig. 2. Azimuthal coordinate of blade midsection’s centre of pressure.

Fig. 3. Main dimensions of the Wind Tunnel sub-grid area (dimensions in mm).

Since the estimation of the correct value of wake blockage was rather difficult, and being the overall calculation domain over-sized by imposing the requirement on the cube of the velocity at the test section, the wake blockage was not considered in these calculations.

In order to allow a full development of the wake, Ferreira et al. [19] placed inlet and outlet boundary conditions respectively 10 diameters upwind and 14 diameters downwind with respect to rotor test section for a wind tunnel CFD simulation. As the aim of the present work is the simulation of a turbine operating in open field conditions and because of the huge domain width necessary to avoid solid blockage (as discussed above), inlet and outlet boundary conditions were placed respectively 37 rotor diameters upwind and 60 rotor diameters downwind with respect to rotor test section.

Two symmetry boundary conditions were used for the two side walls. The circumference around the circular opening, centered on the turbine rotational axis, was set as an interface, thus ensuring the continuity in the flow field.

An unstructured mesh was chosen for the Wind Tunnel sub-grid, in order to reduce engineering time to prepare the CFD simulations.

3.2. Rotor sub-grid

The Rotor sub-grid is the fluid area simulating the revolution of the wind turbine and is therefore characterized by a moving mesh, rotating at the same angular velocity of the turbine. Its location coincides exactly with the circular opening inside the Wind Tunnel sub-grid area and centered on the turbine rotational axis.

Fig. 4 shows the main dimensions and the boundary conditions of the Rotor sub-grid area.

It is good engineering practice to provide that the mesh on both sides of the interface (Rotor sub-grid and Wind Tunnel sub-grid areas) has approximately the same characteristic cell size in order to obtain faster convergence [27].

An isotropic unstructured mesh was chosen for the Rotor subgrid, in order to guarantee the same accuracy in the prediction of rotor’s performance during the revolution - according to the studies of Commings et al. [28] - and also in order to test the prediction capability of a very simple grid. Considering their features of flexibility and adaption capability, unstructured meshes are in fact very easy to obtain, for complex geometries, too, and

Fig. 4. Schema of Rotor sub-grid area (dimensions in mm).
often represent the “first attempt” in order to get a quick response from the CFD in engineering work.

All blade profiles inside the Rotor sub-grid area were enclosed in a control circle of 400 mm diameter. Unlike the interface, it had no physical significance: its aim was to allow a precise dimensional control of the grid elements in the area close to rotor blades, by adopting a first size function operating from the blade profile to the control circle itself and a second size function operating from the control circle to the whole Rotor sub-grid area, ending with grid elements of the same size of the corresponding Wind tunnel subgrid elements. An interior boundary condition was used for control circle borders, thus ensuring the continuity of the cells on both sides of the mesh.

A growth factor of 1.2 was set from the surface of the control circle to Rotor sub-grid, thus expanding grid size from 4 mm to 10 mm , as shown in Fig. 5.

Fig. 5. Rotor sub-grid mesh.

Fig. 6. Size functions applied to control circle elements.

4. The Control circle

Being the area close to the blade profiles, great attention was placed to the control circle. The computational grids were constructed from lower topologies to higher ones, adopting appropriate size functions, in order to cluster grid points near the leading edge and the trailing edge of the blade profile, so as to improve the CFD code capability of determining lift, drag and the separation of the flow from the blade itself.

Two size functions were set inside the control circle, as shown in Fig. 6:

• size function No. 1 started from the leading edge and influenced inner and outer blade profile;
• size function No. 2 started from the blade profile and influenced the whole control circle area, as described in the previous section.

Table 2 summarizes the main features of the mesh close to rotor blade.

The trailing edge thickness was set to about 0.38 mm . The trailing edge itself was meshed using 10 grid points, as shown in Fig. 7, while Fig. 8 displays a view of the control circle grid.

5. Code validation

Before analyzing the models described in the previous sections, a complete validation work based on wind tunnel measurements has been conducted. As shown in Fig. 9, the experimental setup consisted in a classical vertical-bladed Darrieus rotor made of aluminum and carbon fibers, using a NACA 0021 blade profile with a chord length of 85.8 mm , which was tested in Bovisa’s low turbulence facility (Milan).

Table 2
Main features of the mesh close to rotor blade.

Fig. 7. Grid points clustering close to trailing edge.

As can be seen from Fig. 10, a computational domain of rectangular shape has been chosen, having the same wind tunnel external sizes: the wall boundary conditions of the model consisted in two lateral walls spaced 2000 mm apart from the wind tunnel centerline (the wind tunnel measured 4000 mm in width and 3880 mm in height). The rotor axis was placed on the symmetry position of the wind tunnel section. 2D and 3D simulations were performed, in order to take into account the drag effect induced by the spokes, too. Only half of the experimental setup was modeled, due to its vertical symmetry: in this case a symmetry boundary condition was used at the computational domain midsection. Anyway, the geometrical features of the model did not allow other simplifications to be performed. The effect of gravity on the rotor working curves has not been contemplated, being considered not influential for the scope of the validation work.

The main features of the validation model are summarized in Table 3.

The correction due to wind tunnel blockage was not applied (both for experimental measurements and numerical computations), in order to minimize any sources of error due to a wrong

Fig. 8. Control circle grid for NACA 0021 blade section.

Fig. 9. Experimental setup for validation model.
estimation of the blockage of the wind tunnel itself. Furthermore, this choice had the significant advantage of reducing the computational domain, allowing a saving in the total number of mesh elements. The correction of the friction resistive torque due to the bearings was taken into account.

The 2D grid independence of the solutions was verified through simulations performed with grid sizes that varied from about 600,000 to more than 1,000,000 mesh elements. In order to test the code sensitivity to the number of grid points, three unstructured meshes were adopted for the Rotor sub-grid, while the Wind Tunnel sub-grid remained substantially the same.

Fig. 10. Validation model, 3D computational domain.

Table 3
Validation model main features.

After some corrections to take into account spoke drag, the average torque values, measured in the wind tunnel for a 9 m/s9 \mathrm{~m} / \mathrm{s} unperturbed wind speed at test section entrance (and different tip speed ratios), were compared with those obtained from CFD analysis. For more details about the validation procedure and tests, see Ref. [24].

Comparison of the computed results with experimental data showed that prediction obtained using Enhanced Wall Treatment k- ε\varepsilon Realizable turbulence model successfully reproduced most of the flow features associated with the revolution of the tested rotor. In particular, as can be seen in Fig. 11, the numerical code proved able to replicate the shape of the experimental curve and was able to accurately capture the maximum power coefficient tip speed ratio, thus offering a reliable alternative to the development of experimental tests, at least at a first attempt.

The discrepancies between the two curves, roughly constant over rotor operational angular range and corresponding to about one half of two-dimensional power coefficient for optimum tip speed ratio, are due to the combined effects of finite blade length and spokes drag, as described by Raciti Castelli et al. [24].

6. Temporal discretization and convergence criteria

The numerical analysis proposed in the present work kept some common points with the validation model, particularly as far as the blade profile and the rotor radius are concerned. Only the computational domain was enlarged in order to analyze the open field operation of the turbine. The adopted commercial CFD package was Fluent 6.3.26, that implements 2-D Reynolds-averaged Navier-Stokes equations using a finite volume-finite element based solver. The fluid was assumed to be incompressible, being the maximum fluid velocity in the order of 60 m/s60 \mathrm{~m} / \mathrm{s}.

As pointed out by McMullen et al. [29], the calculation of unsteady flows in turbomachinery continues to present a severe challenge to CFD. During VAWT operation, the unsteadiness stems mainly from the relative motion of the rotating blade fields and has a fundamental period which depends both on the rate of rotation and the number of blade passage. For the proposed calculations, the temporal discretization was achieved by imposing a physical time step equal to the lapse of time the rotor takes to make a 1∘1^{\circ} rotation. An improved temporal discretization simulation did not show any significant variation [24].

As a global convergence criterion, each simulation was run until instantaneous torque coefficient values showed a deviation of less than 1%1 \% compared with the equivalent values of the previous period, corresponding to a revolution of 120∘120^{\circ} due to rotor threebladed geometry. Residuals convergence criterion for each physical time step was set to 10−510^{-5}.

The present simulations required about 4 CPU seconds per physical time step. An average of about 25 sub iterations was necessary to make the solution converge at each physical time step. The calculations, performed on an 8 processor, 2.33 GHz clock

Fig. 11. Comparison between wind tunnel measured (red) and numerically determined 2D power curve (blue) for a wind speed at test section entrance of 9 m/s9 \mathrm{~m} / \mathrm{s} (correction due to wind tunnel blockage was not considered). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
frequency computer, required a total CPU time of about 4 days for each simulation.

7. Performance aerodynamic model

The proposed performance analysis is based upon a simplified aerodynamic model, consisting in the analysis of kinematic and dynamic quantities every 4∘4^{\circ} rotor azimuthal position along blades trajectory, as exemplified in Fig. 12.

For each azimuthal position, the main aerodynamic components of actions (such as rotor torque, power and the forces acting on the blade profiles) can be directly determined by the CFD code through the integration of the Navier-Stokes equations in the neighborhood of the wind turbine blade profile. Blade relative velocity is determined as the vector sum of the absolute wind velocity UU and the inverse of blade tangential speed −Ut-U_{t} at blade position, being blade position identified by the coordinate of its centre of pressure, as exemplified in Fig. 13.

Being θ\theta blade azimuthal coordinate and δ\delta the angle between absolute wind velocity at blade position and unperturbed wind direction, the angle γ\gamma between absolute wind velocity at blade position and the inverse of blade translational speed at blade position can be written as:
γ=θ−δ\gamma=\theta-\delta

Fig. 12. Reference azimuthal position along blades trajectory.

Fig. 13. Blade velocity vectors. The following blade characteristic velocities are shown: absolute wind velocity at blade position UU (red), inverse of blade tangential speed at blade position −Ut-U_{t} (blue) and relative wind velocity w at blade position (green). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
where δ\delta is given by:
δ=arctg⁡(Uy/Ux)\delta=\operatorname{arctg}\left(U_{y} / U_{x}\right)
being UxU_{x} and UxU_{x} the components of absolute wind velocity at blade position along xx and yy axes.

The angle between the relative velocity and the tangent to the blade trajectory at each azimuthal position is therefore given by:
α=arctg⁡Usin⁡γUcos⁡γ+(−Ut)\alpha=\operatorname{arctg} \frac{U \sin \gamma}{U \cos \gamma+\left(-U_{t}\right)}
and is named as blade relative angle of attack in the plane of the airfoil cross section.

As reported in Fig. 12, 90 reference azimuthal positions (spaced 4∘4^{\circ} each other) were identified along the circumference representing the trajectory of the blades. For each position and for each temporal time step the UxU_{x} and UyU_{y} components were determined, allowing the calculation of the characteristics blade angles δ,γ\delta, \gamma and θ\theta.

Being the value of blade tangential speed:
Ut=αRrotor U_{t}=\alpha R_{\text {rotor }}
it is possible to determine the values of the relative angle of attack α\alpha from Eq. (5) as a function of both blade azimuthal position and time. For each time step tht_{h}, the vector containing the local values of blade relative angle of attack as a function of azimuthal position jj along blades trajectory can be written as:
Ati=[αti,0−]=[αti,0−;αti,4−;αti,8−;…;αti,356−]A_{t i}=\left[\alpha_{t i, 0^{-}}\right]=\left[\alpha_{t i, 0^{-}} ; \alpha_{t i, 4^{-}} ; \alpha_{t i, 8^{-}} ; \ldots ; \alpha_{t i, 356^{-}}\right]
ranging jj from 0∘0^{\circ} to 356∘356^{\circ}. Some of αti,j\alpha_{t i, j} values are not determinable because of the passage of the rotor blade (there is no flow field at the considered azimuthal position along blades trajectory). These values are registered as “not a number” by the code and are not considered in the calculation.

By composing all AtiA_{t i} vectors during a full rotor revolution, a matrix can be obtained, which is representative of blade relative angles of attack during rotor operation as a function of both angular position and time:
At,j=[αt1,0−αt1,4−αt1,8−αt1,16−…αt1,356−αt2,0−αt2,4−αt2,8−αt2,16−…αt2,356−αt3,0−αt3,4−αt3,8−αt3,16−…αt3,356−………………αtn,4−αtn,8−αtn,16−αtn,0−…αtn,356−]A_{t, j}=\left[\begin{array}{cccccc}\alpha_{t 1,0^{-}} & \alpha_{t 1,4^{-}} & \alpha_{t 1,8^{-}} & \alpha_{t 1,16^{-}} & \ldots & \alpha_{t 1,356^{-}} \\ \alpha_{t 2,0^{-}} & \alpha_{t 2,4^{-}} & \alpha_{t 2,8^{-}} & \alpha_{t 2,16^{-}} & \ldots & \alpha_{t 2,356^{-}} \\ \alpha_{t 3,0^{-}} & \alpha_{t 3,4^{-}} & \alpha_{t 3,8^{-}} & \alpha_{t 3,16^{-}} & \ldots & \alpha_{t 3,356^{-}} \\ \ldots & \ldots & \ldots & \ldots & \ldots & \ldots \\ \alpha_{t n, 4^{-}} & \alpha_{t n, 8^{-}} & \alpha_{t n, 16^{-}} & \alpha_{t n, 0^{-}} & \ldots & \alpha_{t n, 356^{-}}\end{array}\right] being t1−t 1- th the number of time steps contained in a full rotor revolution.

Figs. 14 and 15 display the evolution of blade relative angle of attack as a function of time for two different azimuthal positions along blade trajectory.

Fig. 14. Evolution of the instantaneous values of blade relative angle of attack as a function of time for 40∘40^{\circ} azimuthal position, the red rectangles evidence the number of azimuthal positions which are not considered for the calculation; TSR =2.33=2.33.

The passage of the rotor blade is visible as a peak in the instantaneous value of blade relative angle of attack. The average value of local angle of attack, due to the presence of the rotating turbine inside the flow field at each time step, was determined by eliminating n azimuthal positions ( 4 in the present case, for a total angular spacing of 16∘16^{\circ} ) before and after the peak values (where the flow field is locally perturbed by the passage of the rotor blade) and by calculating the arithmetic mean of the remaining angular values.

The width of the angular sector perturbed by the local passage of rotor blade is a function of both rotor geometrical features (such as blade chord and, therefore, rotor solidity) and spatial discretization of the numerical analysis (such as the chosen number of azimuthal position analyzed along blades trajectory). From Figs. 14 and 15 and Table 4 it can also be seen that, being the distribution of blade relative angles of attack evenly spaced around their mean value (with the exclusion of the perturbed angular sectors), the proposed methodology is independent from the chosen value of nn.

8. Results and discussion

Fig. 16 represents the evolution of rotor average power coefficient, defined as:

Fig. 15. Evolution of the instantaneous values of blade relative angle of attack as a function of time for 96∘96^{\circ} azimuthal position, the red rectangles evidence the number of azimuthal positions which are not considered for the calculation; TSR =2.33=2.33.

Table 4
Resulting average values of local angle of attack for 40∘40^{\circ} azimuthal position as a function of the number of azimuthal positions (both before and after the peak values) which are not considered for the present calculation.

Cp, ave =Pave 1/2ρAbV∞2C_{p, \text { ave }}=\frac{P_{\text {ave }}}{1 / 2 \rho A_{b} V_{\infty}^{2}}
for an incident wind speed of 9 m/s9 \mathrm{~m} / \mathrm{s} and as a function of the tip speed ratio, defined as:
TSR⁡=ωRrotor /V∞\operatorname{TSR}=\omega R_{\text {rotor }} / V_{\infty}
Figs. 17-24 show the distribution of the instantaneous torque coefficient, defined as:
CT(θ)=Trotor (θ)1/2ρAbRrotor V∞2C_{T}(\theta)=\frac{T_{\text {rotor }}(\theta)}{1 / 2 \rho A_{b} R_{\text {rotor }} V_{\infty}^{2}}
and of the relative angle of attack, determined as described in the previous section, as a function of azimuthal position for a single rotor blade and for TSR ranging from 1.44 to 3.30 . Also both positive and negative stall angles for a NACA 0021 airfoil at Re=300,000R e=300,000 are represented.

As can be clearly seen, blade relative angles of attack decrease passing from lower to higher TSR values, because of the raising influence of blade translational speed in the near-blade flow field.

It can also be seen that the maximum torque values are generated during the upwind revolution of the turbine and for azimuthal position where rotor blades are experiencing very high relative angles of attack, even beyond the stall limit, as can be seen in Figs. 25 and 26, representing relative velocity contours and relative pathlines for the angular position (92∘)\left(92^{\circ}\right) of maximum instantaneous torque coefficient (and therefore power, being constant the angular velocity) and for the optimum value of TSR (2.33).

As can be seen from Fig. 20, the azimuthal positions of maximum power extraction along blade trajectory are located

Fig. 16. Evolution of rotor average power coefficient as a function of TSR for an incident wind speed of 9 m/s9 \mathrm{~m} / \mathrm{s}; rotor operation at 8 different TSR values was numerically simulated.

Fig. 17. Evolution of instantaneous torque coefficient and relative blade angle of attack as a function of blade azimuthal position for a single rotor blade; TSR=1.44\mathrm{TSR}=1.44.

Fig. 18. Evolution of instantaneous torque coefficient and relative blade angle of attack as a function of blade azimuthal position for a single rotor blade; TSR=1.68\mathrm{TSR}=1.68.

Fig. 19. Evolution of instantaneous torque coefficient and relative blade angle of attack as a function of blade azimuthal position for a single rotor blade; TSR =2.04=2.04.

Fig. 20. Evolution of instantaneous torque coefficient and relative blade angle of attack as a function of blade azimuthal position for a single rotor blade; TSR=2.33\mathrm{TSR}=2.33.

Fig. 21. Evolution of instantaneous torque coefficient and relative blade angle of attack as a function of blade azimuthal position for a single rotor blade; TSR=2.51\mathrm{TSR}=2.51.

Fig. 22. Evolution of instantaneous torque coefficient and relative blade angle of attack as a function of blade azimuthal position for a single rotor blade; TSR=2.64\mathrm{TSR}=2.64.

Fig. 23. Evolution of instantaneous torque coefficient and relative blade angle of attack as a function of blade azimuthal position for a single rotor blade; TSR=3.09\mathrm{TSR}=3.09.

Fig. 24. Evolution of instantaneous torque coefficient and relative blade angle of attack as a function of blade azimuthal position for a single rotor blade; TSR=3.30\mathrm{TSR}=3.30.

Fig. 25. Relative velocity contours [m/s] and relative pathlines for rotor blade passing through 92∘92^{\circ} azimuthal position (wind is coming from the left) for optimum tip speed ratio (TSR =2.33=2.33 ).
inside the 4th and 5th Cartesian plane octants, as schematized in Fig. 27, showing the distribution of the instantaneous power coefficient, defined as:
Cp(θ)=P(θ)1/2ρAsV∞2C_{p}(\theta)=\frac{P(\theta)}{1 / 2 \rho A_{s} V_{\infty}^{2}}
as a function of azimuthal position for a single rotor blade and for optimum TSR (2.33). This is probably due to the combination of a great energy extraction exerted by the rotor blade (due to the upwind operation of the rotor blade itself) and a relative high lever arm with respect to the rotor axis. Anyway, the distribution of
aerodynamic forces as a function of rotor blade azimuthal position should be further investigated.

As can be seen, although the average single blade power coefficient is quite low, the instantaneous power coefficient locally (close to 92∘92^{\circ} azimuthal position) reaches the Betz’s limit.

Fig. 27 displays the distribution of the instantaneous power coefficient as a function of time for the three rotor blades and for optimum TSR, showing the contribution of each blade to overall rotor performance. Three peaks of the instantaneous power coefficient can be seen, being three periods contained in a full rotor revolution, each developing in an angular sector of:
Ψ=2π/N=2π/3[rad]\Psi=2 \pi / N=2 \pi / 3[\mathrm{rad}]
Once more, although the average rotor power coefficient is lower ( 0.429 ), the instantaneous power coefficient locally (close to 92∘,212∘92^{\circ}, 212^{\circ} and 336∘336^{\circ} azimuthal position) exceeds the Betz’s limit for almost 40∘40^{\circ} of rotor revolution, reaching a maximum value of 0.664 . As a matter of fact, according to Betz’s law, no turbine could capture more than 59.3 percent of the kinetic energy flux associated to wind [30]. Nevertheless, as pointed out by Jonkman [31], among the assumptions of Rankine-Froude actuator disc theory, the static pressure on the boundary of the streamtube portion enclosing the rotor disc should be equal to unperturbed ambient static pressure. This hypothesis is not verified for the present work, as can be seen from Figs. 28 and 29, representing the contours of the differential pressure coefficient (the atmospheric pressure value was set to 0 Pa ), defined as:
Cp(θ)=p(x,y)1/2ρV∞2C_{p}(\theta)=\frac{p(x, y)}{1 / 2 \rho V_{\infty}^{2}}

Fig. 26. Evolution of instantaneous power coefficient as a function of blade azimuthal position for a single rotor blade. The location of maximum power extraction along blade trajectory, inside the 4th and 5th octants, is highlighted; TSR =2.33=2.33.

Fig. 27. Evolution of instantaneous power coefficient as a function of time for the three rotor blades. The areas in red correspond to the overcoming of the Betz’s limit; TSR =2.33=2.33.

Fig. 28. Contours of differential pressure coefficient cp(θ)c_{p}(\theta) close to rotor disc (wind is coming from the left) for 92∘92^{\circ} azimuthal position (maximum local Cp(θ)C_{p}(\theta) ); TSR =2.33=2.33.
in correspondence to the rotor disc for 92∘92^{\circ} azimuthal position (maximum local Cp(θ)C_{p}(\theta) value, exceeding the Betz limit) and 32∘32^{\circ} (minimum local Cp(θ)C_{p}(\theta) value).

A sudden pressure coefficient drop can be seen, concerning the whole rotor disc and the surrounding flow region, especially for
maximum local Cp(θ)C_{p}(\theta) azimuthal position, where the thrust exerted from the rotor to the flow is almost doubled with respect to the case of minimum local Cp(θ)C_{p}(\theta) azimuthal position, as can be seen from Table 5, reporting rotor instantaneous thrust coefficient, defined as:

Fig. 29. Contours of differential pressure coefficient cp(θ)c_{p}(\theta) close to rotor disc (wind is coming from the left) for 32∘32^{\circ} azimuthal position (minimum local Cp(θ)C_{p}(\theta) ); TSR =2.33=2.33.

Table 5
Rotor thrust coefficient for 92∘92^{\circ} and 32∘32^{\circ} azimuthal positions; TSR =2.33=2.33.

C1(θ)=Fa(θ)1/2ρAsV∞2C_{1}(\theta)=\frac{F_{a}(\theta)}{1 / 2 \rho A_{s} V_{\infty}^{2}}

for 92∘92^{\circ} and 32∘32^{\circ} azimuthal positions.
The sudden pressure coefficient drop in the rotor disc region can be better seen from Fig. 30, showing a 3D representation of pressure coefficient contours for the Rotor-sub-grid zone for 92∘92^{\circ} azimuthal position (maximum local Cp(θ)C_{p}(\theta) ) and TSR =2.33=2.33 (local pressure peaks/drops in correspondence of the three rotor blades should not be considered). This phenomenon could be considered the main responsible for the local increase of the instantaneous power coefficient and should be further investigated.

Fig. 30. 3D representation of pressure coefficient contours for the Rotor-sub-grid zone, 92∘92^{\circ} azimuthal position (maximum local Cp(θ)C_{p}(\theta) ); TSR =2.33=2.33. The black arrow represents the direction of unperturbed wind.

9. Conclusions and future work

In this paper, a numerical model (based on the application to CFD of the aerodynamic principles which are currently applied to BE-M theory for rotor performance prediction) for the evaluation of energy performance and aerodynamic forces acting on a straightbladed vertical-axis Darrieus wind turbine has been presented.

A simplified aerodynamic model (based on the analysis of kinematic and dynamic quantities at discrete and fixed rotor azimuthal positions along blades trajectory) was also presented, allowing the correlation between flow geometric characteristics (such as blade angles of attack) and dynamic quantities (such as rotor torque and blade tangential and normal forces).

The results of a 2-D full campaign of simulations were proposed for a classical NACA 0021 three-bladed rotor. Through the analysis of the distribution of the instantaneous torque coefficients and of the relative angles of attack as a function of azimuthal position for a single rotor blade, flow field characteristics were investigated for several values of tip speed ratio, allowing a comparison between rotor operation at optimum and lower CPC_{\mathrm{P}} values.

The obtained results have shown the reduction of blade relative angles of attack passing from lower to higher TSR values, due to the increasing influence of blade translational speed in the near-blade flow field. It has also been shown that the maximum torque values are generated during the upwind revolution of the turbine and for azimuthal positions where rotor blades are experiencing very high relative angles of attack, even beyond the stall limit. The azimuthal positions of maximum power extraction along blade trajectory have been located inside the 4th and 5th octants, probably due to the combination of a great energy extraction exerted by the rotor blade (due to the upwind operation of the rotor blade itself) and a relative high lever arm with respect to the rotor axis. The distribution of aerodynamic forces as a function of rotor blade azimuthal position should be further investigated.

Finally, it has been shown that, although the average rotor power coefficient is rather lower, the instantaneous power coefficient locally exceeds the Betz’s limit three times per rotor revolution. This phenomenon, probably caused by a sudden pressure coefficient drop concerning the whole rotor disc and the surrounding flow region (thus violating the assumptions of Rankine-Froude actuator disc theory) should be further investigated.

References

[1] Glauert H. Airplane propellers, aerodynamic theory, vol. 4. New York: Dover Publication Inc; 1963. Division L, 169-360.
[2] Templin RJ. Aerodynamic performance of vertical-axis wind machines. ASME; 1975. Paper 75-WA/ENER-1.
[3] Strickland, J.H.: The Darrieus turbine: a performance prediction model using multiple streamtube. 1975. SAND 75-0431.
[4] Read S, Sharpe DJ. An extended multiple streamtube theory for vertical axis wind turbines. Department of M.A.P. Engineering Report. United Kingdom: Kingston Polytechnic, Kingston Upon Times; 1980.
[5] Paraschivoiu I. Double-Multiple streamtube model for Darrieus wind turbines. 2185. NASA Conference Publication; May 1981.
[6] Paraschivoiu I, Delclaux F. Double multiple streamtube model with recent Improvements. J Energy 1983;7(3):250-5.
[7] Paraschivoiu I. Double-Multiple streamtube model for Studying vertical-axis wind turbines. J Propulsion Power July-August 1988;4:4.
[8] Mc Coy H, Loth JL. Up- and downwind rotor half Interference model for VAWT. 2nd Terrestrial Energy Systems Conference. Colorado Springs, CO: AIAA; December 1-3, 1981.
[9] Mc Coy H, Loth JL. Optimization of Darrieus turbines with an upwind and downwind momentum model. J Energy 1983;7(4):313-8.
[10] Lanzafame R, Messina M. Power curve control in micro wind turbine design. Energy February 2010;35(2):556-61.
[11] Kishinami K, Taniguchi H, Suzuki J, Ibano H, Kazunou T, Turuhami M, et al. Study on the aerodynamic characteristics of a horizontal axis wind turbine. Energy August-September 2005;30(11-12):2089-100.
[12] Ceyhan O, Ortakaya Y, Korkem B, Sezer-Uzol N, Tuncer IH. Optimization of horizontal axis wind turbines by using BEM theory and genetic algorithm, 5th Ankara International Aerospace Conference. Ankara: METU; 17-19 August, 2009. AIAC-2009-044.
[13] Zhang, J., Numerical modeling of vertical axis wind turbine (VAWT), Master thesis, Department of Mechanical engineering, Technical University of Denmark, December 20, 2004.
[14] Claessens, M.C., The design and testing of airfoils for application in Small vertical axis wind turbines, Master of Science Thesis, Faculty of Aerospace engineering, Delft University of Technology, November 9, 2006.
[15] Paraschivoiu I. Wind turbine design with emphasis on Darrieus concept. Polytechnic International Press, ISBN 2-553-00931-3; 2002. p. 76.
[16] Sheldal, R. E., Klimas, P. C., Aerodynamic characteristics of Seven symmetrical airfoil sections through 180-degree angle of attack for use in aerodynamic analysis of vertical axis wind turbines. 1980. SAND 80-2114, Unlimited Release, UC-60.
[17] Burton T, Sharpe D, Jenkins N, Bossanyi E. Wind energy handbook. John Wiley & Sons, Ltd; 2001.
[18] Vassberg JC, Gopinath AK, Jameson A. Revisiting the vertical-axis windturbine design using advanced computational fluid dynamics. AIAA Paper; 2005:0047. 43rd AIAA ASM, Reno, NV.
[19] Simao Ferreira CJ, Bijl H, van Bussel G, van Kuik G. Simulating dynamic stall in a 2D VAWT: modeling strategy, verification and validation with particle image velocimetry data, the Science of making torque from wind. J Phys Conf Ser 2007;75.
[20] Simao Ferreira CJ, van Bussel G, Scarano F, van Kuik G. 2D PIV Visualization of dynamic stall on a vertical axis wind turbine. AIAA; 2007.
[21] Kumar V, Paraschivoiu M, Paraschivoiu I. Low Reynolds number vertical axis wind turbine for Mars. Wind Eng June 2010;34(4).
[22] Raciti Castelli M, Benini E. Effect of blade inclination angle on a Darrieus wind turbine. Turbo Expo Technical Conference. Glasgow, Scotland, UK: ASME; June 14-18, 2010. GT2010-23332.
[23] D’Alessandro VD, Montefpare S, Ricci R, Secchiaroli A. Unsteady aerodynamics of a Savonius wind rotor: a new computational approach for the simulation of energy performance. Energy August 2010;35(8):3349-63.
[24] Raciti Castelli M, Pavesi G, Battisti I, Benini E, Ardizzon G. Modeling strategy and numerical validation for a Darrieus vertical axis micro-wind turbine. Vancouver, British Columbia, Canada: ASME 2010 International Mechanical Engineering Congress & Exposition; November 12-18, 2010. IMECE201039548.
[25] Raciti Castelli, M., Analisi numerica delle prestazioni di una micro-turbina eolica ad asse verticale modello Darrieus, PhD Thesis (in Italian), Università di Padova, Italy, 2010, pp. 193-194.
[26] Bradshaw P. Experimental fluid Mechanics. Cambridge University Press; 1964.
[27] Fluent Inc., Fluent User’s Manual. 2006. pp. 52, 54, 59, 71, 143.
[28] Cummings RM, Forsythe JR, Morton SA, Squires KD. Computational challenges in high angle of attack flow prediction. Progr Aerosp Sci 2003;39(5): 369-84.
[29] McMullen M, Jameson A, Alonso JJ. Acceleration of convergence to a periodic steady state in turbomachinery flows, 39th AIAA aerospace sciences meeting & exhibit. Reno, NV: AIAA; January 8-11, 2001. 2001-0152.
[30] Betz A. Das Maximum der theoretisch moglichen Ausnutzung des Windes durch Windmotoren. Z das gesamte Turbinenwesen 1920;26:307-9.
[31] Jonkman JM. Modeling of the UAE wind turbine for refinement of Fast_AD. NREL/TP-500-34755; December 2003. p. 7.