Next Article in Journal
Challenges of Establishing Solar Power Stations in Hungary
Previous Article in Journal
An Analysis on the VSC-HVDC Contribution for the Static Voltage Stability Margin and Effective Short Circuit Ratio Enhancement in Hybrid Multi-Infeed HVDC Systems
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

A Numerical Procedure for Variable-Pitch Law Formulation of Vertical-Axis Wind Turbines

Engineering Department, Università della Campania “L.Vanvitelli”, Via Roma 29, I-81031 Aversa, CE, Italy
*
Author to whom correspondence should be addressed.
Energies 2023, 16(1), 536; https://doi.org/10.3390/en16010536
Submission received: 26 November 2022 / Revised: 24 December 2022 / Accepted: 27 December 2022 / Published: 3 January 2023

Abstract

:
A numerical procedure was developed to determine a variable-pitch law that maximized the performance of a vertical-axis wind turbine (VAWT). The methodology was based on the determination, for each blade, of the angle of attack maximizing the stationary aerodynamic efficiency at prescribed azimuthal positions. The angles of attack were determined by means of a panel method with a low computational effort, and the methodology was implemented in Matlab ® software (version R2021a) allowing us to achieve in real time a variable-pitch law suitable for the turbine geometry. The variable pitch law was validated by considering its effect on the torque of a 2D model of an H-Darrieus turbine. U-RANS analyses were carried out with a K ω S S T model and a sliding-mesh technique was used to prescribe the blade motion around the shaft and pitch motion. Results showed how the variable-pitch law delayed the dynamic stall and improved the aerodynamic performance considerably.

1. Introduction

The reduction of pollutant emissions from fossil combustion and the need for a valid alternative to fossil fuels, due to an increasing demand of coal and natural gas resources, have motivated the European Union (EU) to invest into research on renewable energy. EU has established that renewable energy contribution will be 20% higher than the current amount and at least one-third of this production could be possibly derived from wind energy within the 2020s [1,2]. The technological progress achieved in the renewable energies field has been relevant and, in particular, wind power has become the main renewable energy source in many countries.
Wind turbines represent the main technology able to convert kinetic energy extracted from the wind into electric energy [3,4]. Wind turbines are exposed to various environmental conditions during their operating life (i.e., wind/dust, rain, hail, snow, ice formation, ultraviolet exposure). These operating conditions combined with blade tip speeds result in damage to the leading edge due to leading edge erosion (LEE) effect [5]. There are several studies in the literature about the mechanisms, modelling, and possibilities of preventing the surface erosion of wind turbine blades [6,7].
Wind turbines are classified as horizontal-axis wind turbine (HAWT) and vertical-axis wind turbine (VAWT), according to the spatial orientation of the rotating axis. In the last years, an increasing interest in vertical-axis wind turbine has been experienced because this configuration presents some relevant advantages compared to the horizontal-axis layout. Specifically, it is able to work with an omnidirectional asymptotic velocity without a yaw mechanism that generates power loss. This configuration is also more compact and less noisy [8]. Moreover, with the prospective of urban application, the micro-VAWT can represent a valid alternative to photovoltaic devices.
On the other hand, VAWT is less efficient compared to HAWT, due to the presence of unsteady aerodynamic phenomena, e.g., dynamic stall, deep stall, and blade–wake interaction. The steep variation in time of the angle of attack is among the main causes of the dynamic stall phenomenon, which increases the noise and the fatigue life of the turbine components [9,10]. The delay of the dynamic stall phenomenon has been the subject of study in recent years, and different strategies have been considered in literature. Wang et al. [11] proposed a new technology based on flexible blades. According to Huang’s research [12], blade sections are able to change automatically their camber angle to prevent flow separation. These effects were noticed by Bishop [13] and by Swartz et al. [14] when studying flying squirrels and bats, which are able to fly at high angles of attack by changing the shape of their wings during the flapping flight. Furthermore, Ismail and Vijayaraghavan [15] proposed a combination of semicircular inward dimple and a Gurney flap at the surface of a Darrieus turbine. They observed that the adoption of an optimized combination of a Gurney flap and a dimple produced an increase of the tangential component of the aerodynamic force. Another strategy, adopted to delay the dynamic stall, is an application of turbulators [16] or vortex generators [17,18] on the surface of the blades. These technologies energize the boundary layer making the flow turbulent and reducing the size of the laminar bubble with a consequent rapid reattachment of the flow. Even though turbulators and vortex generators are cheap technologies, they present main cons: (i) the right position on the surface of the blades is difficult to find because it is influenced by the variation of the angle of attack; (ii) the turbulent boundary layer creates extra friction drag [19,20].
Another technology often used to reduce the variation in time of the angle of attack is an active pitch system. Paraschivoiu et al. [21] and Rezaeiha et al. [22] studied the effect of a fixed-pitch angle on the aerodynamic performance of a turbine. However, according to that approach, blade sections do not operate at an optimum angle of attack at each azimuthal position. Indeed, turbine performance can be improved further by utilizing a variable-pitch mechanism. Generally, a linkage mechanism is the most used for this purpose. It is characterized by main links that rotate around the main rotational point, and secondary links that rotate around an eccentric rotational point. This way, the blades are controlled in order to avoid the stall and also maintain a favourable angle of attack. Even if a linkage mechanism is able to improve power generation with respect to a fixed-pitch configuration, the blade sections operate at angles of attack which are not optimum in most azimuthal positions [23,24]. Moreover, the determination of link lengths and eccentric rotational point position requires an optimization process, namely, the best combination of variables whose performance improvement is not known a priori. Therefore, this system is not versatile because it is necessary to execute an optimization at any variation of the operating conditions.
Different methodologies to define a variable-pitch law are presented in literature. Hwang et al. [25] proposed an optimization procedure to determine the optimal pitch angle at specific azimuthal positions. Computational fluid dynamics (CFD) was adopted to estimate the aerodynamic performance coupled with a genetic algorithm optimizer. The optimal variable-pitch law improved the power generation by 25% with respect to the one obtained with the cycloidal motion. Paraschivoiu et al. [21] analysed a two-blade H-Darrieus turbine modelling a variable-pitch law with an analytical function whose coefficients were used as optimization process variables. The aerodynamic model used to determine the performance was a double-multiple streamtube (DMS), while the optimization was carried out by using a genetic algorithm. The results showed an increase in the power generation at different tip-speed ratio (TSR). Both methodologies reported above required a high computational effort and should execute at the variation of the operating conditions. Conversely, Staelens et al. [26] proposed three different methodologies, computationally more efficient, to define a variable-pitch law. The first imposed the local angle of attack of airfoils equal to the stall angle value throughout the rotation of the turbine. Although they obtained an increase of efficiency for a higher wind speed, the variation in time of the angle of attack was physically impossible because it presented discontinuity at azimuthal positions 0 ° and 180 ° . The second methodology proposed to replace the value of the angle of attack with the value of the angle of stall when the former exceeded the latter. Although the aerodynamic performance decreased, the variation of the attack angle was smoother than the previous one but it still presented discontinuities. The last methodology proposed a sinusoidal correction function whose amplitude was equal to the maximum difference between the local angle and the static stall angle. This way, the local angle of attack had a continuous variation in time with respect to the cyclical condition. Although these methodologies allowed an improvement of the aerodynamic performance, the blade worked in an unstable condition because the angles of attack were close to the angle of stall, and aerodynamic performance could drop dramatically due to flow separation.
In the present work, a variable-pitch law that makes each blade section, namely, airfoil, able to operate at the conditions of the maximum steady aerodynamic efficiency, is presented. The proposed methodology is computationally efficient compared to equivalent approaches found in the literature. In fact, for eight azimuthal positions, the static polars of the airfoil are calculated using a potential solver, and the effective angles of attack that maximize steady efficiency are determined. The angle by which each blade section must rotate to operate at an effective angle of attack is the pitch angle. A cubic interpolation of pitch angles produces a variable-pitch law. The automatic procedure was implemented in Matlab software.

2. Aerodynamics of VAWT

A schematic representation of a VAWT is reported in Figure 1. A polar reference system ( R , θ ) having its origin in the centre of the shaft is used to represent the current blade position.
The flow field on a blade is characterized by different velocity types: (i) a tangential velocity component due to rotation ( Ω × R ) where Ω is the rotational velocity and R is the radius vector; (ii) an asymptotic wind velocity V ; (iii) an induced velocity v . According to Figure 1, the sum of those velocities returns the relative velocity of a single blade:
V = V ( Ω × R ) + v
As the turbine rotates, both magnitude and direction of the relative velocity change cyclically. Therefore, the angle between the airfoil chord and the relative velocity, namely, the geometric angle of attack α g , follows the same cyclical variation. The geometric angle of attack can be expressed, neglecting for simplicity the small induced velocity v , as a function of azimuthal position ( θ ) as:
α g = arctan sin θ λ + cos θ
where λ is the tip-speed ratio (TSR) defined by:
λ = Ω R V
The parameter that gives knowledge about the three-dimensionality of the flow is the aspect ratio ( A R ) , defined as:
A R = b c
where b is the blade span while c is the chord of the airfoil. By assuming an A R greater than 7.5, three-dimensional effects can be neglecting [27,28]. The aeronautical convention is adopted for α g , i.e., positive values of α g are associated to anticlockwise rotations of the relative velocity with respect to the airfoil chord, as shown in Figure 2. To ensure an effective operating condition of the turbine, the suction side of the i-th airfoil is oriented toward the rotor centre, as shown in Figure 2.
An inverse proportionality relation links the TSR values with the values of α g [10]. Low values of the TSR relate to steep variations of the α g and vice versa. Steep changes in time of α g , and a high frequency oscillation of airfoils (pitching motion) create an unsteady flow field on turbine blades, known as the dynamic stall condition [29]. As reported in the literature [30], the dynamic stall is characterized by the development and interaction of vortical structures of different sizes generated on the airfoil’s suction surface. Blades experience an increase of the lift force for values of the angle of attack far beyond the static stall angle. This phenomenon occurs when the leading edge vortex (LEV) spans over the whole suction surface side. Consequently, when the airfoil approaches its maximum angle of attack, the LEV detaches from the airfoil causing a sudden lift stall and a consequent rapid increase of drag [31,32].
In Figure 3, a schematization of the aerodynamic forces acting on an airfoil is reported. Tangential ( T ) and normal ( N ) forces are related to lift ( L ) and drag ( D ) forces, respectively, by:
N = L cos α g + D sin α g T = L sin α g D cos α g
The tangential force plays a major role in generating power; indeed, T is the only force which is able to generate a torque at the shaft. As the lift contribution is greater than the drag one, the torque is oriented according to the direction of rotation generating power. Conversely, when the lift stall occurs, the lift contribution is negligible and the torque will be in disagreement with the direction of rotation.

3. Case Study

In the present context, a 2D model of an H-Darrieus rotor was analysed. The chosen airfoil was the NACA 0018 with a chord c of 0.1 m. The blades were positioned at the radius equal to 0.5 m from the shaft with an aspect ratio equal to 10.
In order to capture 3D vortical structures, the computational grid required a spatial refinement in both axial and normal directions. Typically, these computational grids are characterized by more than 80 million elements with a high computational effort. In this work the three-dimensional effects could be neglected in the middle plane allowing a much more affordable 2D analysis. The geometrical and operational properties are summarized in Table 1.

4. Methodology

In this section, we address the methodology used to define a variable-pitch law that improves the power efficiency of the turbine.
We considered eight azimuthal positions, equally spaced Δ θ = 45 ° , in which the operating conditions were computed V , R e , M , α . For each operating condition, airfoil drag polars were computed by means of the panel method X-Foil [33]. The effective angle of attack α e f f was chosen in order to maximize the steady aerodynamic efficiency and generate an anticlockwise torque. The pitch angle α p was defined as the difference between the α e f f and α g :
α p = α e f f α g
A cubic spline was used to fit the values of α p defining the variable-pitch law.

4.1. Steady Polar Airfoil

As reported in Section 2, the tangential component of the aerodynamic force produced a torque on the shaft, which could be synchronous with the turbine motion. Therefore, a positive contribution of the tangential force was obtained by reducing the drag contribution and increasing the lift one, as evidenced in Equation (5). Namely, the airfoil should operate at the effective angles of attack which maximize aerodynamic efficiency, defined as:
E m a x = C l ( α e f f ) C d ( α e f f )
where C l and C d are defined as
C l = L ( α e f f ) 1 2 ρ V 2 S C d = D ( α e f f ) 1 2 ρ V 2 S
where L and D are the lift and drag forces of the airfoil, V is the relative speed, and S is the reference surface defined as the product of c for a unit length. The time history of aerodynamic forces is influenced by the variation in time of α g , due to flow unsteadiness. Therefore, the computational overhead deriving from the computation of α e f f at each azimuthal position, shown in Figure 4, was reduced by assuming the flow steady.
Obviously, the α e f f values represent a conservative solution in the dynamic case because, considering the dynamic polar of the airfoil, the maximum aerodynamic efficiency is obtained at angle of attack values greater than those obtained in the static case. In order to find α e f f , the airfoil polars were computed by using the potential solver X-Foil. Airfoil shape, local Reynolds, and Mach numbers were input variables at each azimuthal position adopting the following relations:
R e = ρ V ( θ ) c μ
M = V ( θ ) a
where V is the relative speed that was computed by means of Equation (1) assuming v = 0 , μ is the dynamic viscosity and a is the speed of sound at sea level. In Table 2 are reported the operating conditions for the azimuthal positions.
The operating conditions at the second, third, and fourth azimuthal positions were similarly obtained at the eighth, seventh and sixth azimuthal positions, respectively. By examining Figure 4, the relative velocity is equal for points along the blade trajectory symmetrically placed with respect to the y-axis. In addition, the airfoil is symmetric, whereby only the airfoil polars at the first four azimuthal positions need to be computed. Furthermore, since the relative speed at the fifth position was close to zero, α e f f was imposed equal to zero to avoid an additional increase of the drag. In Figure 5, the airfoil drag polars computed at the first, second, third, and fourth azimuthal positions are reported.
The condition of maximum aerodynamic efficiency ( E m a x ) was obtained graphically by intersecting the airfoil polar with the tangent curve, as shown in Figure 5. In Table 3, the values of α e f f at the azimuthal positions are reported. The sign of α e f f was positive for azimuthal positions in the range 0 ° θ 180 ° , while it was negative in the range 180 ° θ 360 ° . As a result, tangential forces created an anticlockwise torque, which was compliant with the turbine motion.

4.2. Variable-Pitch Law

The geometrical difference between α e f f and α g at eight azimuthal positions is shown in Figure 6. The α g was computed using Equation (2), and the resulting values are reported in Table 4.
As shown in Figure 6, α p represents the angle by which the airfoil has to rotate in order for α g to be exactly equal to the α e f f at each azimuthal position. In Table 5, the values of α p at eight azimuthal positions are reported.
The pitch law was obtained by fitting a cubic spline interpolator through α p . In particular, the Hermite interpolation implemented in Matlab software [34] was used. In order to ensure a unique solution, and to generate a periodic trajectory, at the spline end points, the same value of the first derivative ( 2 2 ) was imposed. The pitch law is shown in Figure 7b, while a comparison between the variation of α e f f and α g as a function of the azimuthal position is illustrated in Figure 7a.
In Figure 8, the layout of the blade orientation is reported. The methodology described above was implemented in Matlab by generating an automatic procedure. Therefore, the user obtains as output the variable-pitch law in real time by providing as input the chosen airfoil, the radius, and the operating conditions of the turbine. Considering the motion of the turbine without the airfoil pitching, a power loss for each blade during the entire revolution was expected. This behaviour was due to the steep variation of α g , as observed in Figure 7a.

5. CFD Simulation

5.1. Computational Domain

The computational domain used to simulate the motion of the rotor is reported in Figure 9.
The height and length of the domain were set equal to 72.5 c and 225 c , respectively, in order to guarantee a full development of the wake [35]. To simulate both blade and pitch motions, a domain partition into both fixed (Figure 9E,F) and rotating regions was required. The rotating regions were further divided into: (i) three circular regions whose centres were at the quarter of the chord of each airfoil, namely, b l a d e   r e g i o n (Figure 9A–C) and (ii) a ring that included the three circular regions, namely, r o t o r (Figure 9D). That way, all the regions were able to rotate around the shaft and, at same time, the circular regions were able to turn around their centres respecting the pitch law. In the solver Ansys Fluent ® [36], the most accurate method for simulating flows in multiple moving reference frames is the sliding mesh model. An interface condition was imposed at the boundaries of the rotating regions, while no slip condition was imposed at the boundaries of the airfoil. At the upstream and downstream of the fixed region, velocity inlet and pressure outlet conditions were assigned; at the top and bottom of the fixed region, symmetry conditions were adopted. The pitch law was developed by means of an ad hoc user-defined function (UDF) [37] compiled within the solver Fluent.

5.2. User-Defined Function

The user-defined functions are C functions which are interpreted or compiled in commercial software Ansys Fluent, and they can be dynamically loaded in the software tool. These functions allow software customisation to simulate a phenomenon of interest by better approximating reality.
In this work, the D E F I N E _ Z O N E _ M O T I O N UDF was adopted to simulate the airfoil pitch motion during the rotation. Namely, this function allows a rigid translation and rotation of computational domain regions either with respect to the basic reference system or with respect to a generic point. The arguments of the UDF are the motion components such as: rotational axis origin vector, rotational speed, rotation axis direction vector, and translational velocity vector. In this work, the circumferences centred at the c / 4 of each airfoil were chosen as the regions of the computational domain. This choice was taken to simplify the components of the pitch motion to be assigned. In that way, at each time step and for each airfoil, the coordinates of points at f r a c c 4 represented the rotation axis origin vector. These were evaluated by:
x c 4 ( t ) = R cos ( θ t + θ 0 ) y c 4 ( t ) = R sin ( θ t + θ 0 ) z c 4 ( t ) = 0
where θ t and θ 0 are the azimuthal positions of each blade at the current time t and at time t = 0, respectively. The values of rotational speed ω b were evaluated at each time step according to the variable-pitch law. In particular, we discretized the variable-pitch law in accordance with the time step used in the simulation ( Δ t ) . The values of ω b were estimated by means of the finite difference method, as reported in the following equation:
ω b = α p t 1 α p t Δ t
where α p Δ t 1 and α p Δ t are the values of the pitch angle of each blade at the previous time ( t 1 ) and at the current step.
In Figure 10, the flow chart of D E F I N E _ Z O N E _ M O T I O N is reported.

5.2.1. Numerical Method

The numerical simulations presented in this paper were performed by adopting the commercial CFD code Fluent. The unsteady Reynolds averaged Navier–Stokes equations (U-RANS) with the k ω S S T turbulence model was used. The equations were spatially discretized using the finite-volume method and the pressure-based solver was adopted in the solution phase. The coupling between pressure and velocity was achieved by using the coupled algorithm. Convective terms were discretized by adopting the second-order upwind scheme, such as the transport term, momentum, and turbulence. Diffusion terms were central-differenced and were always second-order accurate. Cell gradients were computed by using the least-squares cell-based algorithm [24]. Iterative convergence was obtained when all the variable residuals were less than 10 5 in magnitude [38]; generally, 20 iterations on average were required for convergence. A grid convergence solution was obtained when the mean airfoil torque coefficient between two time periods did not change by more than 0.1% [39]. In addition, a temporal discretization was performed using the first-order implicit Euler method that allowed a larger region of stability through diffusive and convective stability constraints [40].

5.2.2. Meshing Strategy

The discretization of the computational domain is a challenging task, especially at low TSR < < 5 , due to the presence of complex flow regimes [38].
In the present context, a hybrid mesh was adopted to partition the computational domain in order to reduce the number of required elements. As shown in Figure 11, the fixed and rotating domains were discretized using a structured quadrilateral and unstructured triangular elements, respectively.
In addition, layered elements were adopted in the near-wall region of the airfoils in order to correctly evaluate the viscous effects and the detachment of the flow vein. The boundary layer thickness was estimated by means of the following equation:
δ = c R e
where the global Reynolds number was computed from Equation (9) assuming the asymptotic wind velocity V . In order to resolve the viscous sublayer, the height of the first layer guaranteed the condition y + 1 [36,41]. In that case, the y + was set equal to 0.6 to ensure that up to six cells were present in the viscous sublayer. At the sliding interface, both quadrilateral and triangular elements had the same edge size to minimize the interpolation error between the two zones [36].

5.2.3. Spatial Resolution Study

To describe the flow structure that originates from the flow field (vortices, wake, etc.), the time step and the mesh size must be small enough [39,42]. Three computational domains ( G I G I I I ) were analysed with different cell densities, as reported in Table 6, to evaluate the effect of spatial resolution on the numerical solution.
This approach allowed a local mesh refinement where vortical structures were expected. In order to reduce the computational overhead, the turbine motion without the pitch law was analysed. As the optimal computational time step was not a priori known, we set it equal to 1 × 10 4 in order to guarantee a unit Courant number. It was evaluated based on the free stream velocity and the minimum element size among the computational domains, i.e., the most critical condition for the computational stability. In addition, this time step guaranteed both convergence with respect to the convective flux and sliding mesh criteria [36].
The results of the three analyses confirmed that, by using of a coarse computational domain, an underestimation of the mean airfoil torque coefficient C m was obtained, due to a poor spatial discretization on the airfoil surfaces, as shown in Figure 12.
The coefficient C M of the total torque was defined as:
C M = M 1 2 ρ V 2 R A
where M is the aerodynamic torque generated by the blade calculated with respect to the rotation axis, and A is the reference surface defined as the product between diameter and height of the turbine.
The coarse computational domain did not provide a good evaluation of detachment phenomena, returning higher values of the adverse gradient pressure that produced a high drop of the torque coefficient. Conversely, by adopting a medium or fine computational domain, the results were extremely close and thus, utilizing a smaller element size did not increase the accuracy of the predicted results. Since the difference between the results of the medium and fine computational domain was less of 5 % , we used the medium computational domain to analyse the motion thereby reducing the computational effort.

5.2.4. Temporal Resolution Study

The choice of a suitable time step plays a main role to obtain a correct description of vortical structures and wake prediction. The time step was chosen in order to guarantee the lowest computational effort without compromising the accuracy of the results [42]. To assure a time independency of the numerical solution, three different simulations using the medium computational domain at Δ t = 10 3 , Δ t = 10 4 , and Δ t = 10 5 were carried out. By assuming a turbine rotation of 17.4533 rad/s for these cases, corresponding turbine rotation steps equal to Δ θ = 1 ° , Δ θ = 0.1 ° , and Δ θ = 0.01 ° were obtained, respectively.
The variation of the torque coefficients as a function of the azimuthal position θ , for the above defined values of Δ t , is shown in Figure 13.
It can be observed that for time steps Δ θ = 1 ° and Δ θ = 0.1 ° the difference between the results became negligible. Hence, a time step Δ θ = 0.1 ° was adopted for all the simulations. One full rotation took about 24 h utilizing a medium computational domain, Δ t = 10 4 , and using a parallel computing architecture based on a 16-core INTEL Xeon Xeon E5 2.9 GHz. Confidently, in order to predict a single torque coefficient value, six full rotations were performed.

6. Results

In this section, the 2-D numerical simulations are reported. In particular, the effects of the variable-pitch law on the VAWT performance are assessed. In Figure 14, the instantaneous torque coefficient obtained with the pitch law (Case 1) is compared to a fixed pitch (Case 0).
As expected, the behaviour obtained in Case 0 was typical of a highly unsteady aerodynamic field characterized by the existence of vortical structures generated by dynamic stall. Instead, in Case 1, the dynamic stall was delayed as it appears by looking at the behaviour of C M in Figure 14, where the sudden loss of C M is not present. Those different behaviours are highlighted by the contours of the pressure and velocity fields shown in Figure 15 and Figure 16, respectively. In both cases, in the range 0 ° < θ < 10 ° , the static pressure fields around the airfoil showed the flow fully attached to the airfoil (Figure 15a,b). Two main reasons caused this behaviour: (i) the airfoil worked at angles of attack far from dynamic stall conditions; (ii) the aerodynamic field in the upstream phase ( 350 ° < θ < 10 ° ) was not influenced by the wake generated by other airfoils, as confirmed by the velocity field reported in Figure 16a,b. Furthermore, pressure fields obtained at 0 ° < θ < 10 ° , represented in Figure 15a,b, showed an increment of the static pressure difference between the pressure side and suction side of the airfoil in Case 1. This caused the increment of C M compared to Case 0, as shown in Figure 14. In the remaining azimuthal range, the phenomena occurring in Case 0 and Case 1 were different from each other. In Case 0, the values of the angle of attack increased in the range 10 ° < θ < 64 ° .
Moreover, a reversed flow occurred at the trailing edge due to the rapid airfoil motion, as shown in Figure 15 and Figure 16. In the same range, the LEV was already formed, delaying the lift stall to an angle greater than that of the static stall. The LEV increased the depression on the suction side raising the difference in static pressure between the surfaces of airfoil. In that way, the lift force increased compared to the static case. The maximum difference of static pressure was obtained at θ = ˜ 64 ° since, as shown in Figure 15d and Figure 16d, the LEV spanned over the whole suction side of the airfoil. At this value of θ , namely θ * , a peak of instantaneous torque coefficient was obtained. For θ > θ * , the angle of attack achieved the dynamic stall angle value and the LEV detached from the airfoil causing a sudden loss of lift, as shown in Figure 14. Since the variation of α g was steep, there was a considerable delay in the flow reattachment. Only when the airfoil returned to the upstream phase did the flow reattach to the airfoil, as shown in Figure 15e–i and Figure 16e–i. Therefore, in that range, the C M values oscillated around the zero value. On the other hand, in Case 1, the flow was attached to the body until the section blade was in the range 10 ° < θ < 170 ° , as shown in Figure 15b–f. Currently, only in the range 170 ° < θ < 210 ° did separation occur at the trailing edge, and a simultaneous growth of LEV was observed. This behaviour is visible in Figure 15f–h and Figure 16f–h. In Case 1, the flow was quasi-steady, and the airfoils worked at a low angle of attack; at the remaining azimuthal positions, the flow reattached, producing an increase of C M values, as shown in Figure 15i and Figure 16i.
As we expected, the maximum value of C M was obtained in Case 0. This behaviour was due to the methodology used to determine the effective angle. In fact, the static polar used underestimated the maximum efficiency condition compared to the dynamic polar. Nevertheless, the improvement in performance achieved adopting the variable pitch was evident when observing the behaviour of the global torque coefficient of the three blades versus the azimuthal positions reported in Figure 17. For Case 0, the average torque coefficient value C M a was equal to 0.16 compared to a C M a value obtained in Case 1 equal to 0.11 with an improvement of 0.27. The average value of the power coefficient C p a was defined as:
C p a = 0 2 π C M g ω V d θ 2 π = 0 2 π C p d θ 2 π .
where C M g is the torque coefficient of the three blades. In Case 0, C p a was equal to 0.026 . As is known in the literature [43,44], the solidity, the Reynolds number, and the cord length identify a class of turbine with the same performance. In particular, Brusca et al. [43] studied the performances of turbines characterized by airfoil section NACA 0018 at low Reynolds numbers, a high solidity and different TSRs. These geometrical and operational properties identified the same class of turbine model presented in this work. In the literature, there are no experimental studies on this class of turbines at low TSRs because they typically do not produce power but absorb it.
Consequentially, in order to validate the performance obtained with the proposed numerical modelling, we compared the C p a obtained in this work with the one obtained by Brusca et al. at the same TSR. The percent error committed was estimated by means of the following equation:
% e r r o r = C p a B r u s c a C p a C p a B r u s c a 100
and it was equal to 12%.
Globally, the turbine without the pitch law did not produce energy but absorbed it, while, in Case 1, the turbine generated an average power equal to 23 W with C p a equal to 0.27 .

7. Conclusions

Vertical-axis wind turbines represent a suitable technology to convert kinetic energy, extracted from the wind, into electric energy. However, the presence of unstable aerodynamic phenomena such as dynamic stall, deep stall, and blade–wave interaction results in a lower efficiency of such turbines. Several methodologies have been proposed in literature to improve the aerodynamic efficiency of wind turbines, but they turn out to be computationally inefficient. In this study, a new efficient methodology to define a variable-pitch law which improves the performance of an H-Darrieus turbine was implemented using an automatic procedure. It was based on the determination of the angles of attack in eight azimuthal positions where the maximum steady aerodynamics efficiency condition was satisfied. In order to validate this methodology, the effects of a variable pitch on the torques of a 2D model of the H-Darrieus turbine were analysed. In particular, a U-RANS analysis with a K ω S S T turbulent model was carried out. The sliding-mesh technique was used to simulate the pitch motion during the cycle rotation.
A comparison of the results obtained by applying the pitch law as opposed to a no-pitch-law condition highlighted how the pitch law actually delayed the dynamic stall improving the aerodynamic performance considerably. The mean value of the power coefficient obtained by implementing the pitch law results significantly increased. Obviously, the global angle of attack obtained with the proposed methodology represents a suboptimal solution in the dynamic case because the dynamic maximum aerodynamic efficiency is obtained at values of the angle of attack greater than the static case. However, this methodology is more versatile than those already present in the literature because it does not require optimization and could be implemented in hardware to modify the pitch law in real time, following different asymptotic wind conditions.

Author Contributions

Conceptualization, C.R. and L.I.; methodology and software, C.R. and D.D.S.; validation, C.R. and L.I.; formal analysis, C.R. and D.D.S.; resources, A.V.; data curation, C.R., D.D.S. and G.P.; writing, C.R. and L.I.; writing and editing, C.R. and L.I.; review, G.P. and A.V.; supervision, A.V. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Wilkes, J.; Moccia, J.; Arapogianni, A.; Dragan, M.; Plytas, N.; Genachte, A.; Guillet, J.; Wilczek, P. The European Offshore Wind Industry Key 2011 Trends and Statistics; Report of European Wind Energy Association; European Wind Energy Association: Brussels, Belgium, 2012. [Google Scholar]
  2. Pryor, S.; Barthelmie, R. Climate change impacts on wind energy: A review. Renew. Sustain. Energy Rev. 2010, 14, 430–437. [Google Scholar] [CrossRef]
  3. Stathopoulos, T.; Alrawashdeh, H.; Al-Quraan, A.; Blocken, B.; Dilimulati, A.; Paraschivoiu, M.; Pilay, P. Urban wind energy: Some views on potential and challenges. J. Wind. Eng. Ind. Aerodyn. 2018, 179, 146–157. [Google Scholar] [CrossRef]
  4. Al-Quraan, A.; Stathopoulos, T.; Pillay, P. Comparison of wind tunnel and on site measurements for urban wind energy estimation of potential yield. J. Wind. Eng. Ind. Aerodyn. 2016, 158, 1–10. [Google Scholar] [CrossRef] [Green Version]
  5. Keegan, M.H.; Nash, D.H.; Stack, M.M. On erosion issues associated with the leading edge of wind turbine blades. J. Phys. D Appl. Phys. 2013, 46, 383001. [Google Scholar] [CrossRef] [Green Version]
  6. Mishnaevsky, L.; Hasager, C.B.; Bak, C.; Tilg, A.M.; Bech, J.I.; Doagou Rad, S.; Fæster, S. Leading edge erosion of wind turbine blades: Understanding, prevention and protection. Renew. Energy 2021, 169, 953–969. [Google Scholar] [CrossRef]
  7. Carraro, M.; Vanna, F.D.; Zweiri, F.; Benini, E.; Heidari, A.; Hadavinia, H. CFD Modeling of Wind Turbine Blades with Eroded Leading Edge. Fluids 2022, 7, 302. [Google Scholar] [CrossRef]
  8. Eriksson, S.; Bernhoff, H.; Leijon, M. Evaluation of different turbine concepts for wind power. Renew. Sustain. Energy Rev. 2008, 12, 1419–1434. [Google Scholar] [CrossRef]
  9. Ferreira, C.S.; Van Kuik, G.; Van Bussel, G.; Scarano, F. Visualization by PIV of dynamic stall on a vertical axis wind turbine. Exp. Fluids 2009, 46, 97–108. [Google Scholar] [CrossRef] [Green Version]
  10. Bos, R. Self-Starting of a Small Urban Darrieus Rotor. Strategies to Boost Performance in Low-Reynolds-Number Flows. Master’s Thesis, TU Delft University, Delft, The Netherlands, 2012. [Google Scholar]
  11. 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. Manag. 2016, 108, 275–286. [Google Scholar] [CrossRef]
  12. Huang, D.; Wu, G. Preliminary study on the aerodynamic characteristics of an adaptive reconfigurable airfoil. Aerosp. Sci. Technol. 2013, 27, 44–48. [Google Scholar] [CrossRef]
  13. Bishop, K.L. The relationship between 3-D kinematics and gliding performance in the southern flying squirrel, Glaucomys volans. J. Exp. Biol. 2006, 209, 689–701. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  14. Swartz, S.; Iriarte-Diaz, J.; Riskin, D.; Tian, X.; Song, A.; Breuer, K. Wing structure and the aerodynamic basis of flight in bats. In Proceedings of the 45th AIAA Aerospace Sciences Meeting and Exhibit, Reno, NV, USA, 8–11 January 2007. [Google Scholar]
  15. Ismail, M.; Vijayaraghavan, K. Investigation of flow over an oscillating airfoil for wind turbine application. In Proceedings of the Canadian Society for Mechanical Engineering (CSME) International Congress, Toronto, ON, Canada, 1–4 June 2014. [Google Scholar]
  16. Gopalarathnam, A.; Broughton, B.A.; McGranahan, B.D.; Selig, M.S. Design of low Reynolds number airfoils with trips. J. Aircr. 2003, 40, 768–775. [Google Scholar] [CrossRef] [Green Version]
  17. Mai, H.; Dietz, G.; Geißler, W.; Richter, K.; Bosbach, J.; Richard, H.; de Groot, K. Dynamic stall control by leading edge vortex generators. J. Am. Helicopter Soc. 2008, 53, 26–36. [Google Scholar] [CrossRef]
  18. Paraschivoiu, I. Wind Turbine Design with Emphasis on Darrieus Concept; Polytechnic International Press: Montreal, QC, Canada, 2002. [Google Scholar]
  19. Ahmed, M.R. Blade sections for wind turbine and tidal current turbine applications–current status and future challenges. Int. J. Energy Res. 2012, 36, 829–844. [Google Scholar] [CrossRef]
  20. Howell, R.; Qin, N.; Edwards, J.; Durrani, N. Wind tunnel and numerical study of a small vertical axis wind turbine. Renew. Energy 2010, 35, 412–422. [Google Scholar] [CrossRef] [Green Version]
  21. Paraschivoiu, I.; Trifu, O.; Saeed, F. H-Darrieus wind turbine with blade pitch control. Int. J. Rotating Mach. 2009, 2009, 505343. [Google Scholar] [CrossRef] [Green Version]
  22. Rezaeiha, A.; Kalkman, I.; Blocken, B. Effect of pitch angle on power performance and aerodynamics of a vertical axis wind turbine. Appl. Energy 2017, 197, 132–150. [Google Scholar] [CrossRef] [Green Version]
  23. Kiwata, T.; Yamada, T.; Kita, T.; Takata, S.; Komatsu, N.; Kimura, S. Performance of a vertical axis wind turbine with variable-pitch straight blades utilizing a linkage mechanism. J. Environ. Eng. 2010, 5, 213–225. [Google Scholar] [CrossRef] [Green Version]
  24. 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–123. [Google Scholar] [CrossRef]
  25. Hwang, I.S.; Lee, Y.H.; Kim, S.J. Optimization of cycloidal water turbine and the performance improvement by individual blade control. Appl. Energy 2009, 86, 1532–1540. [Google Scholar] [CrossRef]
  26. Staelens, Y.; Saeed, F.; Paraschivoiu, I. A straight-bladed variable-pitch VAWT concept for improved power generation. In ASME 2003 Wind Energy Symposium; American Society of Mechanical Engineers: New York, NY, USA, 2003; pp. 146–154. [Google Scholar]
  27. Li, Y.; Calisal, S.M. Three-dimensional effects and arm effects on modeling a vertical axis tidal current turbine. Renew. Energy 2010, 35, 2325–2334. [Google Scholar] [CrossRef]
  28. Kirke, B.K. Evaluation of Self-Starting Vertical Axis Wind Turbines for Stand-Alone Applications. Ph.D. Thesis, Griffith University Australia, South East Queensland, Australia, 1998. [Google Scholar]
  29. Amet, E.; MaÃŽtre, T.; Pellone, C.; Achard, J.L. 2D numerical simulations of blade-vortex interaction in a Darrieus turbine. J. Fluids Eng. 2009, 131, 111103–111115. [Google Scholar] [CrossRef]
  30. Gordon, L.J. Principles of Helicopter Aerodynamics; Cambridge University Press: Cambridge, UK, 2016. [Google Scholar]
  31. Wang, S.; Ingham, D.B.; Ma, L.; Pourkashanian, M.; Tao, Z. Numerical investigations on dynamic stall of low Reynolds number flow around oscillating airfoils. Comput. Fluids 2010, 39, 1529–1541. [Google Scholar] [CrossRef]
  32. Hand, B.; Kelly, G.; Cashman, A. Numerical simulation of a vertical axis wind turbine airfoil experiencing dynamic stall at high Reynolds numbers. Comput. Fluids 2017, 149, 12–30. [Google Scholar] [CrossRef]
  33. Drela, M.; Youngren, H. XFOIL 6.94 User Guide. 2001. Available online: http://web.mit.edu/drela/Public/web/xfoil/ (accessed on 24 December 2022).
  34. MatLab, M. The Language of Technical Computing; The MathWorks, Inc.: Natick, MA, USA, 2012; Available online: http://www.mathworks.com (accessed on 24 December 2022).
  35. Rezaeiha, A.; Kalkman, I.; Blocken, B. CFD simulation of a vertical axis wind turbine operating at a moderate tip speed ratio: Guidelines for minimum domain size and azimuthal increment. Renew. Energy 2017, 107, 373–385. [Google Scholar] [CrossRef] [Green Version]
  36. Fluent, A. 12.0 Theory Guide; Ansys Inc.: Canonsburg, PA, USA, 2009; Volume 5. [Google Scholar]
  37. Fluent, A. 15.0 UDF Manual; Ansys Inc.: Canonsburg, PA, USA, 2013. [Google Scholar]
  38. 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–435. [Google Scholar] [CrossRef]
  39. Balduzzi, F.; Bianchini, A.; Ferrara, G.; Ferrari, L. Dimensionless numbers for the assessment of mesh and timestep requirements in CFD simulations of Darrieus wind turbines. Energy 2016, 97, 246–261. [Google Scholar] [CrossRef]
  40. Ferziger, J.H.; Peric, M. Computational Methods for Fluid Dynamics; Springer Science & Business Media: Berlin, Germany, 2012. [Google Scholar]
  41. Andersson, B.; Andersson, R.; Håkansson, L.; Mortensen, M.; Sudiyo, R.; Van Wachem, B. Computational Fluid Dynamics for Engineers; Cambridge University Press: Cammbridge, UK, 2011. [Google Scholar]
  42. Almohammadi, K.; Ingham, D.; Ma, L.; Pourkashan, M. Computational fluid dynamics (CFD) mesh independency techniques for a straight blade vertical axis wind turbine. Energy 2013, 58, 483–493. [Google Scholar] [CrossRef]
  43. Brusca, S.; Lanzafame, R.; Messina, M. Design of a vertical-axis wind turbine: How the aspect ratio affects the turbine’s performance. Int. J. Energy Environ. Eng. 2014, 5, 333–340. [Google Scholar] [CrossRef] [Green Version]
  44. Eboibi, O.; Danao, L.A.M.; Howell, R.J. Experimental investigation of the influence of solidity on the performance and flow field aerodynamics of vertical axis wind turbines at low Reynolds numbers. Renew. Energy 2016, 92, 474–483. [Google Scholar] [CrossRef]
Figure 1. Two-dimensional view of a three-bladed rotor, showing the velocity triangles.
Figure 1. Two-dimensional view of a three-bladed rotor, showing the velocity triangles.
Energies 16 00536 g001
Figure 2. Positive angle of attack.
Figure 2. Positive angle of attack.
Energies 16 00536 g002
Figure 3. Aerodynamic forces on a VAWT blade cross section.
Figure 3. Aerodynamic forces on a VAWT blade cross section.
Energies 16 00536 g003
Figure 4. Velocity triangles and blade orientation at eight azimuthal positions.
Figure 4. Velocity triangles and blade orientation at eight azimuthal positions.
Energies 16 00536 g004
Figure 5. Airfoil polars at azimuthal positions: (a) first azimuthal position; (b) second azimuthal position; (c) third azimuthal position; (d) fourth azimuthal position.
Figure 5. Airfoil polars at azimuthal positions: (a) first azimuthal position; (b) second azimuthal position; (c) third azimuthal position; (d) fourth azimuthal position.
Energies 16 00536 g005
Figure 6. Pitching motion of the airfoil.
Figure 6. Pitching motion of the airfoil.
Energies 16 00536 g006
Figure 7. (a) Variation of α g and α e f f vs. azimuthal positions; (b) pitch law.
Figure 7. (a) Variation of α g and α e f f vs. azimuthal positions; (b) pitch law.
Energies 16 00536 g007
Figure 8. (left) Layout of the blade orientations during a period of revolution; (right) layout of the blade orientations at the same positions as in Figure 4.
Figure 8. (left) Layout of the blade orientations during a period of revolution; (right) layout of the blade orientations at the same positions as in Figure 4.
Energies 16 00536 g008
Figure 9. Layout of the computational domain (not at scale). Rotating regions: (AD); fixed regions: (E,F); dot: interfaces between domains; square: details of computational domain.
Figure 9. Layout of the computational domain (not at scale). Rotating regions: (AD); fixed regions: (E,F); dot: interfaces between domains; square: details of computational domain.
Energies 16 00536 g009
Figure 10. Flow chart of D E F I N E _ Z O N E _ M O T I O N (UDF).
Figure 10. Flow chart of D E F I N E _ Z O N E _ M O T I O N (UDF).
Energies 16 00536 g010
Figure 11. Details of the discretized computational domain: (a) interface rotor and blade region (I in Figure 9); (b) boundary layer (II in Figure 9); (c) interface rotor and fixed domain (III in Figure 9); (d) interface rotor and fixed domain (IV in Figure 9).
Figure 11. Details of the discretized computational domain: (a) interface rotor and blade region (I in Figure 9); (b) boundary layer (II in Figure 9); (c) interface rotor and fixed domain (III in Figure 9); (d) interface rotor and fixed domain (IV in Figure 9).
Energies 16 00536 g011
Figure 12. Spatial independence.
Figure 12. Spatial independence.
Energies 16 00536 g012
Figure 13. Temporal independence.
Figure 13. Temporal independence.
Energies 16 00536 g013
Figure 14. Instantaneous torque coefficient versus azimuthal position: case without pitch law (Case 0) and with pitch law (Case 1).
Figure 14. Instantaneous torque coefficient versus azimuthal position: case without pitch law (Case 0) and with pitch law (Case 1).
Energies 16 00536 g014
Figure 15. Pressure field around the no-pitch and variable-pitch blades at different azimuthal positions.
Figure 15. Pressure field around the no-pitch and variable-pitch blades at different azimuthal positions.
Energies 16 00536 g015
Figure 16. Relative velocity field around the no-pitch and variable-pitch blades at different azimuthal positions.
Figure 16. Relative velocity field around the no-pitch and variable-pitch blades at different azimuthal positions.
Energies 16 00536 g016
Figure 17. Torque coefficient versus azimuthal position: case without pitch law ( C M m e d i u m = 0.16 ) and with pitch law ( C M m e d i u m = + 0.11 ).
Figure 17. Torque coefficient versus azimuthal position: case without pitch law ( C M m e d i u m = 0.16 ) and with pitch law ( C M m e d i u m = + 0.11 ).
Energies 16 00536 g017
Table 1. Geometrical and operational properties of the VAWT under study.
Table 1. Geometrical and operational properties of the VAWT under study.
ParameterValue
Number of blades, n3
Radius, R (m)0.5
Span length, b (m)1
Swept area, A (m2)1
Chord/radius ratio, c R 0.1
AirfoilNACA 0018      
Solidity, σ ( N c R ) 0.6
Tip-speed ratio (TSR), λ 1.25
Asymptotic wind velocity, V   ( m s ) 7
Global Reynolds number ( ρ V c ν ) 1 × 10 5
Rotational speed, Ω   ( rad s ) 17.4533
Table 2. Operating conditions at the eight azimuthal positions.
Table 2. Operating conditions at the eight azimuthal positions.
θ V R e M   
First position 0 ° 15.73 1.1 × 10 5 0.046
Second position 45 ° 14.54 1.0 × 10 5 0.04
Third position 90 ° 11.19 7.6 × 10 4 0.03
Fourth position 135 ° 6.23 4.2 × 10 4 0.02
Fifth position 180 ° 1.7 1.7 × 10 4 0.005
Sixth position 225 ° 6.23 4.2 × 10 4 0.02
Seventh position 270 ° 11.19 7.6 × 10 4 0.03
Eighth position 315 ° 14.54 1.1 × 10 5 0.04
Table 3. Values of α e f f at the eight azimuthal positions.
Table 3. Values of α e f f at the eight azimuthal positions.
θ α e f f
First position 0 ° 7 °
Second position 45 ° 7 °
Third position 90 ° 6.5 °
Fourth position 135 ° 4.5 °
Fifth position 180 ° 0 °
Sixth position 225 ° 4.5 °
Seventh position 270 ° 6.5 °
Eighth position 315 ° 7 °
Table 4. Values of α g at the eight azimuthal positions.
Table 4. Values of α g at the eight azimuthal positions.
θ α g
First position 0 ° 0 °
Second position 45 ° 19.89 °
Third position 90 ° 38.73 °
Fourth position 135 ° 52.6 °
Fifth position 180 ° 0 °
Sixth position 225 ° 52.65 °
Seventh position 270 ° 38.73 °
Eighth position 315 ° 19.89 °
Table 5. Values of α p at the eight azimuthal positions.
Table 5. Values of α p at the eight azimuthal positions.
θ α p
First position 0 ° 7 °
Second position 45 ° 12.90 °  
Third position 90 ° 32.23 °
Fourth position 135 ° 48.15 °
Fifth position 180 ° 0 °
Sixth position 225 ° 48.15 °
Seventh position 270 ° 32.23 °
Eighth position 315 ° 12.90 °
Table 6. Details of the computational domains.
Table 6. Details of the computational domains.
Node Numbers G I G II G III
Airfoil208416416
Layers404040
Inlet206412824
Outlet206412824
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Rainone, C.; De Siero, D.; Iuspa, L.; Viviani, A.; Pezzella, G. A Numerical Procedure for Variable-Pitch Law Formulation of Vertical-Axis Wind Turbines. Energies 2023, 16, 536. https://doi.org/10.3390/en16010536

AMA Style

Rainone C, De Siero D, Iuspa L, Viviani A, Pezzella G. A Numerical Procedure for Variable-Pitch Law Formulation of Vertical-Axis Wind Turbines. Energies. 2023; 16(1):536. https://doi.org/10.3390/en16010536

Chicago/Turabian Style

Rainone, Cinzia, Danilo De Siero, Luigi Iuspa, Antonio Viviani, and Giuseppe Pezzella. 2023. "A Numerical Procedure for Variable-Pitch Law Formulation of Vertical-Axis Wind Turbines" Energies 16, no. 1: 536. https://doi.org/10.3390/en16010536

APA Style

Rainone, C., De Siero, D., Iuspa, L., Viviani, A., & Pezzella, G. (2023). A Numerical Procedure for Variable-Pitch Law Formulation of Vertical-Axis Wind Turbines. Energies, 16(1), 536. https://doi.org/10.3390/en16010536

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop