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

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 −Cp distribution along the suction surface.


Introduction
Panel methods are numerical techniques which solve potential flow problems around complex geometries by employing singularity distributions, or panels, to represent an object's surface and wake. Although panel methods are generally limited to the solution of potential flow problems (i.e. inviscid, incompressible and irrotational flow), there are several boundary layer models which extend the scope of analysable problems to include those involving viscous flow. 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 Navier-Stokes (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 flow problems. A further problem with RANS solvers is the significant amount of time that is generally required for mesh generation, while in panel methods, discretisation is confined to the object surface, allowing for fast, convenient modification 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 difficulties presented by the rotation of the rotor, fluctuations in wind speed and direction, rotor-tower 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 sacrifice of aerodynamic efficiency 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 flow field deviates from that predicted by potential flow and separation becomes unavoidable, thus requiring a solution of the boundary layer. For RANS solvers, this necessitates a very fine mesh resolution (y ? \ 1) near the wall, which is very computationally expensive. For potential flow solvers, this boundary layer solution is typically based on empirical viscous corrections.
In inviscid flow, there is no boundary layer and the flow is everywhere tangent to the body surface, which may therefore be represented by a streamline. In viscous flow, however, the no-slip condition states that the fluid at the body surface has zero relative velocity, resulting in a reduction in the flow rate around the body. Two methods of modelling this reduction in flow rate, first 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 flows 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 airfoil's 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 configuration 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.24-0.3), these airfoils are commonly used at the mid-span and inboard stations of wind turbine blades. Results have been verified against wind tunnel test data taken from studies performed at the VELUX wind tunnel in Østbirk, Denmark [6], and at TU Delft's 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]. Potential theory is based on the assumption of an irrotational flow of an ideal fluid (i.e. inviscid and incompressible, expressed by the continuity equation, r Á u From irrotationality, where r Â u * ¼ 0, u * 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 Laplace's equation, which is the governing equation for potential flow. Two solutions to Laplace's equation are: 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 field to simulate the flow field 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 This is called the Dirichlet boundary condition. The geometry may now be discretised into N P surface panels comprising source and doublet distributions and N W wake panels containing doublet distributions. The boundary condition is specified on each panel at a collocation point ( Fig. 2), and the perturbation potential at each collocation point is calculated by summing the influences of all body and wake panels: The integrals are evaluated by defining 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: and /ðx; y; zÞ ¼ Àl 4p where r ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi x À n ð Þ 2 þ y À n ð Þ 2 þz 2 q : The N P source strengths in (5) are known from boundary condition (4), and the remaining unknowns, namely, the N P ? N W 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 N P , and it is now possible to solve the N P linear algebraic equations for l k .
The induced velocity field on the object's 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 field (Eq. 9): The boundary layer For a steady, 2-D, incompressible flow, the Navier-Stokes equations, which are the equations of motion for viscous flow, reduce to the following form: Within the boundary layer, these equations may be simplified 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.
Before going any further, it is first 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 defined as the distance by which the inviscid surface streamline must be displaced so as to account for the reduction in flow rate in the boundary layer (Eq. 15). Similarly, the momentum thickness, h, is defined 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 coefficient are given by Eqs. (17) and (18), where U e is the flow speed outside of the boundary layer, as predicted by inviscid analysis, and s w 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 defined above, we are able to obtain von Kármán's momentum integral equation: 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 flow separates from the body surface. The methods adopted herein are Thwaites' method for laminar boundary layers [9] and Head's 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) U e m By plotting l and H against k, Thwaites found that the right hand side of Eq. (22) is very well approximated by: Substituting (23) into (22), and solving for h 2 , we get: which may, in turn, be used to calculate k. Thwaites found that laminar separation occurs at k \ -0.09. Due to the difficulty in modelling laminar separation bubbles, it is assumed that at the laminar separation point the flow 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.

Head's method
Head's method considers the entrainment of the free stream flow into the turbulent boundary layer, and defines entrainment velocity, E, as the rate of increase in volumetric flow rate, Q, within the boundary layer: Substituting in displacement thickness (15) and introducing a turbulent shape factor, h , the entrainment velocity may be rewritten: The dimensionless entrainment velocity, E/U e , is determined using the following empirical formulae: Finally, Head adopts the Ludwieg-Tillman skin-friction law to solve for c f : where and a 2nd order Runge-Kutta method is employed to solve ordinary differential Eqs. (19) and (28) for h and H 1 . Turbulent separation is assumed to occur when H 1 \ 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 flow 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 insignificant viscous stresses, and is therefore taken to be a potential flow region. As the separated wake shape is not known a priori, an initial wake geometry is specified according to an empirical relationship between the airfoil's thickness-to-chord ratio (a 'wake length factor' is introduced, defined 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: where n * i and t * i are, respectively, the normal and tangential vectors for the ith wake panel.
As mentioned, the induced velocity field on the object's 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 Maskew's 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 coefficient at the ith panel is now given by: 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 García [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 flow field 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 field 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 author's knowledge, this method of modelling separated flow has not been presented in the literature. The initial wake geometry is specified 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 first linearly interpolate C p,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 z i n to z i n?1 , as calculated by the nth iteration of the secant method: where DC p,i = C p,sep -C p,i , i.e. the difference between the pressure coefficient 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 x n out þ1 sep ¼ x n out sep . This interaction between the inviscid and viscous solvers is known as direct coupling, and is addressed in the following section.

Viscous-inviscid coupling
As mentioned above, direct coupling of the inviscid and viscous solvers refers to the case where the inviscid velocity field U e,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 socalled Goldstein singularity. This is caused by a breakdown in the inviscid-to-viscous hierarchy, where the inviscid flow 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 flow equations into one big system of nonlinear equations, which is solved simultaneously by employing the Newton-Raphson method. This fully simultaneous coupling, while robust, is computationally expensive due to full matrix inversion. More recently, fully simultaneous coupled viscous-inviscid 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 flow. 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 [20][21][22][23].
Despite the aforementioned convergence problems with direct coupling, the current method has been found to satisfactorily simulate severely separated flows 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 modified or adapted for the analysis of a wide range of fluid 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.

Results
In order to assess the validity of the current method's results for the case of a rectangular wing in uniform flow, computed pressure distribution data at mid-span are compared with experimental and numerical data taken from Bertagnolio's ''Wind Turbine Airfoil Catalogue'' [8]. The airfoils analysed are NACA 63-430, FFA-W3-301, FFA-W3-241, and DU 91-W2-250, which, due to their high thickness-to-chord ratios (0.24-0.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 Delft's 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 between 6.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 flow 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 -C p 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 -C p 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 -C p 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 flow separation, and the associated difficulty in modelling such separation, is evident from an analysis of the -C p plots. Once again, the referenced numerical data are found to overpredict the -C p 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 -C p suction peak by more than 82 % at a = 16.3°, and with XFOIL overpredicting the -C p suction peak by more than 95 %, while EllipSys2D is seen to underpredict the -C p 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 -C p 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 flow 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 -C p 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 the -C p 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 final airfoil analysed is the DU 91-W2-250 airfoil, for which t/c = 0.25. Like the FFA-W3-241 airfoil, the flow 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 -C p distribution fairly accurately, with the -C p peak overpredicted by around 10 %. Bear in mind that the DU 91-W2-250 airfoil produces a pronounced -C p peak, even at low angles of attack, and that beside the slightly overpredicted -C p 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 -C p 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 results predicted by XFOIL, EllipSys2D, and the current method.

Discussion
For mildly separated flows, the -C p plots (Figs. 6,7,8,9) show excellent agreement between all data sets. In the presence of severe flow separation, however, it is evident that the referenced numerical methods, XFOIL and Ellip-sys2D, 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 flow separation.
An obvious disparity between the experimental data and that of the current method, which must be addressed, is that for mildly separated flows, 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 flow in this region being ''dominated by the boundary layer behaviour''. He defends the assumption

Conclusions
The flow fields 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 significantly overestimate the -C p 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://crea 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.