The numerical analysis of wind turbine airfoils at high angles of attack

Abstract

Despite decades of research, the accurate numerical simulation of severely separated flows is still one of the major problems in computational fluid dynamics. In this paper, a viscous-coupled 3D panel method is proposed for the aerodynamic analysis of wind turbine airfoils at high angles of attack. The Hess-Smith panel method is adopted for inviscid calculations, and an empirically based boundary layer analysis is performed in order to determine the separation point. The separated thick wake is then modelled as an extension of the surface geometry along which a constant pressure distribution is assumed. The wake geometry is determined iteratively, and an outer iterative loop is run to update the location of the separation point. The validity of the current method's results is assessed by comparison with experimental and numerical results for several high thickness wind turbine airfoils, namely NACA 63-430, FFA-W3-301, FFA-W3-241, and DU 91-W2-250. At low angles of attack, pressure data predicted by the current method show excellent agreement with the experimental data, as well as the referenced numerical data. At higher angles of attack, the current method shows reasonable agreement with the experimental data, while the referenced numerical data significantly overestimate the -C ^sub p^ distribution along the suction surface.

http://crossmark.crossref.org/dialog/?doi=10.1007/s40095-015-0197-6&domain=pdf

Web End = http://crossmark.crossref.org/dialog/?doi=10.1007/s40095-015-0197-6&domain=pdf

Web End = http://crossmark.crossref.org/dialog/?doi=10.1007/s40095-015-0197-6&domain=pdf

Web End = http://crossmark.crossref.org/dialog/?doi=10.1007/s40095-015-0197-6&domain=pdf

Web End = http://crossmark.crossref.org/dialog/?doi=10.1007/s40095-015-0197-6&domain=pdf

Web End = Int J Energy Environ Eng (2016) 7:112 DOI 10.1007/s40095-015-0197-6

http://crossmark.crossref.org/dialog/?doi=10.1007/s40095-015-0197-6&domain=pdf

Web End = http://crossmark.crossref.org/dialog/?doi=10.1007/s40095-015-0197-6&domain=pdf

Web End = http://crossmark.crossref.org/dialog/?doi=10.1007/s40095-015-0197-6&domain=pdf

Web End = http://crossmark.crossref.org/dialog/?doi=10.1007/s40095-015-0197-6&domain=pdf

Web End = ORIGINAL RESEARCH

http://crossmark.crossref.org/dialog/?doi=10.1007/s40095-015-0197-6&domain=pdf

Web End = The numerical analysis of wind turbine airfoils at high angles of attack

Bryan Nelson1 Jen-Shiang Kouh1

Received: 11 August 2015 / Accepted: 7 November 2015 / Published online: 30 November 2015 The Author(s) 2015. This article is published with open access at Springerlink.com

Abstract Despite decades of research, the accurate numerical simulation of severely separated ows is still one of the major problems in computational uid dynamics. In this paper, a viscous-coupled 3D panel method is proposed for the aerodynamic analysis of wind turbine airfoils at high angles of attack. The HessSmith panel method is adopted for inviscid calculations, and an empirically based boundary layer analysis is performed in order to determine the separation point. The separated thick wake is then modelled as an extension of the surface geometry along which a constant pressure distribution is assumed. The wake geometry is determined iteratively, and an outer iterative loop is run to update the location of the separation point. The validity of the current methods results is assessed by comparison with experimental and numerical results for several high thickness wind turbine airfoils, namely NACA 63-430, FFA-W3-301, FFA-W3-241, and DU 91-W2-250. At low angles of attack, pressure data predicted by the current method show excellent agreement with the experimental data, as well as the referenced numerical data. At higher angles of attack, the current method shows reasonable agreement with the experimental data, while the referenced numerical data signicantly overestimate the -Cp distribution along the suction surface.

Keywords Wind turbine Panel method Boundary layer

analysis Separated wake

Introduction

Panel methods are numerical techniques which solve potential ow problems around complex geometries by employing singularity distributions, or panels, to represent an objects surface and wake. Although panel methods are generally limited to the solution of potential ow problems(i.e. inviscid, incompressible and irrotational ow), there are several boundary layer models which extend the scope of analysable problems to include those involving viscous ow. Furthermore, the fast computation times of panel methods, in the order of seconds, make them an excellent CFD tool at the preliminary design stage. On the other hand, Reynolds-averaged NavierStokes (RANS) solvers with their long computation times, typically in the order of hours, are less suitable at the preliminary design stage, and as such are usually relegated to post-design analysis of ow problems. A further problem with RANS solvers is the signicant amount of time that is generally required for mesh generation, while in panel methods, discretisation is conned to the object surface, allowing for fast, convenient modication of geometric parameters. The panel method adopted for this research is that pioneered by Hess and Smith [1], which discretises the surface geometry into a number of quadrilateral panels comprising constant strength source/doublet distributions, and the wake is modelled as a sheet of doublet panels extending downstream. An open-source version of this 3D panel method, released by Filkovic [2], who based his program on the Fortran code provided by Katz and Plotkin [3], is currently being adapted for the analysis of a full, 3D, rotating wind turbine blade.

The aerodynamic analysis of wind turbine blades, however, presents a number of challenges. In addition to the difculties presented by the rotation of the rotor,

& Bryan Nelson [email protected]

1 Department of Engineering Science and Ocean Engineering, National Taiwan University, No. 1, Sec. 4, Roosevelt Rd.,, Taipei 10617, Taiwan

123

2 Int J Energy Environ Eng (2016) 7:112

uctuations in wind speed and direction, rotortower interaction, and aeroelasticity, to name a few, the CFD analysis of wind turbine blades is further complicated by high bending moments close to the blade root, which necessitate a sacrice of aerodynamic efciency as the thickness-to-chord ratio is increased to improve structural integrity. A direct result is that, as the blade cross section becomes more circular approaching the root, the ow eld deviates from that predicted by potential ow and separation becomes unavoidable, thus requiring a solution of the boundary layer. For RANS solvers, this necessitates a very ne mesh resolution (y? \ 1) near the wall, which is very computationally expensive. For potential ow solvers, this boundary layer solution is typically based on empirical viscous corrections.

In inviscid ow, there is no boundary layer and the ow is everywhere tangent to the body surface, which may therefore be represented by a streamline. In viscous ow, however, the no-slip condition states that the uid at the body surface has zero relative velocity, resulting in a reduction in the ow rate around the body. Two methods of modelling this reduction in ow rate, rst described by Lighthill [4], are by a normal displacement of the inviscid surface streamline by a small distance, known as the boundary layer displacement thickness, d*, or by the addition of a transpiration velocity at the real body surface, which approximates the effective displacement body surface. For ows with more severe separation, the free shear layers bounding the separated wake are modelled as a pair of free vortex sheets [5], the geometries of which

are iteratively determined according to the induced velocities. This is commonly referred to as the double wake method.

The current method is based on the assumption that there is a constant pressure distribution along the section of the airfoils surface between the separation point and the trailing edge, henceforth referred to as the separation region, and, furthermore, that this constant pressure distribution extends throughout the separated wake and may therefore be applied along the streamlines bounding the wake, henceforth referred to as the separation streamlines. All that remains, then, is to determine, iteratively, the conguration of the separation streamlines so as to produce a constant pressure distribution beyond the separation point. In essence, the separated wake is treated as an extension of the body surface, and the separation region is disregarded. This method has been applied to the analysis of several wind turbine airfoils, namely NACA 63-430, FFA-W3-301, FFA-W3-241, and DU 91-W2-250 (Fig. 1), at various angles of attack between 6 and 21. Due to their high thickness-to-chord ratios (0.240.3), these airfoils are commonly used at the mid-span and inboard stations of wind turbine blades. Results have been veried against wind tunnel test data taken from studies performed at the VELUX wind tunnel in stbirk, Denmark [6], and at TU Delfts low-speed low-turbulence wind tunnel [7]. These experimental results, as well as corresponding numerical results, have been included in a comprehensive catalogue of airfoil data compiled by Bertagnolio et al. [8].

Fig. 1 Wind turbine airfoils. a NACA 63-430, b FFA-W3-301, c FFA-W3-241, d DU 91-W2-250

123

Int J Energy Environ Eng (2016) 7:112 3

Fig. 2 Body and wake panels, and collocation points

Fig. 3 Local coordinate system and eld point P

(Fig. 2), and the perturbation potential at each collocation point is calculated by summing the inuences of all body and wake panels:

XNPk11 4p

Numerical formulation

Potential ow

Potential theory is based on the assumption of an irrotational ow of an ideal uid (i.e. inviscid and incompressible, expressed by the continuity equation, r u

Zbody ln

r

* dS X

NW

l1

1 r

1 4p

Zwake ln

*

r

1 r

dS

X

NP

k1

1 4p

Zbody

r r

dS 0: 5

The integrals are evaluated by dening local coordinate systems for each panel, as shown in Fig. 3, such that the boundary integral equations to be solved for a constant strength source distribution and a constant strength doublet distribution are, respectively:

/x; y; z

r 4p

* 0).

* may be expressed as the gradient of some scalar r/. This scalar is called the

velocity potential, and substituting this into the continuity equation gives us Laplaces equation,

r2/ 0 1

which is the governing equation for potential ow. Two solutions to Laplaces equation are:

/

r4pr 2 and

/

l 4p

From irrotationality, where r u

* 0, u

ZS

dSr 6

and

/x; y; z

l 4p

ZS

zdSr3 ; 7

where r

x n

3

which are, respectively, the potentials for a source of strength r and doublet of strength l. The panel method takes advantage of these solutions by using their induced velocity eld to simulate the ow eld around arbitrary solid bodies. The zero normal velocity boundary condition is, in this case, indirectly imposed by setting the total potential equal to the free stream velocity potential, which is only valid for

r n

2z2

q :

The NP source strengths in (5) are known from boundary condition (4), and the remaining unknowns, namely, the

NP ? NW doublet strengths in the body and wake surface panels, are solved by imposing the Morino Kutta condition, which states that the strength of the wake doublet is equal to the difference between the adjacent upper and lower body surface doublets at the trailing edge. By expressing the wake doublets in terms of the unknown surface doublets, the number of unknowns is reduced to NP, and it is now possible to solve the NP linear algebraic equations for lk.

The induced velocity eld on the objects surface (i.e. the induced velocities at each of the panels collocation points) is determined by simply differentiating the doublet strengths with respect to the local coordinates (Eq. 8), and is in turn used to solve the pressure eld (Eq. 9):

o on

1 r

2 y n

* U1: 4

This is called the Dirichlet boundary condition. The geometry may now be discretised into NP surface panels comprising source and doublet distributions and NW wake panels containing doublet distributions. The boundary condition is specied on each panel at a collocation point

123

4 Int J Energy Environ Eng (2016) 7:112

u

where Ue is the ow speed outside of the boundary layer, as predicted by inviscid analysis, and sw is the shear stress at the body surface (i.e. at y = 0).

Combining and manipulating the boundary layer equations, and substituting in the parameters dened above, we are able to obtain von Krmns momentum integral equation:dhdx

hUe 2 H

* locali

olioxlocali 8

Cp;i 1

2i

U21

: 9

The boundary layer

For a steady, 2-D, incompressible ow, the NavierStokes equations, which are the equations of motion for viscous ow, reduce to the following form:

q u ou

ox v

dUe dx

ou oy

opox l

10

o2u

ox2

o2u

oy2

: 11

Within the boundary layer, these equations may be simplied by an order of magnitude analysis to obtain Eqs. (12) and (13), which, together with the 2D continuity equation (Eq. 14), comprise the boundary layer equations.

q u ou

ox v

q u ov

ox v

ov oy

opoy l

o2v

ox2

o2v

oy2

1

2 cf 19

where the shape factor, H d h, has been introduced. By

supplementing Eq. (19) with additional relations, it is possible to determine the point at which the ow separates from the body surface. The methods adopted herein are Thwaites method for laminar boundary layers [9] and Heads method for turbulent boundary layers [10].

Thwaites method

Thwaites introduced two new parameters, l and k, shown below, and rewrote the momentum integral Eq. (19) in terms of these new parameters (Eq. 22)

l

Ueh

m

ou oy

opox l

o2u

oy2 12

opoy 0 13

ou ox

cf

2 20

k

h2 m

dUedx 21

ovoy 0 14

Before going any further, it is rst necessary to introduce a few parameters which play an important role in the description of the boundary layer. As mentioned, the boundary layer displacement thickness, d*,

is dened as the distance by which the inviscid surface streamline must be displaced so as to account for the reduction in ow rate in the boundary layer (Eq. 15). Similarly, the momentum thickness, h, is dened as the distance by which the inviscid streamline must be displaced so as to account for the reduction in total momentum (Eq. 16). Finally, the shear stress in the boundary layer and the skin friction coefcient are given by Eqs. (17) and (18),

d Z

dh2dx 2l k 2 H

22 By plotting l and H against k, Thwaites found that the right hand side of Eq. (22) is very well approximated by:

2l k 2 H

0:45 6k 23 Substituting (23) into (22), and solving for h2, we get:

h2 0:45

m U6e

Ue m

Z

x

0 1

1 dy 15

u Ue

0 U5edx 24

which may, in turn, be used to calculate k. Thwaites found that laminar separation occurs at k \ -0.09. Due to the difculty in modelling laminar separation bubbles, it is assumed that at the laminar separation point the ow immediately reattaches to the surface, and the laminar separation point is taken as the point of transition from a laminar to a turbulent boundary layer.

Heads method

Heads method considers the entrainment of the free stream ow into the turbulent boundary layer, and denes entrainment velocity, E, as the rate of increase in volu-metric ow rate, Q, within the boundary layer:

h Z

uUe 1

1 dy 16

u Ue

0

s l

ouoy 17

cf

sw12 qU2e

18

123

Int J Energy Environ Eng (2016) 7:112 5

dQdx ; 25

where

Q x

Z

d

E

udy: 26

Substituting in displacement thickness (15) and introducing a turbulent shape factor, H1 d dh, the entrainment

velocity may be rewritten:

E

0

!; 32

*

i

*

i

u

n

u

t

*

i

*

i

ddx UehH1: 27

The dimensionless entrainment velocity, E/Ue, is

determined using the following empirical formulae: 1

Ue

where n

* i and t

ddx UehH1 0:0306 H1 3

0:617; 28

where H1 3:3 0:8234 H 1:1

1:287 for H 1:6 and H1 3:3 1:5501 H 0:6778

3:064 for H [ 1:6

29

Finally, Head adopts the LudwiegTillman skin-friction law to solve for cf:

cf 0:246 10 0:678HRe 0:268h; 30 where

Reh

Ueh

m 31 and a 2nd order RungeKutta method is employed to solve ordinary differential Eqs. (19) and (28) for h and H1.

Turbulent separation is assumed to occur when H1 \ 3.3, as per Wauquiez [11].

The separated wake

The double wake method

Once the turbulent separation point has been determined, the next step is to model the separated wake. In the 1970s, Maskew and Dvorak [5] developed a method for modelling the ow around airfoils at high angles of attack, where the separated wake is treated as a region of constant total pressure bounded by a pair of free shear layers, which are modelled by vortex panels. Consequently, this approach is often referred to as the double wake method.

The wake between the two trailing shear layers is a region with low vorticity and insignicant viscous stresses, and is therefore taken to be a potential ow region. As the separated wake shape is not known a priori, an initial wake geometry is specied according to an empirical relationship between the

airfoils thickness-to-chord ratio (a wake length factor is introduced, dened as WF : 0.081 9 c/t ? 1.1) and the extent of separation (represented by the wake height, illustrated in Fig. 4). The induced velocities on the wake panels are computed, and the panel positions are then iteratively adjusted until the normal velocity on each wake panel is equal to zero. According to Robinson [12], the panels are rotated according to the following equation:

Dhi arctan

* i are, respectively, the normal and tangential vectors for the ith wake panel.

As mentioned, the induced velocity eld on the objects surface is determined by simply differentiating the doublet strengths with respect to the local coordinates, as per Eq. (8). On the other hand, the calculation of off-body induced velocities, as required by Maskews wake relaxation algorithm, is a more complicated and laborious process, which is not necessary within the current method, and as such is not covered here. (The reader may refer to Katz and Plotkin [3],p. 284 and p. 287 for velocities induced by, respectively, constant strength source and doublet distributions.)

A further complication that arises with this method is the need to offset the source distributions along the airfoil area between the two free shear layers by adding a pressure jump to the panel stations situated on this separation region. The pressure coefcient at the ith panel is now given by:

Cp;i 1

u2i

U21

:

kH DH

; 33

where DH is the increase in total pressure inside the separation bubble, being equal to zero all over the airfoil except in the separated wake region. The calculation of DH is covered by Ramos Garca [13].

The thick wake method (current method)

The displacement body model, where the inviscid surface streamline is offset from the actual body surface by the boundary layer displacement thickness d*, is already an established method of approximating the real viscous ow

Fig. 4 Separated wake description

1

2 qU21

123

6 Int J Energy Environ Eng (2016) 7:112

eld around an object. The current method extends this idea to include not just the boundary layer, but also the separated wake. As per the above method, the separated wake is treated as a region of constant total pressure, such that the pressure along the separation region is equal to that along the separation streamlines, which is set equal to the pressure at the separation point. The aforementioned problems of calculating the off-body induced velocity eld and the pressure jump in the separated wake are thus avoided by disregarding the separation region and treating the separation streamlines bounding the separated wake as an extension of the displacement body streamline (as illustrated in Fig. 5). This means that the wake surface may be modelled as per the airfoil surface, namely, by employing source distributions to simulate the thickness of the separated wake, hence the name: thick wake method. To the authors knowledge, this method of modelling separated ow has not been presented in the literature.

The initial wake geometry is specied as per Fig. 4. The wake panel positions are then adjusted according to the local pressure distribution [easily obtained via Eqs. (8) and (9)], where the desired result is a constant pressure distribution beyond the separation point. As the pressure values are known only at the panel collocation points, it is necessary to rst linearly interpolate Cp,i to determine the pressure values at the panel vertices. The wake shape may now be adjusted by moving the ith panel vertex vertically from position zin to

zin?1, as calculated by the nth iteration of the secant method:

zn1i zni

zni zn 1i DCnp;i DCn 1p;i

DCnp;i; 34

where DCp,i = Cp,sep - Cp,i, i.e. the difference between the pressure coefcient at the separation point and that at the ith wake panel vertex. Once the desired pressure distribution has been attained, an outer iterative loop is run to update the location of the separation point. Convergence is obtained when the predicted separation point is the same for two successive runs of this outer loop, i.e. when xnout1sep xnoutsep. This interaction between the inviscid and

viscous solvers is known as direct coupling, and is addressed in the following section.

Viscousinviscid coupling

As mentioned above, direct coupling of the inviscid and viscous solvers refers to the case where the inviscid velocity eld Ue,i is adopted as an input parameter for the boundary layer solution, and a viscous correction, in the form of either a boundary layer displacement or transpiration velocity, is in turn incorporated into the inviscid solver. This is certainly the most straightforward and intuitive coupling method, but the solution usually breaks down at separation due to the so-called Goldstein singularity. This is caused by a breakdown in the inviscid-to-viscous hierarchy, where the inviscid ow is no longer dominant.

Drela [14] was able to overcome these hierarchal problems by treating both the velocity and the boundary layer thickness as unknowns, and then merging the viscous and inviscid ow equations into one big system of nonlinear equations, which is solved simultaneously by employing the NewtonRaphson method. This fully simultaneous coupling, while robust, is computationally expensive due to full matrix inversion. More recently, fully simultaneous coupled viscousinviscid solvers employing the double wake method have been applied to yacht sail systems [15], the unsteady case of pitching airfoils in stall conditions [16], and the quasi-3D (2D panel method with rotational effects from Coriolis and centrifugal forces) cases of, initially, rotating airfoils [17] and, later, full wind turbine rotors [18].

To reduce the computational cost of the fully simultaneous coupling method, Veldman [19] introduced the concept of the interaction law, which gives an approximation of the outer inviscid ow. The interaction law and the boundary layer equations are solved simultaneously, and the solution is then input into the true inviscid solver, as per direct coupling. Most of the recent work on this quasi-simultaneous method has been focused on investigating and deriving suitable interaction laws [2023].

Despite the aforementioned convergence problems with direct coupling, the current method has been found to satisfactorily simulate severely separated ows around a number of thick airfoils. By modelling the separated wake as an extension of the displacement body, the thick wake model is able to avoid the Goldstein singularity, allowing for a fast, robust CFD tool, for which the computational expense of simultaneously solving the inviscid and boundary layer equations is no longer necessary. Furthermore, the simplicity and intuitiveness of the current method mean that it may easily be modied or adapted for the analysis of a wide range of uid dynamics problems. This should facilitate the achievement of the long-term objective of this study, which is to develop a 3D panel method for the analysis of a full, 3D, rotating wind turbine blade.

Fig. 5 Illustration of separation streamlines as extension of displacement body streamline

123

Int J Energy Environ Eng (2016) 7:112 7

Table 1 Summary of results of airfoil analyses

Airfoil Angle of attack Separation region1 xsepc

Error (XF)

minC

p;XF

minCp;exp 1

Error (ES)

minC

p;ES

minCp;exp 1

Error (CM)

minC

p;CM

minCp;exp 1

NACA 63-430 6.5 0.38 0.08 0.14 010.6 0.46 0.19 0.26 -0.0716.3 0.6 0.43 0.52 0.0920.6 0.67 0.53 0.51 0.15 FFA-W3-301 6.5 0.19 0.12 0.1 0.059.7 0.39 0.13 0.12 0.0816.3 0.67 0.82 0.91 0.1520.2 0.76 0.95 -0.51 0.24 FFA-W3-241 6.7 0.08 0.03 0.01 0.019.9 0.25 0.07 0.04 -0.0617.9 0.67 0.49 0.51 0.0421.4 0.79 0.5 0.42 0.06 DU 91-W2-250 7.7 0.1 0.11 0.04 09.7 0.19 0.14 0.1 -0.0211.7 0.3 0.22 0.25 0.1215.2 0.56 0.37 0.12 0.12

Results

In order to assess the validity of the current methods results for the case of a rectangular wing in uniform ow, computed pressure distribution data at mid-span are compared with experimental and numerical data taken from Bertagnolios Wind Turbine Airfoil Catalogue [8]. The airfoils analysed are NACA 63-430, FFA-W3-301, FFAW3-241, and DU 91-W2-250, which, due to their high thickness-to-chord ratios (0.240.3), are commonly used at the mid-span and inboard stations of wind turbine blades. Ref. [8] includes wind tunnel results from the VELUX wind tunnel in stbirk, Denmark (NACA and FFA airfoils), and TU Delfts low-speed low-turbulence wind tunnel (DU airfoil), as well as numerical results predicted by XFOIL [24], a very popular viscous-coupled panel method, and EllipSys2D [25], a 2D RANS solver. The results for these airfoils at various angles of attack between6.5 and 21.4 are summarised in Table 1 and Figs. 6, 7, 8 and 9.

NACA 63-430

Due to the high thickness ratio of this airfoil (t/c = 0.3), even at low angles of attack, the ow is still observed to separate quite severely. At an angle of attack of a = 6.5

(Fig. 6a), the separation point predicted by the current method is located at x/c = 0.62, or in other words, the separation region covers 38 % chordwise of the upper surface. As the angle of attack increases, so does the

severity of the separation, with a predicted separation region of 46 % chord at a = 10.6 (Fig. 6b), 60 % chord at a = 16.3 (Fig. 6c), and 67 % chord at a = 20.6 (Fig. 6d).

For all tested angles-of-attack, the referenced numerical data are found to overpredict the -Cp distribution on the suction surface, and at higher angles of attack, there is very poor agreement between the XFOIL and EllipSys2D data and the experimental data, with the -Cp suction peak overpredicted by more than 45 % at a = 16.3, and by more than 50 % at a = 20.6. Comparatively, the current method is seen to correlate far more closely with the experimental data, with excellent agreement at a = 6.5, and a slightly underpredicted suction peak at a = 10.6. Even at a = 16.3 and a = 20.6, the -Cp suction peaks are overpredicted by just 8 and 15 %, respectively, which is clearly a vast improvement over the referenced numerical data.

FFA-W3-301

Although the thickness of the FFA-W3-301 airfoil is equal to that of the previous airfoil (t/c = 0.3), the separation at low angles of attack is far less severe, with a predicted separation region of just 19 % chord at a = 6.5 (Fig. 7a).

At higher angles of attack, however, the separation exceeds that of the previous airfoil, with a predicted separation region of 67 % chord at a = 16.3 (Fig. 7c), and as much as 76 % at a = 20.2 (Fig. 7d).

The effect of this extensive ow separation, and the associated difculty in modelling such separation, is

123

8 Int J Energy Environ Eng (2016) 7:112

Fig. 6 Comparison of pressure distribution data for NACA 63-430. a a = 6.5, b a = 10.6, c a = 16.3, d a = 20.6

evident from an analysis of the -Cp plots. Once again, the referenced numerical data are found to overpredict the -Cp distribution on the suction surface, and at higher angles of attack, there is, again, very poor agreement between the XFOIL and EllipSys2D data and the experimental data, with both models overpredicting the -Cp suction peak by more than 82 % at a = 16.3, and with XFOIL overpredicting the -Cp suction peak by more than 95 %, while EllipSys2D is seen to underpredict the -Cp suction peak by more than 50 %, at a = 20.2. Comparatively, the current method is seen to correlate far more closely with the experimental data, with marginally overpredicted results at a = 6.5, and a slightly underpredicted suction peak at a = 9.7. At a = 16.3 and a = 20.6, the -Cp suction peaks are overpredicted by 15 and 24 %, respectively, which is still a vast improvement over the referenced numerical data.

FFA-W3-241

The FFA-W3-241 airfoil is considerably thinner than the previous two airfoils (t/c = 0.24), and so ow separation should be less severe than that of the previous cases. This is certainly true at low angles of attack, with a predicted separation region of just 7 % chord at a = 6.7 (Fig. 8a), and 25 % chord at a = 9.9 (Fig. 8b), but at higher angles of attack, the extent of the separation is very similar to that of the FFA-W3-301 airfoil, with a predicted separation region of 67 % chord at a = 17.9 (Fig. 8c), and 79 % at a = 21.4 (Fig. 8d).

Due to the low separation at low angles of attack, the referenced numerical data are found to predict the -Cp distribution fairly accurately. However, as the angle of attack increases, and with it the extent of the separation region, the referenced numerical data are seen to deviate substantially from the experimental data, overpredicting

123

Int J Energy Environ Eng (2016) 7:112 9

Fig. 7 Comparison of pressure distribution data for FFA-W3-301. a a = 6.5, b a = 9.7, c a = 16.3, d a = 20.2

the -Cp peak by 51 % at a = 17.9, and by 54 % at a = 21.4. On the other hand, the current method shows excellent agreement with the experimental data for all tested angles of attack, with a maximum error of just 7 % at a = 21.4.

DU 91-W2-250

The nal airfoil analysed is the DU 91-W2-250 airfoil, for which t/c = 0.25. Like the FFA-W3-241 airfoil, the ow at low angles of attack is almost fully attached, with a separation region of just 10 % chord at a = 7.7 (Fig. 9a), and 19 % chord at a = 9.7 (Fig. 9b). Even at higher angles of attack, the extent of the separation is less severe than that of the other three airfoils, with a predicted separation region of 30 % chord at a = 11.7 (Fig. 9c), and 56 % at a = 15.2 (Fig. 9d).

Due to the low separation at low angles of attack, the referenced numerical data are found to predict the -Cp distribution fairly accurately, with the -Cp peak overpredicted by around 10 %. Bear in mind that the DU 91-W2-250 airfoil produces a pronounced -Cp peak, even at low angles of attack, and that beside the slightly overpredicted -Cp peak, there is excellent agreement between the numerical and experimental results. There is also a slightly improved correlation between the referenced numerical data and experimental data at higher angles of attack, with a maximum deviation of 37 % at a = 15.2. The current method once again shows excellent agreement with the experimental data at low angles of attack, and only a slight overprediction of the -Cp suction peak at higher angles of attack, with a maximum deviation of 14 % at a = 15.2.

All results are summarised in Table 1, where the abbreviations XF, ES, and CM, respectively, denote

123

10 Int J Energy Environ Eng (2016) 7:112

Fig. 8 Comparison of pressure distribution data for FFA-W3-241. a a = 6.7, b a = 9.9, c a = 17.9, d a = 21.4

results predicted by XFOIL, EllipSys2D, and the current method.

Discussion

For mildly separated ows, the -Cp plots (Figs. 6, 7, 8, 9) show excellent agreement between all data sets. In the presence of severe ow separation, however, it is evident that the referenced numerical methods, XFOIL and Ellipsys2D, are not able to accurately predict the location of the separation point, and consequently are found to drastically overpredict the pressure distributions along the airfoils suction surfaces. On the other hand, the current method is found to far more accurately predict the separation point, providing a considerably improved correlation with experimental data.

The most likely reason for this dramatic improvement is the marked transition at the separation point from adverse pressure gradient to constant pressure distribution. As the separation region encroaches on the suction surface, there is a distinct increase in the upstream pressure gradient, which further increases the likelihood of upstream ow separation.

An obvious disparity between the experimental data and that of the current method, which must be addressed, is that for mildly separated ows, it appears that the assumption of constant pressure in the separation region is not entirely realistic, and that there is in fact a slight pressure gradient in this region. The author contends that the pressures in this region are typically very close to zero, and therefore do not contribute significantly to the total forces induced on the airfoil. Rawlinson-Smith [26], who employs the shear layer model, explains that this disparity is due to the ow in this region being dominated by the boundary layer behaviour. He defends the assumption

123

Int J Energy Environ Eng (2016) 7:112 11

Fig. 9 Comparison of pressure distribution data for DU 91-W2-250. a a = 7.7, b a = 9.7, c a = 11.7, d a = 15.2

of constant pressure in the separated wake as being within engineering accuracy.

Conclusions

The ow elds around several wind turbine airfoils, namely NACA 63-430, FFA-W3-301, FFA-W3-241, and DU 91-W2-250, have been calculated using the proposed viscous-coupled 3D panel method. Results have been compared with available experimental data, as well as with numerical data from XFOIL, another viscous-coupled panel method, and EllipSys2D, a RANS-based solver. At low angles of attack, pressure data predicted by the current method show excellent agreement with the experimental data, as do the referenced numerical data. At higher angles of attack, the current method shows reasonable agreement with the experimental data, while the referenced numerical data signicantly overestimate the -Cp distribution along

the suction surface. These discrepancies between the referenced numerical data and experimental data, and the considerably improved correlation achieved by the current method, all serve to highlight the importance of accurately predicting the separation point. Taking into account the considerably improved performance of the current method, combined with the fast computation time and time saved by reducing the dimensionality of the discretisation (compared with RANS solvers), the author believes that the current method could be a powerful tool for the design of wind turbine blades.

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/

Web End =http://crea http://creativecommons.org/licenses/by/4.0/

Web End =tivecommons.org/licenses/by/4.0/ ), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

123

12 Int J Energy Environ Eng (2016) 7:112

References

1. Hess, J.L., Smith, A.M.O.: Calculation of potential ow about arbitrary bodies. Progr. Aerosp. Sci. 8, 1138 (1966)

2. Filkovic, D.: Graduate Work. Faculty of Mechanical Engineering and Naval Architecture, University of Zagreb, Croatia (2008)3. Katz, J., Plotkin, A.: Low-Speed Aerodynamics. McGraw-Hill International Editions, Singapore (1991)

4. Lighthill, M.J.: On displacement thickness. J. Fluid Mech 4(4), 383392 (1958)

5. Maskew, B., Dvorak, F.A.: The prediction of CL max using a separated ow model. J. Am. Helicopter Soc. 23(2), 28 (1977)

6. Fuglsang, P., Antoniou, I., Dahl, K.S., Madsen, H.A.: Tunnel Tests of the FFA-W3-241, FFA-W3-241 and NACA 63-430 Airfoils. Ris-R-1041(EN), Ris National Laboratory, Denmark (1998)

7. Timmer, W.A., van Rooy, R.P.J.O.M.: Wind tunnel results for a 25 % thick wind turbine blade airfoil. In: Proceedings of EUWEC, Lbeck-Travemnde, Germany (1993)

8. Bertagnolio, F., Srensen, N., Johansen, J., Fuglsang, P.: Wind Turbine Airfoil Catalogue. Ris-R-1280(EN), Ris National Laboratory, Denmark (2001)

9. Thwaites, B.: Approximate calculation of the laminar boundary layer. Aeronaut. Q. 1, 245280 (1949)

10. Head. M. R.: Entrainment in the turbulent boundary layer. Technical report, Aeronautical Research Council R&M 3152 (1960)

11. Wauquiez, C.: Shape optimization of low speed airfoils using MATLAB and automatic differentiation. Licentiates thesis, Department of Numerical Analysis and Computing Science, Royal Institute of Technology, Sweden (2000)

12. Robinson, D.E.: Implementation of a separated ow panel method for wall effects on nite swept wings. M. S. thesis, Department of Aeronautics and Astronautics, Massachusetts Institute of Technology, USA (1988)

13. Ramos Garca, N.: Unsteady viscousinviscid interaction technique for wind turbine airfoils. Ph.D. thesis, Department of Mechanical Engineering, Technical University of Denmark, Denmark (2011)

14. Drela, M., Giles, M.B.: Viscousinviscid analysis of transonic and low Reynolds number airfoils. AIAA J. 25(10), 13471355 (1987)

15. Veiga, A. E. L.: The analysis of partially separated ow on sail systems using a sectional method. Ph.D. thesis, School of Engineering Sciences, University of Southampton, England (2006)

16. Riziotis, V.A., Voutsinas, S.G.: Dynamic stall modelling on airfoils based on strong viscousinviscid interaction coupling. Int.J. Numer. Methods Fluids 56, 185208 (2008)17. Ramos Garca, N., Srensen, J.N., Shen, W.Z.: A strong viscous inviscid interaction model for rotating airfoils. Wind Energy 17(12), 19571984 (2013). doi:http://dx.doi.org/10.1002/we.1677

Web End =10.1002/we.1677

18. Ramos Garca, N.: Development of a three-dimensional viscous inviscid coupling method for wind turbine computations. In: Proceedings of ICOWES2013, Denmark, pp. 6981 (2013)

19. Veldman, A.E.P., Carter, J.P.F.: New quasi-simultaneous method to calculate interacting boundary layers. AIAA J. 19, 7985 (1981)

20. Veldman, A.E.P.: A simple interaction law for viscousinviscid interaction. J. Eng. Math. 65, 367383 (2009)

21. Bijleveld, H.A.: Application of a quasi-simultaneous interaction method for the prediction of three-dimensional aerodynamic ow over wind turbine blades. Ph.D. thesis, Rijksuniversiteit Groningen, The Netherlands (2013)

22. Smith, L.: An interactive boundary layer modelling methodology for aerodynamic ows. MEng thesis, Department of Mechanical and Aeronautical Engineering, University of Pretoria, South Africa (2012)

23. Sturdza, P., Suzuki, Y., Martins-Rivas, H., Rodriguez, D.: A quasi-simultaneous interactive boundary-layer model for a Cartesian euler solver. In: 50th AIAA Aerospace Sciences Meeting, Nashville, Tennessee (2012)

24. Drela, M.: XFOIL: an analysis and design system for low Reynolds number airfoils. In: Conference on Low Reynolds Number Airfoil Aerodynamics, University of Notre Dame (1989)

25. Srensen, N.N.: General Purpose Flow Solver Applied to Flow over Hills. Ris-R-827(EN), Ris National Laboratory, Denmark (1995)

26. Rawlinson-Smith, R.: Computational study of stalled wind turbine rotor performance. Ph.D. thesis, School of Mechanical Engineering, Craneld Institute of Technology, UK (1991)

123