- Email: [email protected]
Analysis of a free vortex wake model for the study of the rotor and near wake flow of a vertical axis wind turbine
Renewable Energy 87 (2016) 552e563
Contents lists available at ScienceDirect
Renewable Energy journal homepage: www.elsevier.com/locate/renene
Analysis of a free vortex wake model for the study of the rotor and near wake flow of a vertical axis wind turbine ~o Ferreira, G.J.W. van Bussel G. Tescione*, C.J. Sima Delft University Wind Energy Research Institute, TUDelft, Kluyverweg 1, 2629 HS, Delft, The Netherlands
a r t i c l e i n f o
a b s t r a c t
Article history: Received 16 February 2015 Accepted 1 October 2015 Available online xxx
The 3D dynamics of the near wake of a vertical axis wind turbine is simulated by a numerical model based on a free vortex wake/panel method. The model is first verified with respect to discretization sensitivity, stability and robustness. Results show little sensitivity to the discretization, requiring a medium spatial and temporal grid to reach convergence. Stability issues are addressed, identifying dependency of the simulation parameters and validity of the results. The in-rotor flow shows little or no effect of instabilities and the model proves to be an effective analysis tool for the flow at the blade level. Wake dynamics, especially the tip vortex evolution, is more subject to stability issues which have to be taken into consideration. The numerical results are then validated with experimental data to assess accuracy. The model is able to capture all the important dynamics in the near wake of a vertical axis wind turbine with good quantitative evaluation of the flowfield. Discrepancies in the comparison are discussed and limitations in the use of the model are presented. © 2015 Elsevier Ltd. All rights reserved.
Keywords: Vertical axis wind turbine Vortex element methods Panel method Numerical aerodynamics Wake aerodynamics Vorticity dynamics Verification and validation
1. Introduction Future trends of wind energy industry show the development of large turbines in offshore wind farms. As reported by the European Wind Energy Association [23], the latest generation of offshore wind turbines has rotor diameters up to 170 m and rated capacity up to 7 MW, while 10e20 MW machines are currently in the development phase. Offshore wind farm development trends towards bigger farms built further from te coast and in deeper waters. Going beyond the current limit of 50 m of water depth for fixed substructure concepts will allow to expand the market in deeper seas, as the remaining 66% of the North Sea, the Mediterranean and the Atlantic basins in Europe, the coast of North America, China or Japan. Within these waters, floating support structures are shown to be more economical [24]. The development of large floating wind farms of multi-MW turbines has made Vertical Axis Wind Turbines (VAWTs) recently gain new popularity in the wind energy community, as their use
Acronyms: BWI, BladeeWake Interaction; CPAN, chordwise panels; GPU, Graphic Processing Unit; H-VAWT, H-shaped VAWT; HAWT, Horizontal Axis Wind Turbine; ROT, simulated rotations; SPAN, spanwise panels; TSR, Tip Speed Ratio; VAWT, Vertical Axis Wind Turbine. * Corresponding author. E-mail address: [email protected] (G. Tescione). http://dx.doi.org/10.1016/j.renene.2015.10.002 0960-1481/© 2015 Elsevier Ltd. All rights reserved.
can potentially alleviate some of the new challenges [2,12,13]. VAWTs are mechanically simpler machines than Horizontal Axis Wind Turbines (HAWTs), having no yaw, and often no pitch, mechanisms. In an offshore environment, where operation and maintenance have a high weight in the Cost of Energy, the presence of fewer moving parts has a great value. Secondly, VAWTs show higher potentials of scalability, overcoming some of the limitations of HAWTs in this aspect. Moreover the possibility to place the generator at sea level or even under water helps stability and would decrease the size and cost of the floating support structure. Different research programmes in Europe and North America [20,22], are exploring the potential of floating large VAWTs. This renewed interest finds however a fragmented research field. After its peak during the 80's, VAWT research was discontinued for almost 15 years, during the 90's and early years of 2000. The commercial success of the propeller type horizontal axis turbines drove most of the research effort of the wind energy community, with few isolated researchers focusing on the vertical axis turbines, and their use limited to small turbines for urban environment. VAWTs missed the extensive aerodynamic research characterizing HAWTs, resulting in a lack of a solid phenomenological background and an underdevelopment of proper models to assist analysis and design [3,16]. The aerodynamics of vertical axis turbines differs from that of the horizontal axis devices in several fundamental aspects. The
G. Tescione et al. / Renewable Energy 87 (2016) 552e563
List of symbols
az0(q) az0j90 m s f f∞ Vy ! u GTE g gS gT
Dq q c D H M
perceived angle of attack of the middle section of the blade az0(q) at the most upwind position (q ¼ 90 ) doublet singularity strength source singularity strength velocity potential undisturbed velocity potential rotational velocity vorticity vector circulation on the trailing edge segment vortex lament strength shed vorticity (wake filament) trailing vorticity (wake filament) azimuth step azimuth angle blade chord rotor diameter rotor height Mach number
N ! n R Rec std(V) ! s t ! u ! ui V∞ VAVG V (q) Vx Vy Vz XST x y z
553
number of bodies in N-body problem normal vector (to body/wake surface) rotor radius chord-based Reynolds number standard deviation of velocity trailing edge segment vector time velocity vector (general form) induced velocity vector (general form) free-stream velocity averaged absolute velocity instantaneous absolute velocity stream-wise velocity component cross-stream velocity component vertical velocity component limit of stability (in D) for the wake stream-wise extension stream-wise spatial coordinate cross-stream spatial coordinate vertical spatial coordinate
3. Numerical model
energy conversion process is based on the variation of the blade's circulation along its rotation, resulting in a different creation of the wake. Unsteadiness is an inherent part of the turbine operation and its wake is characterized by regions of concentrated vorticity interacting and evolving in coherent structures [15,21]. Accurate wake dynamics prediction is fundamental in VAWT design process as part of the rotor operates in its own wake. Moreover in wind farm design the far wake evolution is of high interest. Low-fidelity engineering models developed for HAWT have been adapted for VAWT design and analysis without considering the inherent differences, resulting in poor predictions [15e17]. Moreover they lack the details of the wake development. Higherorder, grid-based numerical methods tend to underestimate the role and the evolution of the vortical structures due to their inherent numerical dissipation. Vortex element methods have been used in VAWT research for their ability to retain the strength of the vortical structures and to give details of the wake evolution [1,4,11,19,25,26].
The model used for this numerical investigation is 3D, unsteady and based on a panel method coupled with a free vortex wake. Aerodynamic panel methods have been known since the 1960s and extensively used in the aeronautic field as analysis and design tools.1 The model is based on the assumption of an incompressible and potential (inviscid, irrotational) flow, with the exception of vorticity singularities introduced in the wake elements and on the body. Given the relative low Mach number (M < 0.3) and high Reynolds numbers (O(Rec)) experienced by VAWTs such assumptions are considered to be reasonable everywhere in the flow outside the boundary layer. The potential flow assumption allows for the use of a Boundary Integral Equation for the velocity potential (f, the velocity being ! u ¼ Vf) expressed in one of its general forms as:
2. Approach
f¼
The present paper investigates the use of an unsteady free vortex-wake model developed at TUDelft [5] for the analysis of the aerodynamics of the wake of a VAWT. First, the model is introduced by the mathematical formulation, the numerical implementation used and the simulated VAWT operation. Then, a verification of the model is performed addressing convergence and stability of the solution. Geometrical discretizations of the paneled body and of the vortex wake and temporal discretizations are related to the perceived angle of attack of the blade section at the symmetry plane, to assess sensitivity of the results. The effect of the initial transient solution originated by the impulsive start of the simulation is also determined and guidelines on the duration of such transient are given. Stability issues concerning the evolution in the wake downwind the turbine are also addressed to evaluate the applicability of the model. Finally, the results of the verified model are validated against the experimental measurements presented and discussed by the present authors in Ref. [21].
3.1. Mathematical formulation
1 4p
Z body
1 v 1 1 sþ m dS þ r vn r 4p
Z wake
v 1 mdS þ f∞ vn r (1)
where s and m are the singularities on the body and wake surface (S) and r is the distance between the singularity and the location of the potential evaluation. The present model uses a coupled formulation of sources ðsÞ and doublets ðmÞ on the actual body geometry for the singularity mix [10]. Here the doublet strength is taken as independent on each panel and equals the exterior surface value of perturbation potential while the source strength equals its ! ! ! normal derivative ( n ,Vf ¼ n , u ) to fulfill the zero-normalvelocity boundary condition. To solve for the non-uniqueness of the solution the Kutta condition, that the pressure difference vanishes at the trailing edge, is used, which drives the entire solution
1 For a review on the history of the methods and overview of the different formulations the reader is referred to Hess [6,25], and the book from Katz & Plotkin [9].
554
G. Tescione et al. / Renewable Energy 87 (2016) 552e563
by determining the circulation distribution. The wake is modeled as a curved lattice of straight vortex filaments. The wake contribution to the velocity potential in Equation (1) (the second integral) is shown as a doublets influence; the actual model shifts from doublets to vortex filaments ðgÞ after few time steps. This is done as the use of a vortex lattice allows for a better handling of deformations of the wake geometry. The equivalence of the doublets and vortex distribution is well known [7]. Given the incompressibility assumption the time dependency of the solution follows straightforwardly from the time-dependent boundary conditions, hence the unsteady nature of the model. This is of primary importance given the inherent unsteadiness of VAWT aerodynamics [15,21]. The use of vorticity to describe the flow brings important advantages when associated to a Lagrangian approach and modifies the governing equation to a simplified form. The further hypothesis of no viscous diffusion and incompressible flow lead to the following formulation for the transport equation of the vorticity:
D! u ! ð! u $VÞ u ¼ 0 Dt
(2)
However the non-linearity of the original NaviereStokes equations is maintained in the vorticity formulation as the resulting equation describes a non-linear transport phenomenon, since the ! transport velocity ð u Þ depends itself on the transport quantity (the vorticity ! u ). The vorticity field induces a velocity field which is responsible for the evolution of the Lagrangian particles. The recovery of the velocity from the vorticity field is obtained with the integral approach. The Poisson equation for the velocity ! (V2 u ¼ V ! u ) is solved by means of the Green's function for the Laplacian. This leads to the Biot-Savart law for the velocity induc! tion ð u i Þ by a vorticity field ð! u Þ:
1 ! ! u i ð x ; tÞ ¼ 4p
Z !! ! ! u ð y ; tÞ ð x y Þ ! dy ! !3 x y D
(3)
which is the basis for the wake-to-wake and wake-to-body induction calculation. The advantage of a vorticity formulation rather than a velocity formulation is that vorticity has a compact support, so the integration is limited to the domain D where vorticity is nonzero while the velocity field can be evaluated on a larger domain.
3.2. Numerical implementation The numerical model is implemented in hybrid MATLAB®/ CUDA®-C code for an easier use with an Object Oriented Programming approach to increase flexibility. The code is, indeed,
capable of handling multi-body objects in 2D or 3D of generic geometries and different arbitrary (and time-dependent) motions [5]. The rotor blades are discretized with flat panels of constant source and doublet strengths. The geometric discretization is refined towards the leading and trailing edges and towards the tips of the blades to account for curvature and higher pressure (velocity) gradients. No deformation of the body is accounted for, so the blades are subject to a rigid rotation around the vertical axis. This allows to build in a pre-processing phase the Aerodynamic Influence Coefficients of each blade panel to each other, and to use matrix multiplication during the simulation to compute the bodyto-body induction. The wake is modeled as a curved lattice of straight vortex filaments (see Fig. 1). The strength of the wake elements released each time step is derived by the linearized Kutta condition at the blade ! trailing edge segments (along the coordinate s ). Spatial derivatives of the circulation between adjacent trailing edge segments ( DG !) Ds determine the trailing vorticity ðgtÞ of the wake filaments perpendicular to the trailing edge. Temporal derivatives (DG Dt ) determine the shed vorticity ðgSÞ of the wake filaments parallel to the trailing edge. Wake elements are emitted only by the trailing edges, no wake elements are emitted in the spanwise direction from the tip of the blades. The new wake marker position is defined by calculating the local kinematic velocity at the trailing edge. A free-wake analysis is used where the new position of each wake marker is an integral part of the problem. This is essential in the cases of VAWT as the rotor operates partly in its own wake. For each time step the induction field of all the wake filaments and of the body panels is solved on the location of the wake markers (the filaments' endpoints) by a direct summation. The filaments have a viscous core with a Rankine-type behavior to ensure desingularization of Equation (3). No viscous diffusion of the filament vorticity is currently included in the model. As the endpoints of the filaments are used as markers, vortex stretching and tilting is auto! matically considered (the ð! u ,VÞ u term in Equation (2)). Because of its completely Lagrangian formulation the model requires no grid creation and avoids numerical dissipation of vorticity. An explicit second order Adam-Bashforth method is used for the time marching scheme. No wake cut-off is used: all wake elements are retained for the whole simulation. The resulting N-body problem (with an O(N2) computational cost and an OðNÞ memory access) is very suitable for computational parallelization, especially with the use of Graphic Processing Unit (GPU) -based systems [18]. The most computational expensive functions, the wake-to-wake and waketo-body induction calculations, have been implemented in CUDA®-C and run parallel on a NVIDIA® Quadro 4000 GPU (256 cores, 243GFlops in double precision and 89.6 GB/s memory
Fig. 1. Schematics of the vortex lattice wake discretization and the different contributions to the wake vorticity.
G. Tescione et al. / Renewable Energy 87 (2016) 552e563
Fig. 2. 3D rendering of the VAWT model used in the experiments (dimensions in mm).
bandwidth). This brings to a speed up of the calculations up to 120 with respect to a pure MATLAB® execution for a mid-refined simulation. This allows for higher refined and/or long-time simulations. 3.3. VAWT geometry and operation A model H-shaped VAWT (H-VAWT) of a previous experimental campaign [21], is used for the current calculations. It consists of 2 straight blades of 1 m height (H) with a NACA0018 airfoil of 6 cm chord (c). The blades rotate at constant speed ðUÞ at a radius (R) of 0.5 m from the central axis, and the attachment points are located at 0.4c from the leading edge. Fig. 2 shows the 3D drawing of the real turbine used in the experiment. A simplified geometry consisting of the turbine's blades only has been used in the numerical analysis. No central shaft and no supporting struts are included in the computation. The results presented refer to the higher Tip Speed Ratio (TSR) of 4.5, where dynamic stall can be excluded. Moreover the incoming wind ðV∞ Þ is supposed to be uniform without any turbulence added. 3.4. Discretization and sensitivity analysis A sensitivity analysis has been conducted to determine the influence of the temporal and spatial discretization of the model. Moreover, since the model is unsteady, an initial transient develops before the flowfield reaches a converged periodic solution. An effect similar to the starting vortex of an airfoil is observed with the starting wake rolling up behind the rotor and the induction field at the rotor level gradually building up as more wake is released. The required minimum wake length for a converged periodic solution is also evaluated here. Four independent parameters have been considered: the azimuth step ðdqÞ, for the temporal discretization; the number of simulated rotations (ROT), for the effect of the starting wake and the build-up of the induction field; the number of spanwise panels (SPAN) and chordwise panels (CPAN), for the spatial discretization. A simplified one-factor-at-a-time approach is used, where for each
555
study only one of the four parameters is varied while the others are kept at a baseline value. The perceived angle of attack of the section at the middle of the blade (az0(q)) is used as the evaluation output. This is determined from the induced velocity at 6 points distributed around the quarter chord point, and represents a more detailed choice than using an integral parameter such as the thrust or power coefficient. Fig. 3 shows the sensitivity analysis for the azimuth step, ranging from Ref. dq ¼ 10 (blue line) to symb:deltatheta ¼ 1 (black line), for a standard simulation (30 panels chordwise, 30 panels spanwise, 15 rotations). In the plot two parts with different behavior can be identified: the upwind part of the blade rotation (0+ q 180+ ), characterized by a smooth sinusoidal-like distribution of az0(q), and the downwind part of the rotation (180 q 360 ), where multiple bwi cause a noisy, relatively constant, distribution of az0(q). The four lines in the figure show overall small deviations. Apart from some small regions, differences in Ref. az0(q) are bounded to 0.2+ maximum. Since there is no sensible difference between the Dq ¼ 5 and the Dq ¼ 2 , while the computational cost increases by almost ten times, Dq ¼ 5 is chosen for the following simulations. Fig. 4 shows convergence analysis for the simulated rotations, up to 15, for a standard simulation (30 panels chordwise, 30 panels spanwise, Dq ¼ 5 ). In the upwind portion of the rotor, convergence is reached at 7 rotations, after that the differences in Ref. az0(q) are bounded at 0.1. In the downwind part the convergence is less monotonic with high fluctuations due to the intense bwi which differs for each realization. However a good degree of convergence is also reached at 7 rotations with variations in calculated angles of attack around 0.2 . The influence of the spatial discretization is investigated by refining the number of panels on the blade, both in chordwise and spanwise direction. Increasing the number of spanwise panels increases also the number of elements in the wake resulting in a significant increase in computational time; while increasing the number of chordwise panels has a smaller effect. Fig. 5 shows convergence analysis for the chordwise panels. The chordwise discretization seems to be of influence, with the first refinement (from 16 to 30 cpan) being of higher impact than the second (from 30 to 60 cpan), and 30 chordwise panels already giving a good convergence. Fig. 6 shows the effect of a refinement in the number of spanwise panels. A first refinement from 16 to 30 spanwise panels decreases the az0(q) at q ¼ 90 of circa 0.1 (with an increase of computational time of 3). A similar decrease in the calculated angle of attack is achieved by a further refinement from 30 to 60 panels (with an increase of computational time of 8.5). Doubling again the number of spanwise panels leads the az0(q) to higher values. A value of 60 spanwise panels seems to be the most accurate solution (within the cases considered), while reducing to 30 panels lowers the accuracy minimally with a considerable gain in computational time. 3.5. Numerical stability The sensitivity analysis shows a similar non-monotonic behavior for the spanwise spatial or temporal discretizations. When the discretization is refined (doubling the number of spanwise panels or halving the azimuth step), az0(q) first decreases to then increase again (Fig. 7 shows this behavior for the blade at the most upwind position). While the initial decrease is considered to be caused by the reduction of the discretization error, the subsequent increase is attributed to the effect of a numerical instability of the model and analyzed in the present paragraph. The flowfields show localized regions of high values of velocity and velocity gradients. Such regions start to appear around 1.5
556
G. Tescione et al. / Renewable Energy 87 (2016) 552e563
Fig. 3. Effect of the azimuth step ðDqÞ on the angle of attack at the middle of the blade (az0(q)) across the azimuth (q).
diameters downstream the turbine shaft in the center of the wake, with small differences depending on the discretization parameters. The study of the wake geometry highlighted pronounced deformations in the central part of the wake. These are the results of strong close interactions between wake elements and blade panels and mutually between vortex filaments. A higher refined temporal or spatial grid results in a higher probability of such instabilities
which tend to be amplified by the Lagrangian approach and the explicit time scheme. A simulation where the body influence had been excluded led to a very limited improvement in stability, revealing how the vortex filaments have a dominant effect. At the rotor level the discretization error reduction is limited and a refinement of the spatial or temporal grid does not lead to a better solution. At the wake level the solution is stable and rapidly
Fig. 4. Effect of the simulated number of rotations (rot) on the angle of attack at the middle of the blade (az0(q)) across the azimuth (q).
G. Tescione et al. / Renewable Energy 87 (2016) 552e563
557
Fig. 5. Effect of the number of chordwise panel (cpan) on the angle of attack at the middle of the blade (az0(q)) across the azimuth (q).
converged until a certain point, following the motion of the intense structures of the tip vortices, after which high gradients and fluctuations of velocity are present. The stability of the wake has been analyzed by evaluating the absolute value of the velocity (according to Equation (3)) on the median line of the wake (y ¼ 0,z ¼ 0) for each azimuth position of the final simulated rotation, and calculating both the averaged value
and standard deviation. Fig. 8 shows the velocity profile for one of the simulated cases (RB1, Dq ¼ 5 ,span ¼ 30,cpan ¼ 30,rot ¼ 8). The plots show the upwind (A) and downwind (C) blade passages (with high temporal fluctuations due to the position of the blades), the interactions with the vortex filaments inside the rotor area (B) and in the near wake (D), and the point where the vortex filaments of the tip vortices intersect the median line causing high gradients and
Fig. 6. Effect of the number of spanwise panel (span) on the angle of attack at the middle of the blade (az0(q)) across the azimuth (q).
558
G. Tescione et al. / Renewable Energy 87 (2016) 552e563 Spanwise Panels 15
30
60
120 Spanwise Panels
α
90
[deg]
8.20
8.15
8.10
8.05
Δθ
10
Δθ [deg]
5
2
1
Fig. 7. Variation of the perceived angle of attack at the middle of the blade in the most upwind position (az0 90 at the refinement of the spatial (span) and temporal (Dq) discretization.
Fig. 8. Absolute velocity profiles, average and standard deviation on the wake median line (y ¼ 0,z ¼ 0) for the RB1 run (Dq¼5+,span ¼ 30,cpan ¼ 30,rot ¼ 8). Upper plot shows the velocity profiles for the different blade azimuth (V(q), thin gray lines) and the velocity profile averaged across one full rotation (VAVG, thick black line). Lower plot shows the standard deviation (std(V), thick red line) across one full rotation. All velocities are made non-dimensional with the free-stream velocity.
fluctuations of velocity (E). Table 1 presents the results of the stability analysis for all the 10 different runs analyzed in Section 3.4 for the sensitivity analysis. The results are grouped by the 3 parameters defining the discretization (Dq,span,cpan). The limit of stability (XST) has been measured as the location where the standard deviation of the absolute velocity over one rotation starts to peak and is expressed in rotor diameters (D) from the turbine shaft. The table shows how a coarser blade paneling or higher azimuth steps lead to a more stable wake up to 1.5e1.7 diameter downwind of the turbine shaft, while the finer paneling or smaller time step reduce the stable part of the wake down to 1.1e1.3 diameter downwind. More work is needed, and undergoing, to improve the
Table 1 Limit of stability for the different runs, effect of the discretization parameters. Effect of Dq
Effect of span
Effect of cpan
Run
Dq
XST
Run
Span
XST
Run
Cpan
XST
RB5 RB1 RB3 RB6
10+ 5+ 2 1
1.45 1.38 1.24 1.06
RB8 RB1 RB2 RB7
15 30 60 120
1.56 1.38 1.40 1.32
RB10 RB1 RB4 RB9
16 30 60 120
1.71 1.38 1.35 1.12
convergency and stability of the further part of the wake solution.
4. Model validation From the results of the sensitivity and stability analysis (Sections 3.4 and 3.5 respectively), the model used for the validation has a 5 azimuth step, 30 panels chordwise, 30 panels spanwise and a 7 rotations transient, in order to reach a trade-off between convergence, stability and computational cost. The solution is considered to be stable up to 1.4 diameters downwind the turbine shaft in the mid-plane. On the workstation used2 the computational time is circa 30 min. The experimental campaign [21] provided phase-locked 2D velocity fields at the turbine mid-span horizontal plane and 3D velocity fields at 7 vertical planes spanning in the cross-stream direction (y/ D ¼ 0.5, 0.4, 0.2, 0, 0.2, 0.4, 0.5) in the rotor area and till 3 diameters downstream. The reader is referred to [21] for a detailed phenomenological explanation of the nature of the vortical
2 intel i7-950 with 8 cores @ 3 GHz, 12 GB RAM, and NVIDIA® Quadro 4000 GPU with 256 cores, 243GFlops in double precision and 89.6 GB/s memory bandwidth.
G. Tescione et al. / Renewable Energy 87 (2016) 552e563
559
Fig. 9. Contour levels of streamwise velocity (Vx/V∞) in the mid-plane (z ¼ 0) for the experimental (left) and numerical (right) results.
structures in the vawt wake that would help understanding the following analysis, the same terminology and conventions are also used. 4.1. Horizontal plane The first comparison between the simulations and the experimental data is about the velocity field in the horizontal mid-plane. Figs. 9 and 10 show the streamwise and cross-stream velocity components (Vx and Vy respectively) for the horizontal mid-plane (z ¼ 0). The experimental data are shown on the left and the numerical data on the right. The velocities are made non-dimensional by the free stream V∞ and the experimental and numerical contour plots share the same contour levels for an easy comparison. The streamwise velocities (Fig. 9) are very well predicted for the upwind (Dx 0:5) and in-rotor (0:5 Dx 0:5) area with the exclusion of the region just downstream of the turbine shaft (Dx 0; Dy z0). The experimental results show a dominant viscous wake of the rotating cylinder while in the numerical results this is absent since the shaft is not simulated. The numerical solution of the downwind part of the wake (Dx 0:5) presents a good comparison in terms of wake extension for the leeward side (Dy 0) and a lower expansion for the windward side (Dy 0). The latter is partially due to the absence of the rotating shaft wake. Moreover, the numerical results show a more extended region of higher wake deficit (VV∞x ¼ 0:2) as well as higher acceleration of the flow outside
the wake (VV∞x ¼ 1:1). This may be explained with a higher strength of the vortical structures at the boundaries of the wake and it is also partly direct consequence of the absence of any decay of vorticity in the numerical model. It should be noted that the numerical results are instantaneous phase-locked fields, while the experimental results are phase-locked average. The process of averaging introduces a smoothing effect of the near-field inductions of vortical structures due to a certain degree of phase unsteadiness [21]. The cross-stream velocity plots (Fig. 10) show a good quantitative agreement of the contours in all the shown domain. In the rotor area (0:5 Dx 0:5) the vortex elements follow closely the position of the blades wake, as well as the vortical structures downwind the rotor (Dx 0:5) in the leeward side (Dy 0). Some differences are present in the velocity contours and in the positions of the structures in the windward side (Dy 0). 4.2. Vertical planes Figs. 11e13 show the streamwise, cross-stream and vertical velocity components (Vx, Vy, and Vz respectively) for the vertical planes (y/D ¼ 0.5, 0.4, 0.2, 0, 0.2, 0.4, 0.5). As for the previous case, the velocities are made non-dimensional by the free stream V∞ and the same contour levels are used for the experimental (left) and numerical (right) results. Observing the contour plots an overall good agreement can be seen between experimental and numerical results. Contour levels
Fig. 10. Contour levels of cross-stream velocity (Vy/V∞) in the mid-plane (z ¼ 0) for the experimental (left) and numerical (right) results.
560
G. Tescione et al. / Renewable Energy 87 (2016) 552e563
Fig. 11. Contour levels of streamwise velocity (Vx/V∞) in the vertical planes (y/D ¼ 0.5, 0.4, 0.2, 0, 0.2, 0.4, 0.5) for the experimental (left) and numerical (right) results.
are well predicted for the streamwise velocity (Fig. 11) and especially for the cross-stream and vertical velocity (Figs. 12 and 13), with some local discrepancies. Also the dynamics of the tip vortices is correctly captured with their inboard motion in the middle part of the wake, already observed in experiments and simulations [8,14,15,21]. The most noticeable difference is the vertical expansion of the
windward side of the wake (Dy 0) where the experiments showed a pronounced vertical expansion, higher than the one in the leeward side (Dy 0), while the numerical solution shows a reduced expansion comparable to the one on the leeward side. A similar result was observed for the horizontal mid-plane. A second mismatch is in the local induction of the tip vortices. The numerical flowfield shows stronger local induction, observable
G. Tescione et al. / Renewable Energy 87 (2016) 552e563
561
Fig. 12. Contour levels of cross-stream velocity (Vy/V∞) in the vertical planes (y/D ¼ 0.5, 0.4, 0.2, 0, 0.2, 0.4, 0.5) for the experimental (left) and numerical (right) results.
in Fig. 11 in the higher accelerated flow outside the wake (VV∞x ¼ 1:1) or in Fig. 13 in the higher contours associated to the vortical structures. The numerical results suggest stronger released tip vortices which conserve their strength as they travel downstream, while the experimental results have shown a fast decay of the vortical structures as soon as the upwind-released vortex interacts with the downwind-released one [21].
Also the absence of the turbine shaft in the Dy ¼ 0 plane is noticed by the absence of its wake as well as the absence of the two struts wakes in the Dy ¼ 0:5 (windward) and Dy ¼ 0:5 (leeward) planes. The cross-stream and vertical velocity fields give more insight in the wake dynamics not correctly predicted by the numerical model and explain some of the differences found in the horizontal mid-
562
G. Tescione et al. / Renewable Energy 87 (2016) 552e563
Fig. 13. Contour levels of vertical velocity (Vz/V∞) in the vertical planes (y/D ¼ 0.5, 0.4, 0.2, 0, 0.2, 0.4, 0.5) for the experimental (left) and numerical (right) results.
plane. The overall behavior of the wake is well captured. The horizontal expansion in the cross-stream velocity contours (yellow regions of positive values in the windward side and blue regions of negative values in the leeward side) and the vertical motion in the vertical velocity contours (yellow regions of positive value in the windward and leeward side and blue regions of negative values in the central part of the wake) are similar. Some regions present significant quantitative discrepancies: the windward side of the
wake downstream the rotor (Dy 0 and Dx 0:5) shows significant lower cross-stream (Fig. 12) and vertical (Fig. 13) velocities. The reason is to be found again in the higher strength of the vortical structures, opposed to the experimental case (Rec ¼ 170000), causing a stronger local induction velocity field. The observed reduced expansion of the wake of the windward side both in the horizontal and in the vertical planes confirms this. Also the vertical inboard motion of the central part of the wake observed in the
G. Tescione et al. / Renewable Energy 87 (2016) 552e563
experiment is reduced in the numerical analysis. A last observation concerns the numerical stability of the wake in the vertical planes. In Section 3.5 the analysis focused on the mid-plane, showing a wake numerically stable until at least Dx ¼ 1:4. Figs. 12 and 13 show signs of numerical instabilities (represented by high local values of induced velocities and velocity gradients) developing around the tip vortices at smaller downwind positions. Such instabilities, though of a rather localized nature, are also partially responsible of the local differences between numerical and experimental results observed in the present analysis.
563
comparison between a numerical model without diffusion scheme and an experimental model run at a low Reynolds number. Moreover the absence of the struts in the simulation may be responsible for an initial overshoot of the tip vortex strength and the absence of the turbine shaft for a part of the wake asymmetry. It is to believe that simulations of higher Reynolds number cases, comprising the struts and the turbine shaft will see their accuracy increased, making the model a valuable tool for the analysis of multi-MW VAWTs. References
5. Conclusions A free vortex wake model has been verified and validated with experimental measurements of a vawt wake. The model has been analyzed in terms of sensitivity of the output to the simulation parameters. For the flow inside the rotor (for loads evaluations of a single turbine), a medium discretization (30 chordwise panels, 30 spanwise panels, 5 of azimuth step) gives a converged solution in terms of angles of attack experienced by the blade. The initial transient to reach the periodic converged solution, for the present application, is around 7 rotations. On a normal workstation equipped with a gpu such simulation leads to a computational time of 30 min. Any further grid refinement, by doubling the blade panels or halving the azimuth step, or any added simulated rotation leads to a marginal change in the angle of attack. The simulation has also been analyzed in terms of numerical stability. The tip vortex evolution represents the most critical factor for the stability of the model. Given the high local induction of these vortical structures and the Lagrangian formulation of the model, the wake is subject to strong local deformations leading to high velocity gradients. A second-order time scheme and the desingularization of the viscous vortex core help controlling such effect. Moreover the grid refinement, while reducing the discretization error by a limited amount, decrease the stability of the model. The effects of such instabilities are rather local and do not travel upwind, so the rotor and the near-wake flows are less affected. With the parameters used in the present analysis, the numerical solution in the wake is considered stable up to 1.4 diameters downstream the turbine shaft in the symmetry plane. The numerical results have been compared to the measurements of a previous experimental activity on a 2-bladed hvawt, running at a tsr of 4.5 and a chord-based Reynolds number of 170000. The model is capable of capturing all the important dynamics in the wake of vawt and give also good quantitative evaluation of the flowfield. Comparisons show very good prediction of the flowfield upwind and inside the rotor area. The dynamics of the vortical structures in the horizontal plane is very well captured as well as the vertical motion of the tip vortices with the already observed central contraction and lateral expansion. The results in the near wake (0:5 Dx 1:5) show good prediction of the wake flow in the leeward side but tend to underestimate the wake horizontal and vertical expansion in the windward side compared to the experimental results, decreasing the high asymmetry observed in the experiments. Also the local induction of the vortical structures in the far part of the wake (both of the horizontal plane and of the tip vortices) is higher in the numerical results. The reasons for such discrepancies are to be found in the
rimentale et Nume rique du De crochage Dynamique [1] P.L. Beaudet, Etude Expe (Ph.D. thesis), Universit de sur une Eolienne a Axe Vertical de Forte Solidite Poitiers, 2014. [2] M. Borg, M. Collu, F.P. Brennan, Offshore floating vertical axis wind turbines: advantages, disadvantages, and dynamics modelling state of art, Mar. Offshore Renew. Energy (2012) 33e46. [3] M. Borg, A. Shires, M. Collu, Offshore floating vertical axis wind turbines, dynamics modelling state of the art, Part I Aerodyn. Renew. Sustain. Energy Rev. 39 (2014) 1214e1225. [4] P. Deglaire, Analytical Aerodynamic Simulation Tools for Vertical Axis Wind Turbines (Ph.D. thesis), Uppsala University, 2010. [5] K.R. Dixon, The Near Wake Structure of a Vertical Axis Wind Turbine (Master's thesis), Delft University of Technology, 2008. [6] J.L. Hess, Panel methods in computational fluid dynamics, Annu. Rev. Fluid Mech. 22 (1990) 255e274. [7] J.L. Hess, A.M.O. Smith, Calculation of potential flow about arbitrary bodies, Prog. Aeronaut. Sci. 8 (1967) 1e138. [8] C. Hofemann, C.J. Sim~ ao Ferreira, K. Dixon, G.J. van Bussel, G. van Kuik, F. Scarano, 3D stereo piv study of tip vortex evolution on a vawt, in: Proceeding of EWEC, Brussels, 2008. [9] J. Katz, A. Plotkin, Low-speed Aerodynamics: from Wing Theory to Panel Methods, McGraw-Hill, New York, 1991. [10] L. Morino, C.C. Kuo, Subsonic potential aerodynamics for complex configurations: a general theory, AAIA J. 12 (2) (1974) 191e197. [11] J.C. Murray, M. Barone, The development of cactus: a wind and marine turbine performance simulation code, in: 49th AIAA Aerospace Sciences Meeting, Orlando, FL, 2011. [12] P.J. Musgrove, Wind energy conversion: recent progress and future prospects, Sol. Wind Technol. 4 (1) (1987) 37e49. [13] S. Peace, Another approach to wind, Mech. Eng. 126 (6) (2004) 28e31. [14] F. Scheurich, T.M. Fletcher, R.E. Brown, Simulating the aerodynamic performance and wake dynamics of a vertical-axis wind turbine, Wind Energy 14 (2) (2011) 159e177. ~o Ferreira, The Near Wake of the VAWT, 2D and 3D Views of the [15] C.J. Sima VAWT Aerodynamics (Ph.D. thesis), Delft University of Technology, 2009. [16] C.J. Sim~ ao Ferreira, H. Aagard Madsen, M. Barone, B. Roscher, P. Deglaire, I. Arduin, Comparison of aerodynamic models for vertical axis wind turbines, J. Phys. Conf. Ser. 524 (1) (2014). ~o Ferreira, F. Scheurich, Demonstrating that power and instantaneous [17] C.J. Sima loads are decoupled in a vertical-axis wind turbine, Wind Energy 17 (2014) 385e396. [18] M.J. Stock, A. Gharakhani, Toward efficient gpu-accelerated n-body simulations, in: 46th AIAA Aerospace Sciences Meeting and Exhibit, January 2008. [19] J.H. Strickland, B.T. Webster, T. Nguyen, A vortex model of the darrieus turbine: an analytical and experimental study, J. Fluids Eng. 101 (1979) 500e505. [20] H.J. Sutherland, D.E. Berg, T.D. Ashwill, A Retrospective of Vawt Technology, Technical report, Sandia National Laboratories, January 2012. [21] G. Tescione, D. Ragni, C. He, C.J. Sim~ ao Ferreira, G.J.W. van Bussel, Near wake flow analysis of a vertical axis wind turbine by stereoscopic particle image velocimetry, Renew. Energy 70 (2014) 47e61. [22] L. Vita, U.S. Paulsen, T.F. Pedersen, H.A. Madsen, F. Rasmussen, Deep wind: a novel floating wind turbine concept, Wind Int. 6 (4) (2010) 29e31. [23] Vv.Aa, Wind in Our Sails. The Coming of Europe's Offshore Wind Energy Industry, Technical report, European Wind Energy Association (EWEA), November 2011. [24] Vv.Aa, Deep Water. The Next Step for Offshore Wind Energy, Technical report, European Wind Energy Association (EWEA), July 2013. [25] D.J. Willis, An Unsteady, Accelerated, High Order Panel Method with Vortex Particle Wakes (Ph.D. thesis), Massachusetts Institute of Technology, 2006. ~o Ferreira, Wake modelling of a vawt in [26] A. Zanon, P. Giannattasio, C.J. Sima dynamic stall: impact on the prediction of flow and induction fields, Wind Energy 18 (11) (November 2015) 1855e1874, http://dx.doi.org/10.1002/ we.1793.
Contents lists available at ScienceDirect
Renewable Energy journal homepage: www.elsevier.com/locate/renene
Analysis of a free vortex wake model for the study of the rotor and near wake flow of a vertical axis wind turbine ~o Ferreira, G.J.W. van Bussel G. Tescione*, C.J. Sima Delft University Wind Energy Research Institute, TUDelft, Kluyverweg 1, 2629 HS, Delft, The Netherlands
a r t i c l e i n f o
a b s t r a c t
Article history: Received 16 February 2015 Accepted 1 October 2015 Available online xxx
The 3D dynamics of the near wake of a vertical axis wind turbine is simulated by a numerical model based on a free vortex wake/panel method. The model is first verified with respect to discretization sensitivity, stability and robustness. Results show little sensitivity to the discretization, requiring a medium spatial and temporal grid to reach convergence. Stability issues are addressed, identifying dependency of the simulation parameters and validity of the results. The in-rotor flow shows little or no effect of instabilities and the model proves to be an effective analysis tool for the flow at the blade level. Wake dynamics, especially the tip vortex evolution, is more subject to stability issues which have to be taken into consideration. The numerical results are then validated with experimental data to assess accuracy. The model is able to capture all the important dynamics in the near wake of a vertical axis wind turbine with good quantitative evaluation of the flowfield. Discrepancies in the comparison are discussed and limitations in the use of the model are presented. © 2015 Elsevier Ltd. All rights reserved.
Keywords: Vertical axis wind turbine Vortex element methods Panel method Numerical aerodynamics Wake aerodynamics Vorticity dynamics Verification and validation
1. Introduction Future trends of wind energy industry show the development of large turbines in offshore wind farms. As reported by the European Wind Energy Association [23], the latest generation of offshore wind turbines has rotor diameters up to 170 m and rated capacity up to 7 MW, while 10e20 MW machines are currently in the development phase. Offshore wind farm development trends towards bigger farms built further from te coast and in deeper waters. Going beyond the current limit of 50 m of water depth for fixed substructure concepts will allow to expand the market in deeper seas, as the remaining 66% of the North Sea, the Mediterranean and the Atlantic basins in Europe, the coast of North America, China or Japan. Within these waters, floating support structures are shown to be more economical [24]. The development of large floating wind farms of multi-MW turbines has made Vertical Axis Wind Turbines (VAWTs) recently gain new popularity in the wind energy community, as their use
Acronyms: BWI, BladeeWake Interaction; CPAN, chordwise panels; GPU, Graphic Processing Unit; H-VAWT, H-shaped VAWT; HAWT, Horizontal Axis Wind Turbine; ROT, simulated rotations; SPAN, spanwise panels; TSR, Tip Speed Ratio; VAWT, Vertical Axis Wind Turbine. * Corresponding author. E-mail address: [email protected] (G. Tescione). http://dx.doi.org/10.1016/j.renene.2015.10.002 0960-1481/© 2015 Elsevier Ltd. All rights reserved.
can potentially alleviate some of the new challenges [2,12,13]. VAWTs are mechanically simpler machines than Horizontal Axis Wind Turbines (HAWTs), having no yaw, and often no pitch, mechanisms. In an offshore environment, where operation and maintenance have a high weight in the Cost of Energy, the presence of fewer moving parts has a great value. Secondly, VAWTs show higher potentials of scalability, overcoming some of the limitations of HAWTs in this aspect. Moreover the possibility to place the generator at sea level or even under water helps stability and would decrease the size and cost of the floating support structure. Different research programmes in Europe and North America [20,22], are exploring the potential of floating large VAWTs. This renewed interest finds however a fragmented research field. After its peak during the 80's, VAWT research was discontinued for almost 15 years, during the 90's and early years of 2000. The commercial success of the propeller type horizontal axis turbines drove most of the research effort of the wind energy community, with few isolated researchers focusing on the vertical axis turbines, and their use limited to small turbines for urban environment. VAWTs missed the extensive aerodynamic research characterizing HAWTs, resulting in a lack of a solid phenomenological background and an underdevelopment of proper models to assist analysis and design [3,16]. The aerodynamics of vertical axis turbines differs from that of the horizontal axis devices in several fundamental aspects. The
G. Tescione et al. / Renewable Energy 87 (2016) 552e563
List of symbols
az0(q) az0j90 m s f f∞ Vy ! u GTE g gS gT
Dq q c D H M
perceived angle of attack of the middle section of the blade az0(q) at the most upwind position (q ¼ 90 ) doublet singularity strength source singularity strength velocity potential undisturbed velocity potential rotational velocity vorticity vector circulation on the trailing edge segment vortex lament strength shed vorticity (wake filament) trailing vorticity (wake filament) azimuth step azimuth angle blade chord rotor diameter rotor height Mach number
N ! n R Rec std(V) ! s t ! u ! ui V∞ VAVG V (q) Vx Vy Vz XST x y z
553
number of bodies in N-body problem normal vector (to body/wake surface) rotor radius chord-based Reynolds number standard deviation of velocity trailing edge segment vector time velocity vector (general form) induced velocity vector (general form) free-stream velocity averaged absolute velocity instantaneous absolute velocity stream-wise velocity component cross-stream velocity component vertical velocity component limit of stability (in D) for the wake stream-wise extension stream-wise spatial coordinate cross-stream spatial coordinate vertical spatial coordinate
3. Numerical model
energy conversion process is based on the variation of the blade's circulation along its rotation, resulting in a different creation of the wake. Unsteadiness is an inherent part of the turbine operation and its wake is characterized by regions of concentrated vorticity interacting and evolving in coherent structures [15,21]. Accurate wake dynamics prediction is fundamental in VAWT design process as part of the rotor operates in its own wake. Moreover in wind farm design the far wake evolution is of high interest. Low-fidelity engineering models developed for HAWT have been adapted for VAWT design and analysis without considering the inherent differences, resulting in poor predictions [15e17]. Moreover they lack the details of the wake development. Higherorder, grid-based numerical methods tend to underestimate the role and the evolution of the vortical structures due to their inherent numerical dissipation. Vortex element methods have been used in VAWT research for their ability to retain the strength of the vortical structures and to give details of the wake evolution [1,4,11,19,25,26].
The model used for this numerical investigation is 3D, unsteady and based on a panel method coupled with a free vortex wake. Aerodynamic panel methods have been known since the 1960s and extensively used in the aeronautic field as analysis and design tools.1 The model is based on the assumption of an incompressible and potential (inviscid, irrotational) flow, with the exception of vorticity singularities introduced in the wake elements and on the body. Given the relative low Mach number (M < 0.3) and high Reynolds numbers (O(Rec)) experienced by VAWTs such assumptions are considered to be reasonable everywhere in the flow outside the boundary layer. The potential flow assumption allows for the use of a Boundary Integral Equation for the velocity potential (f, the velocity being ! u ¼ Vf) expressed in one of its general forms as:
2. Approach
f¼
The present paper investigates the use of an unsteady free vortex-wake model developed at TUDelft [5] for the analysis of the aerodynamics of the wake of a VAWT. First, the model is introduced by the mathematical formulation, the numerical implementation used and the simulated VAWT operation. Then, a verification of the model is performed addressing convergence and stability of the solution. Geometrical discretizations of the paneled body and of the vortex wake and temporal discretizations are related to the perceived angle of attack of the blade section at the symmetry plane, to assess sensitivity of the results. The effect of the initial transient solution originated by the impulsive start of the simulation is also determined and guidelines on the duration of such transient are given. Stability issues concerning the evolution in the wake downwind the turbine are also addressed to evaluate the applicability of the model. Finally, the results of the verified model are validated against the experimental measurements presented and discussed by the present authors in Ref. [21].
3.1. Mathematical formulation
1 4p
Z body
1 v 1 1 sþ m dS þ r vn r 4p
Z wake
v 1 mdS þ f∞ vn r (1)
where s and m are the singularities on the body and wake surface (S) and r is the distance between the singularity and the location of the potential evaluation. The present model uses a coupled formulation of sources ðsÞ and doublets ðmÞ on the actual body geometry for the singularity mix [10]. Here the doublet strength is taken as independent on each panel and equals the exterior surface value of perturbation potential while the source strength equals its ! ! ! normal derivative ( n ,Vf ¼ n , u ) to fulfill the zero-normalvelocity boundary condition. To solve for the non-uniqueness of the solution the Kutta condition, that the pressure difference vanishes at the trailing edge, is used, which drives the entire solution
1 For a review on the history of the methods and overview of the different formulations the reader is referred to Hess [6,25], and the book from Katz & Plotkin [9].
554
G. Tescione et al. / Renewable Energy 87 (2016) 552e563
by determining the circulation distribution. The wake is modeled as a curved lattice of straight vortex filaments. The wake contribution to the velocity potential in Equation (1) (the second integral) is shown as a doublets influence; the actual model shifts from doublets to vortex filaments ðgÞ after few time steps. This is done as the use of a vortex lattice allows for a better handling of deformations of the wake geometry. The equivalence of the doublets and vortex distribution is well known [7]. Given the incompressibility assumption the time dependency of the solution follows straightforwardly from the time-dependent boundary conditions, hence the unsteady nature of the model. This is of primary importance given the inherent unsteadiness of VAWT aerodynamics [15,21]. The use of vorticity to describe the flow brings important advantages when associated to a Lagrangian approach and modifies the governing equation to a simplified form. The further hypothesis of no viscous diffusion and incompressible flow lead to the following formulation for the transport equation of the vorticity:
D! u ! ð! u $VÞ u ¼ 0 Dt
(2)
However the non-linearity of the original NaviereStokes equations is maintained in the vorticity formulation as the resulting equation describes a non-linear transport phenomenon, since the ! transport velocity ð u Þ depends itself on the transport quantity (the vorticity ! u ). The vorticity field induces a velocity field which is responsible for the evolution of the Lagrangian particles. The recovery of the velocity from the vorticity field is obtained with the integral approach. The Poisson equation for the velocity ! (V2 u ¼ V ! u ) is solved by means of the Green's function for the Laplacian. This leads to the Biot-Savart law for the velocity induc! tion ð u i Þ by a vorticity field ð! u Þ:
1 ! ! u i ð x ; tÞ ¼ 4p
Z !! ! ! u ð y ; tÞ ð x y Þ ! dy ! !3 x y D
(3)
which is the basis for the wake-to-wake and wake-to-body induction calculation. The advantage of a vorticity formulation rather than a velocity formulation is that vorticity has a compact support, so the integration is limited to the domain D where vorticity is nonzero while the velocity field can be evaluated on a larger domain.
3.2. Numerical implementation The numerical model is implemented in hybrid MATLAB®/ CUDA®-C code for an easier use with an Object Oriented Programming approach to increase flexibility. The code is, indeed,
capable of handling multi-body objects in 2D or 3D of generic geometries and different arbitrary (and time-dependent) motions [5]. The rotor blades are discretized with flat panels of constant source and doublet strengths. The geometric discretization is refined towards the leading and trailing edges and towards the tips of the blades to account for curvature and higher pressure (velocity) gradients. No deformation of the body is accounted for, so the blades are subject to a rigid rotation around the vertical axis. This allows to build in a pre-processing phase the Aerodynamic Influence Coefficients of each blade panel to each other, and to use matrix multiplication during the simulation to compute the bodyto-body induction. The wake is modeled as a curved lattice of straight vortex filaments (see Fig. 1). The strength of the wake elements released each time step is derived by the linearized Kutta condition at the blade ! trailing edge segments (along the coordinate s ). Spatial derivatives of the circulation between adjacent trailing edge segments ( DG !) Ds determine the trailing vorticity ðgtÞ of the wake filaments perpendicular to the trailing edge. Temporal derivatives (DG Dt ) determine the shed vorticity ðgSÞ of the wake filaments parallel to the trailing edge. Wake elements are emitted only by the trailing edges, no wake elements are emitted in the spanwise direction from the tip of the blades. The new wake marker position is defined by calculating the local kinematic velocity at the trailing edge. A free-wake analysis is used where the new position of each wake marker is an integral part of the problem. This is essential in the cases of VAWT as the rotor operates partly in its own wake. For each time step the induction field of all the wake filaments and of the body panels is solved on the location of the wake markers (the filaments' endpoints) by a direct summation. The filaments have a viscous core with a Rankine-type behavior to ensure desingularization of Equation (3). No viscous diffusion of the filament vorticity is currently included in the model. As the endpoints of the filaments are used as markers, vortex stretching and tilting is auto! matically considered (the ð! u ,VÞ u term in Equation (2)). Because of its completely Lagrangian formulation the model requires no grid creation and avoids numerical dissipation of vorticity. An explicit second order Adam-Bashforth method is used for the time marching scheme. No wake cut-off is used: all wake elements are retained for the whole simulation. The resulting N-body problem (with an O(N2) computational cost and an OðNÞ memory access) is very suitable for computational parallelization, especially with the use of Graphic Processing Unit (GPU) -based systems [18]. The most computational expensive functions, the wake-to-wake and waketo-body induction calculations, have been implemented in CUDA®-C and run parallel on a NVIDIA® Quadro 4000 GPU (256 cores, 243GFlops in double precision and 89.6 GB/s memory
Fig. 1. Schematics of the vortex lattice wake discretization and the different contributions to the wake vorticity.
G. Tescione et al. / Renewable Energy 87 (2016) 552e563
Fig. 2. 3D rendering of the VAWT model used in the experiments (dimensions in mm).
bandwidth). This brings to a speed up of the calculations up to 120 with respect to a pure MATLAB® execution for a mid-refined simulation. This allows for higher refined and/or long-time simulations. 3.3. VAWT geometry and operation A model H-shaped VAWT (H-VAWT) of a previous experimental campaign [21], is used for the current calculations. It consists of 2 straight blades of 1 m height (H) with a NACA0018 airfoil of 6 cm chord (c). The blades rotate at constant speed ðUÞ at a radius (R) of 0.5 m from the central axis, and the attachment points are located at 0.4c from the leading edge. Fig. 2 shows the 3D drawing of the real turbine used in the experiment. A simplified geometry consisting of the turbine's blades only has been used in the numerical analysis. No central shaft and no supporting struts are included in the computation. The results presented refer to the higher Tip Speed Ratio (TSR) of 4.5, where dynamic stall can be excluded. Moreover the incoming wind ðV∞ Þ is supposed to be uniform without any turbulence added. 3.4. Discretization and sensitivity analysis A sensitivity analysis has been conducted to determine the influence of the temporal and spatial discretization of the model. Moreover, since the model is unsteady, an initial transient develops before the flowfield reaches a converged periodic solution. An effect similar to the starting vortex of an airfoil is observed with the starting wake rolling up behind the rotor and the induction field at the rotor level gradually building up as more wake is released. The required minimum wake length for a converged periodic solution is also evaluated here. Four independent parameters have been considered: the azimuth step ðdqÞ, for the temporal discretization; the number of simulated rotations (ROT), for the effect of the starting wake and the build-up of the induction field; the number of spanwise panels (SPAN) and chordwise panels (CPAN), for the spatial discretization. A simplified one-factor-at-a-time approach is used, where for each
555
study only one of the four parameters is varied while the others are kept at a baseline value. The perceived angle of attack of the section at the middle of the blade (az0(q)) is used as the evaluation output. This is determined from the induced velocity at 6 points distributed around the quarter chord point, and represents a more detailed choice than using an integral parameter such as the thrust or power coefficient. Fig. 3 shows the sensitivity analysis for the azimuth step, ranging from Ref. dq ¼ 10 (blue line) to symb:deltatheta ¼ 1 (black line), for a standard simulation (30 panels chordwise, 30 panels spanwise, 15 rotations). In the plot two parts with different behavior can be identified: the upwind part of the blade rotation (0+ q 180+ ), characterized by a smooth sinusoidal-like distribution of az0(q), and the downwind part of the rotation (180 q 360 ), where multiple bwi cause a noisy, relatively constant, distribution of az0(q). The four lines in the figure show overall small deviations. Apart from some small regions, differences in Ref. az0(q) are bounded to 0.2+ maximum. Since there is no sensible difference between the Dq ¼ 5 and the Dq ¼ 2 , while the computational cost increases by almost ten times, Dq ¼ 5 is chosen for the following simulations. Fig. 4 shows convergence analysis for the simulated rotations, up to 15, for a standard simulation (30 panels chordwise, 30 panels spanwise, Dq ¼ 5 ). In the upwind portion of the rotor, convergence is reached at 7 rotations, after that the differences in Ref. az0(q) are bounded at 0.1. In the downwind part the convergence is less monotonic with high fluctuations due to the intense bwi which differs for each realization. However a good degree of convergence is also reached at 7 rotations with variations in calculated angles of attack around 0.2 . The influence of the spatial discretization is investigated by refining the number of panels on the blade, both in chordwise and spanwise direction. Increasing the number of spanwise panels increases also the number of elements in the wake resulting in a significant increase in computational time; while increasing the number of chordwise panels has a smaller effect. Fig. 5 shows convergence analysis for the chordwise panels. The chordwise discretization seems to be of influence, with the first refinement (from 16 to 30 cpan) being of higher impact than the second (from 30 to 60 cpan), and 30 chordwise panels already giving a good convergence. Fig. 6 shows the effect of a refinement in the number of spanwise panels. A first refinement from 16 to 30 spanwise panels decreases the az0(q) at q ¼ 90 of circa 0.1 (with an increase of computational time of 3). A similar decrease in the calculated angle of attack is achieved by a further refinement from 30 to 60 panels (with an increase of computational time of 8.5). Doubling again the number of spanwise panels leads the az0(q) to higher values. A value of 60 spanwise panels seems to be the most accurate solution (within the cases considered), while reducing to 30 panels lowers the accuracy minimally with a considerable gain in computational time. 3.5. Numerical stability The sensitivity analysis shows a similar non-monotonic behavior for the spanwise spatial or temporal discretizations. When the discretization is refined (doubling the number of spanwise panels or halving the azimuth step), az0(q) first decreases to then increase again (Fig. 7 shows this behavior for the blade at the most upwind position). While the initial decrease is considered to be caused by the reduction of the discretization error, the subsequent increase is attributed to the effect of a numerical instability of the model and analyzed in the present paragraph. The flowfields show localized regions of high values of velocity and velocity gradients. Such regions start to appear around 1.5
556
G. Tescione et al. / Renewable Energy 87 (2016) 552e563
Fig. 3. Effect of the azimuth step ðDqÞ on the angle of attack at the middle of the blade (az0(q)) across the azimuth (q).
diameters downstream the turbine shaft in the center of the wake, with small differences depending on the discretization parameters. The study of the wake geometry highlighted pronounced deformations in the central part of the wake. These are the results of strong close interactions between wake elements and blade panels and mutually between vortex filaments. A higher refined temporal or spatial grid results in a higher probability of such instabilities
which tend to be amplified by the Lagrangian approach and the explicit time scheme. A simulation where the body influence had been excluded led to a very limited improvement in stability, revealing how the vortex filaments have a dominant effect. At the rotor level the discretization error reduction is limited and a refinement of the spatial or temporal grid does not lead to a better solution. At the wake level the solution is stable and rapidly
Fig. 4. Effect of the simulated number of rotations (rot) on the angle of attack at the middle of the blade (az0(q)) across the azimuth (q).
G. Tescione et al. / Renewable Energy 87 (2016) 552e563
557
Fig. 5. Effect of the number of chordwise panel (cpan) on the angle of attack at the middle of the blade (az0(q)) across the azimuth (q).
converged until a certain point, following the motion of the intense structures of the tip vortices, after which high gradients and fluctuations of velocity are present. The stability of the wake has been analyzed by evaluating the absolute value of the velocity (according to Equation (3)) on the median line of the wake (y ¼ 0,z ¼ 0) for each azimuth position of the final simulated rotation, and calculating both the averaged value
and standard deviation. Fig. 8 shows the velocity profile for one of the simulated cases (RB1, Dq ¼ 5 ,span ¼ 30,cpan ¼ 30,rot ¼ 8). The plots show the upwind (A) and downwind (C) blade passages (with high temporal fluctuations due to the position of the blades), the interactions with the vortex filaments inside the rotor area (B) and in the near wake (D), and the point where the vortex filaments of the tip vortices intersect the median line causing high gradients and
Fig. 6. Effect of the number of spanwise panel (span) on the angle of attack at the middle of the blade (az0(q)) across the azimuth (q).
558
G. Tescione et al. / Renewable Energy 87 (2016) 552e563 Spanwise Panels 15
30
60
120 Spanwise Panels
α
90
[deg]
8.20
8.15
8.10
8.05
Δθ
10
Δθ [deg]
5
2
1
Fig. 7. Variation of the perceived angle of attack at the middle of the blade in the most upwind position (az0 90 at the refinement of the spatial (span) and temporal (Dq) discretization.
Fig. 8. Absolute velocity profiles, average and standard deviation on the wake median line (y ¼ 0,z ¼ 0) for the RB1 run (Dq¼5+,span ¼ 30,cpan ¼ 30,rot ¼ 8). Upper plot shows the velocity profiles for the different blade azimuth (V(q), thin gray lines) and the velocity profile averaged across one full rotation (VAVG, thick black line). Lower plot shows the standard deviation (std(V), thick red line) across one full rotation. All velocities are made non-dimensional with the free-stream velocity.
fluctuations of velocity (E). Table 1 presents the results of the stability analysis for all the 10 different runs analyzed in Section 3.4 for the sensitivity analysis. The results are grouped by the 3 parameters defining the discretization (Dq,span,cpan). The limit of stability (XST) has been measured as the location where the standard deviation of the absolute velocity over one rotation starts to peak and is expressed in rotor diameters (D) from the turbine shaft. The table shows how a coarser blade paneling or higher azimuth steps lead to a more stable wake up to 1.5e1.7 diameter downwind of the turbine shaft, while the finer paneling or smaller time step reduce the stable part of the wake down to 1.1e1.3 diameter downwind. More work is needed, and undergoing, to improve the
Table 1 Limit of stability for the different runs, effect of the discretization parameters. Effect of Dq
Effect of span
Effect of cpan
Run
Dq
XST
Run
Span
XST
Run
Cpan
XST
RB5 RB1 RB3 RB6
10+ 5+ 2 1
1.45 1.38 1.24 1.06
RB8 RB1 RB2 RB7
15 30 60 120
1.56 1.38 1.40 1.32
RB10 RB1 RB4 RB9
16 30 60 120
1.71 1.38 1.35 1.12
convergency and stability of the further part of the wake solution.
4. Model validation From the results of the sensitivity and stability analysis (Sections 3.4 and 3.5 respectively), the model used for the validation has a 5 azimuth step, 30 panels chordwise, 30 panels spanwise and a 7 rotations transient, in order to reach a trade-off between convergence, stability and computational cost. The solution is considered to be stable up to 1.4 diameters downwind the turbine shaft in the mid-plane. On the workstation used2 the computational time is circa 30 min. The experimental campaign [21] provided phase-locked 2D velocity fields at the turbine mid-span horizontal plane and 3D velocity fields at 7 vertical planes spanning in the cross-stream direction (y/ D ¼ 0.5, 0.4, 0.2, 0, 0.2, 0.4, 0.5) in the rotor area and till 3 diameters downstream. The reader is referred to [21] for a detailed phenomenological explanation of the nature of the vortical
2 intel i7-950 with 8 cores @ 3 GHz, 12 GB RAM, and NVIDIA® Quadro 4000 GPU with 256 cores, 243GFlops in double precision and 89.6 GB/s memory bandwidth.
G. Tescione et al. / Renewable Energy 87 (2016) 552e563
559
Fig. 9. Contour levels of streamwise velocity (Vx/V∞) in the mid-plane (z ¼ 0) for the experimental (left) and numerical (right) results.
structures in the vawt wake that would help understanding the following analysis, the same terminology and conventions are also used. 4.1. Horizontal plane The first comparison between the simulations and the experimental data is about the velocity field in the horizontal mid-plane. Figs. 9 and 10 show the streamwise and cross-stream velocity components (Vx and Vy respectively) for the horizontal mid-plane (z ¼ 0). The experimental data are shown on the left and the numerical data on the right. The velocities are made non-dimensional by the free stream V∞ and the experimental and numerical contour plots share the same contour levels for an easy comparison. The streamwise velocities (Fig. 9) are very well predicted for the upwind (Dx 0:5) and in-rotor (0:5 Dx 0:5) area with the exclusion of the region just downstream of the turbine shaft (Dx 0; Dy z0). The experimental results show a dominant viscous wake of the rotating cylinder while in the numerical results this is absent since the shaft is not simulated. The numerical solution of the downwind part of the wake (Dx 0:5) presents a good comparison in terms of wake extension for the leeward side (Dy 0) and a lower expansion for the windward side (Dy 0). The latter is partially due to the absence of the rotating shaft wake. Moreover, the numerical results show a more extended region of higher wake deficit (VV∞x ¼ 0:2) as well as higher acceleration of the flow outside
the wake (VV∞x ¼ 1:1). This may be explained with a higher strength of the vortical structures at the boundaries of the wake and it is also partly direct consequence of the absence of any decay of vorticity in the numerical model. It should be noted that the numerical results are instantaneous phase-locked fields, while the experimental results are phase-locked average. The process of averaging introduces a smoothing effect of the near-field inductions of vortical structures due to a certain degree of phase unsteadiness [21]. The cross-stream velocity plots (Fig. 10) show a good quantitative agreement of the contours in all the shown domain. In the rotor area (0:5 Dx 0:5) the vortex elements follow closely the position of the blades wake, as well as the vortical structures downwind the rotor (Dx 0:5) in the leeward side (Dy 0). Some differences are present in the velocity contours and in the positions of the structures in the windward side (Dy 0). 4.2. Vertical planes Figs. 11e13 show the streamwise, cross-stream and vertical velocity components (Vx, Vy, and Vz respectively) for the vertical planes (y/D ¼ 0.5, 0.4, 0.2, 0, 0.2, 0.4, 0.5). As for the previous case, the velocities are made non-dimensional by the free stream V∞ and the same contour levels are used for the experimental (left) and numerical (right) results. Observing the contour plots an overall good agreement can be seen between experimental and numerical results. Contour levels
Fig. 10. Contour levels of cross-stream velocity (Vy/V∞) in the mid-plane (z ¼ 0) for the experimental (left) and numerical (right) results.
560
G. Tescione et al. / Renewable Energy 87 (2016) 552e563
Fig. 11. Contour levels of streamwise velocity (Vx/V∞) in the vertical planes (y/D ¼ 0.5, 0.4, 0.2, 0, 0.2, 0.4, 0.5) for the experimental (left) and numerical (right) results.
are well predicted for the streamwise velocity (Fig. 11) and especially for the cross-stream and vertical velocity (Figs. 12 and 13), with some local discrepancies. Also the dynamics of the tip vortices is correctly captured with their inboard motion in the middle part of the wake, already observed in experiments and simulations [8,14,15,21]. The most noticeable difference is the vertical expansion of the
windward side of the wake (Dy 0) where the experiments showed a pronounced vertical expansion, higher than the one in the leeward side (Dy 0), while the numerical solution shows a reduced expansion comparable to the one on the leeward side. A similar result was observed for the horizontal mid-plane. A second mismatch is in the local induction of the tip vortices. The numerical flowfield shows stronger local induction, observable
G. Tescione et al. / Renewable Energy 87 (2016) 552e563
561
Fig. 12. Contour levels of cross-stream velocity (Vy/V∞) in the vertical planes (y/D ¼ 0.5, 0.4, 0.2, 0, 0.2, 0.4, 0.5) for the experimental (left) and numerical (right) results.
in Fig. 11 in the higher accelerated flow outside the wake (VV∞x ¼ 1:1) or in Fig. 13 in the higher contours associated to the vortical structures. The numerical results suggest stronger released tip vortices which conserve their strength as they travel downstream, while the experimental results have shown a fast decay of the vortical structures as soon as the upwind-released vortex interacts with the downwind-released one [21].
Also the absence of the turbine shaft in the Dy ¼ 0 plane is noticed by the absence of its wake as well as the absence of the two struts wakes in the Dy ¼ 0:5 (windward) and Dy ¼ 0:5 (leeward) planes. The cross-stream and vertical velocity fields give more insight in the wake dynamics not correctly predicted by the numerical model and explain some of the differences found in the horizontal mid-
562
G. Tescione et al. / Renewable Energy 87 (2016) 552e563
Fig. 13. Contour levels of vertical velocity (Vz/V∞) in the vertical planes (y/D ¼ 0.5, 0.4, 0.2, 0, 0.2, 0.4, 0.5) for the experimental (left) and numerical (right) results.
plane. The overall behavior of the wake is well captured. The horizontal expansion in the cross-stream velocity contours (yellow regions of positive values in the windward side and blue regions of negative values in the leeward side) and the vertical motion in the vertical velocity contours (yellow regions of positive value in the windward and leeward side and blue regions of negative values in the central part of the wake) are similar. Some regions present significant quantitative discrepancies: the windward side of the
wake downstream the rotor (Dy 0 and Dx 0:5) shows significant lower cross-stream (Fig. 12) and vertical (Fig. 13) velocities. The reason is to be found again in the higher strength of the vortical structures, opposed to the experimental case (Rec ¼ 170000), causing a stronger local induction velocity field. The observed reduced expansion of the wake of the windward side both in the horizontal and in the vertical planes confirms this. Also the vertical inboard motion of the central part of the wake observed in the
G. Tescione et al. / Renewable Energy 87 (2016) 552e563
experiment is reduced in the numerical analysis. A last observation concerns the numerical stability of the wake in the vertical planes. In Section 3.5 the analysis focused on the mid-plane, showing a wake numerically stable until at least Dx ¼ 1:4. Figs. 12 and 13 show signs of numerical instabilities (represented by high local values of induced velocities and velocity gradients) developing around the tip vortices at smaller downwind positions. Such instabilities, though of a rather localized nature, are also partially responsible of the local differences between numerical and experimental results observed in the present analysis.
563
comparison between a numerical model without diffusion scheme and an experimental model run at a low Reynolds number. Moreover the absence of the struts in the simulation may be responsible for an initial overshoot of the tip vortex strength and the absence of the turbine shaft for a part of the wake asymmetry. It is to believe that simulations of higher Reynolds number cases, comprising the struts and the turbine shaft will see their accuracy increased, making the model a valuable tool for the analysis of multi-MW VAWTs. References
5. Conclusions A free vortex wake model has been verified and validated with experimental measurements of a vawt wake. The model has been analyzed in terms of sensitivity of the output to the simulation parameters. For the flow inside the rotor (for loads evaluations of a single turbine), a medium discretization (30 chordwise panels, 30 spanwise panels, 5 of azimuth step) gives a converged solution in terms of angles of attack experienced by the blade. The initial transient to reach the periodic converged solution, for the present application, is around 7 rotations. On a normal workstation equipped with a gpu such simulation leads to a computational time of 30 min. Any further grid refinement, by doubling the blade panels or halving the azimuth step, or any added simulated rotation leads to a marginal change in the angle of attack. The simulation has also been analyzed in terms of numerical stability. The tip vortex evolution represents the most critical factor for the stability of the model. Given the high local induction of these vortical structures and the Lagrangian formulation of the model, the wake is subject to strong local deformations leading to high velocity gradients. A second-order time scheme and the desingularization of the viscous vortex core help controlling such effect. Moreover the grid refinement, while reducing the discretization error by a limited amount, decrease the stability of the model. The effects of such instabilities are rather local and do not travel upwind, so the rotor and the near-wake flows are less affected. With the parameters used in the present analysis, the numerical solution in the wake is considered stable up to 1.4 diameters downstream the turbine shaft in the symmetry plane. The numerical results have been compared to the measurements of a previous experimental activity on a 2-bladed hvawt, running at a tsr of 4.5 and a chord-based Reynolds number of 170000. The model is capable of capturing all the important dynamics in the wake of vawt and give also good quantitative evaluation of the flowfield. Comparisons show very good prediction of the flowfield upwind and inside the rotor area. The dynamics of the vortical structures in the horizontal plane is very well captured as well as the vertical motion of the tip vortices with the already observed central contraction and lateral expansion. The results in the near wake (0:5 Dx 1:5) show good prediction of the wake flow in the leeward side but tend to underestimate the wake horizontal and vertical expansion in the windward side compared to the experimental results, decreasing the high asymmetry observed in the experiments. Also the local induction of the vortical structures in the far part of the wake (both of the horizontal plane and of the tip vortices) is higher in the numerical results. The reasons for such discrepancies are to be found in the
rimentale et Nume rique du De crochage Dynamique [1] P.L. Beaudet, Etude Expe (Ph.D. thesis), Universit de sur une Eolienne a Axe Vertical de Forte Solidite Poitiers, 2014. [2] M. Borg, M. Collu, F.P. Brennan, Offshore floating vertical axis wind turbines: advantages, disadvantages, and dynamics modelling state of art, Mar. Offshore Renew. Energy (2012) 33e46. [3] M. Borg, A. Shires, M. Collu, Offshore floating vertical axis wind turbines, dynamics modelling state of the art, Part I Aerodyn. Renew. Sustain. Energy Rev. 39 (2014) 1214e1225. [4] P. Deglaire, Analytical Aerodynamic Simulation Tools for Vertical Axis Wind Turbines (Ph.D. thesis), Uppsala University, 2010. [5] K.R. Dixon, The Near Wake Structure of a Vertical Axis Wind Turbine (Master's thesis), Delft University of Technology, 2008. [6] J.L. Hess, Panel methods in computational fluid dynamics, Annu. Rev. Fluid Mech. 22 (1990) 255e274. [7] J.L. Hess, A.M.O. Smith, Calculation of potential flow about arbitrary bodies, Prog. Aeronaut. Sci. 8 (1967) 1e138. [8] C. Hofemann, C.J. Sim~ ao Ferreira, K. Dixon, G.J. van Bussel, G. van Kuik, F. Scarano, 3D stereo piv study of tip vortex evolution on a vawt, in: Proceeding of EWEC, Brussels, 2008. [9] J. Katz, A. Plotkin, Low-speed Aerodynamics: from Wing Theory to Panel Methods, McGraw-Hill, New York, 1991. [10] L. Morino, C.C. Kuo, Subsonic potential aerodynamics for complex configurations: a general theory, AAIA J. 12 (2) (1974) 191e197. [11] J.C. Murray, M. Barone, The development of cactus: a wind and marine turbine performance simulation code, in: 49th AIAA Aerospace Sciences Meeting, Orlando, FL, 2011. [12] P.J. Musgrove, Wind energy conversion: recent progress and future prospects, Sol. Wind Technol. 4 (1) (1987) 37e49. [13] S. Peace, Another approach to wind, Mech. Eng. 126 (6) (2004) 28e31. [14] F. Scheurich, T.M. Fletcher, R.E. Brown, Simulating the aerodynamic performance and wake dynamics of a vertical-axis wind turbine, Wind Energy 14 (2) (2011) 159e177. ~o Ferreira, The Near Wake of the VAWT, 2D and 3D Views of the [15] C.J. Sima VAWT Aerodynamics (Ph.D. thesis), Delft University of Technology, 2009. [16] C.J. Sim~ ao Ferreira, H. Aagard Madsen, M. Barone, B. Roscher, P. Deglaire, I. Arduin, Comparison of aerodynamic models for vertical axis wind turbines, J. Phys. Conf. Ser. 524 (1) (2014). ~o Ferreira, F. Scheurich, Demonstrating that power and instantaneous [17] C.J. Sima loads are decoupled in a vertical-axis wind turbine, Wind Energy 17 (2014) 385e396. [18] M.J. Stock, A. Gharakhani, Toward efficient gpu-accelerated n-body simulations, in: 46th AIAA Aerospace Sciences Meeting and Exhibit, January 2008. [19] J.H. Strickland, B.T. Webster, T. Nguyen, A vortex model of the darrieus turbine: an analytical and experimental study, J. Fluids Eng. 101 (1979) 500e505. [20] H.J. Sutherland, D.E. Berg, T.D. Ashwill, A Retrospective of Vawt Technology, Technical report, Sandia National Laboratories, January 2012. [21] G. Tescione, D. Ragni, C. He, C.J. Sim~ ao Ferreira, G.J.W. van Bussel, Near wake flow analysis of a vertical axis wind turbine by stereoscopic particle image velocimetry, Renew. Energy 70 (2014) 47e61. [22] L. Vita, U.S. Paulsen, T.F. Pedersen, H.A. Madsen, F. Rasmussen, Deep wind: a novel floating wind turbine concept, Wind Int. 6 (4) (2010) 29e31. [23] Vv.Aa, Wind in Our Sails. The Coming of Europe's Offshore Wind Energy Industry, Technical report, European Wind Energy Association (EWEA), November 2011. [24] Vv.Aa, Deep Water. The Next Step for Offshore Wind Energy, Technical report, European Wind Energy Association (EWEA), July 2013. [25] D.J. Willis, An Unsteady, Accelerated, High Order Panel Method with Vortex Particle Wakes (Ph.D. thesis), Massachusetts Institute of Technology, 2006. ~o Ferreira, Wake modelling of a vawt in [26] A. Zanon, P. Giannattasio, C.J. Sima dynamic stall: impact on the prediction of flow and induction fields, Wind Energy 18 (11) (November 2015) 1855e1874, http://dx.doi.org/10.1002/ we.1793.