Investigation of the three-dimensional flow past a flatback wind turbine airfoil at high angles of attack

Flatback airfoils are airfoils with a blunt trailing edge. They are currently commonly used in the inboard part of large wind turbine blades, as they offer a number of aerodynamic, structural, and aeroelastic benefits. However, the flow past them at high angles of attack (AoA) has received relatively little attention until now. This is important because they usually operate at high AoA at the inboard part of Wind Turbine blades. The present investigation uses Reynolds averaged Navier–Stokes (RANS) and hybrid RANSþ large eddy simulation predictions to analyze the flow in question. The numerical results are validated against previously published wind tunnel experiments. The analysis reveals that to successfully simulate this flow, the spanwise extent of the computational domain is crucial, more so than the selection of the modeling approach. Additionally, a low-drag regime observed at angles of attack before stall is identified and analyzed in detail. Finally, the complex interaction between the three-dimensional separated flow beyond maximum lift (stall cells) with the vortex shedding from the blunt trailing edge is revealed. VC 2021 Author(s). All article content, except where otherwise noted, is licensed under a Creative Commons Attribution (CC BY) license (http:// creativecommons.org/licenses/by/4.0/). https://doi.org/10.1063/5.0055822

Investigation of the three-dimensional flow past a flatback wind turbine airfoil at high angles of attack INTRODUCTION As wind turbines grow in size in order to reduce the levelized cost of energy (LCOE), their blades become more slender and more flexible, posing a significant challenge for designers: How can we increase aerodynamic performance and improve blade aeroelastic behavior with limited chord length?
One of the answers to this technical issue is the use of flatback (FB) airfoils at the inboard part of the blades. 1,2 FB airfoils are airfoils with finite trailing edge (TE) thickness. They provide a number of advantages compared to thin TE airfoils of the same thickness. They have improved performance in terms of maximum lift and also provide reduced surface roughness sensitivity, 3 which is crucial for wind turbines. The use of FB airfoils in blades also increases blade flapwise stiffness and reduces blade weight significantly. 1 From a structural point of view, FB profiles offer more space for blade spars to go through, as they have larger cross-sectional area. Unsurprisingly, these benefits come at the cost of increased drag, caused by the low pressure at the blunt TE. To reduce the drag, a number of TE flow control devices have been proposed in the past. [4][5][6] Given these characteristics, most of the FB airfoil profile studies to date have focused either on the aerodynamic/aeroacoustic performance of the profiles [6][7][8][9] or on the wake characteristics of the flow, [10][11][12][13][14] almost entirely at a ¼ 0 or low angles of attack (AoA). On the contrary, the separated flow over FB profiles has received limited attention up until now. According to Ref. 5, the separated flow past thick FB airfoils is three-dimensional (3D) and coherent structures known as stall cells form. Stall cells are large-scale 3D structures of separated flow that consist of two counter-rotating vortices [15][16][17][18][19] that form on the suction side of airfoils experiencing TE stall. 20 Additionally, the literature suggests that in high Re number experiments and simulations, 13,[21][22][23][24] reduced drag values are observed just before stall. This phenomenon has also received practically no research interest. In fact, while the wake structure downstream of cylinders has received significant attention, this is not the case for the wake of streamlined bodies 25 and even less so for FB airfoils at high Reynolds numbers. 26 However, secondary instabilities in the wake of bluff bodies can be used to implement efficient flow control strategies [27][28][29][30] and their knowledge is crucial. These are knowledge gaps that need to be covered, since FB profiles are now commonly used in wind turbine design. An improved understanding of their characteristics is necessary to gain the greatest benefits from their use.
The combination of the three-dimensional unsteady bluff body wake from the blunt TE with the three-dimensional unsteady Stall Cell (SC) flow makes the numerical investigation of this flow a challenge for computational fluid dynamics (CFD). More traditional methods like Reynolds averaged Navier-Stokes (RANS) solvers can provide only reasonable force predictions for FB airfoils and they tend to nonphysically damp the flow. 31 If wake characteristics are of interest, then higher fidelity approaches like detached eddy simulations (DES), are more suitable. 10,11,32 In general, mesh resolution and the size of the domain are crucial factors. 8,[33][34][35] Until now most investigations at high Reynolds numbers concern low aspect ratio (AR) studies, 8,[33][34][35][36][37][38][39] although it is generally agreed that high AR simulations are required if the wake characteristics are to be captured. SC simulations also need to be of sufficiently high AR, 17,18,[40][41][42] if the phenomenon needs to be captured.
The present study, a continuation of the work presented in Refs. 10, 11, and 32, has a dual aim. On one hand, to identify a suitable numerical approach for the study of the complex flow in question and on the other, to use the validated high-fidelity simulations to analyze the flow and its characteristics. More specifically, we investigate the significance of the flow modeling approach and the spanwise length of the computational domain. Additionally, we intend to identify the low-drag regime at angles of attack just before stall and to investigate the complex wake at AoA beyond stall.
This report continuous with a description of the "Methodology" section, followed by the "Results" section. In the latter, first the effect of the computational domain width is investigated and then the flow past a FB airfoil at high AoA is examined in detail. The paper closes with the "Discussion and Conclusions" section.

METHODOLOGY Numerical framework
The computational tool used throughout the paper is MaPFlow. MaPFlow is an unsteady (URANS) solver developed at the National Technical University of Athens. 12 It is a cell centered CFD Solver that can use both structured and unstructured grids. It is capable of solving compressible flows, as well as fully incompressible flows using the artificial compressibility method. 43 Additionally, incompressible flows with small compressibility effects are feasible using low-Mach preconditioning. In all cases, the convective fluxes are discretized using the approximate Riemann solver of Roe. 44 For the reconstruction of the flow field, a second-order piecewise linear interpolation scheme is used. 12 The viscous fluxes are discretized using a central second-order scheme. The Venkatakrishnan limiter 45 is utilized when need (in cases of shock waves) nevertheless; in this work, no limiting is employed, since the flow is nearly incompressible, and erroneous activation of flux limiting in the smooth region of the flow can degrade the overall accuracy. Finally, the Green-Gauss approach is used to evaluate gradients of the flow at cell centers.
MaPFlow can handle both steady and unsteady flows. Time integration is achieved in an implicit manner permitting large Courant-Friedrichs-Lewy (CFL) numbers. The unsteady calculations use a second-order time accurate scheme combined with the dual time-stepping technique to facilitate convergence. 46 MaPFlow is parallelized using the MPI library in a multiblock fashion in which each processor solves a partition of the original computational domain. MaPFlow's parallel performance has been investigated in High Performance Computing (HPC) platforms where the scalability of the code up to thousands of processors is verified.
Turbulence closures implemented on MaPFlow include the one equation turbulence model of Spalart-Allmaras (SA) 47 as well as the two-equation turbulence model of Menter [k-x Shear Stress Transport (SST)]. 48 Regarding higher fidelity models improved delayed DES (IDDES) as defined in Ref. 49 is employed. In the present study, limited 2D and 3D unsteady RANS with the SA turbulence model are reported and the main focus is on 3D simulations with the IDDES. RANS remain the industry standard when it comes to CFD simulations and are hence included in this investigation. The computational cost of wall resolved large eddy simulation remains prohibitive at this Reynolds number and aspect ratio, hence the use of IDDES was selected as a compromise.

Computational grids
In the present study, 3D simulations using the unsteady RANS approach with the SA turbulence model and the IDDES variant are performed. For IDDES simulations, the generation of a suitable computational grid is critical. On the one hand, in the near-wall region, the grid must be coarse enough to ensure that it will be handled by the RANS part of the model while in the wake the grid must be refined accordingly. Following the authors' previous studies 10,11 in this work, a 12 Â 10 6 cell grid is generated.
In the near airfoil region, a grid consisting of hexahedral cells is generated, while in the rest of the domain, an unstructured grid composed of tetrahedral cells is employed. The hexahedral structured-like region extends 0.1 chords around the airfoil and consists of 100 points in the normal-to-the-wall direction while the airfoil is discretized using 430 nodes.
In this region, the yþ value is approximately 0.1 which makes the near wall region appropriate for RANS simulations. In order to achieve these yþ values the first cell center is located 1 Â 10 À6 nondimensional units from the wall while in the spanwise and chordwise direction the cell is two orders of magnitude bigger. This type of spacing ensures that the IDDES model in the boundary layer (BL) region is switched off. More specifically, the switch from RANS to DDES brand is accomplished by altering the distance from the wall used in the original model according to where l RANS is the distance from the wall (original SA model) and l LES ¼ C DES WD where C DES is an empirical constant and W a low Reynolds number correction. The local length scale D is defined as where C w is an empirical constant, d w is the distance for the wall, h wn is the mesh step in the wall normal direction, and h max ¼ maxðDx; Dy; DzÞ. Since the mesh is constructed in the BL region to have a large AR l DDES ¼ l RANS is ensured. In the spanwise direction, a constant spacing was of Dz ¼ 0:008 (125 spanwise locations) is used, based on the findings of Ref. 11. In the wake region and up to five chords downstream, an unstructured grid is composed of tetrahedral cells. Beyond that point and toward

ARTICLE
scitation.org/journal/phf the computational boundary, the grid gradually coarsens. Regarding the spanwise size of the wing used in the numerical simulations four different AR, namely, 0:2; 0:5; 1; and 2. All the computational grids are identical 2D-wise, while only the length of the wing in the spanwise direction varies. Symmetry conditions were applied at the side boundaries of the computational domain. It is noted here that a rational choice would be that of periodic lateral boundaries; nevertheless, the numerical results are compared with measurements of a finite AR wing, and therefore, symmetry conditions at the sides were preferred. It is noted that it was not attempted to simulate the wind tunnel sidewall boundary layer as no data were available regarding its characteristics. This is a no trivial task, 50 out of the scope of the present investigation. Overall, the symmetry conditions were selected as a suitable compromise.
Regarding the far-field boundary, it was located 100 chords away from the wing to minimize the influence of the external boundary conditions on the simulations. 51 All the simulations consider the flow fully turbulent. Finally, due to the thick TE, the flow was considered unsteady in all cases.

Validation
The present work builds on previous studies, 10,11 where the numerical framework was validated. A comparison between experimental stereo particle image velocimetry (PIV), hot-wire anemometry, and pressure data at a ¼ 0 can be found in Ref. 10, while the grid dependence of results is discussed in Ref. 11. The numerical predictions for the dominant St number are close to the experimental results, but lower, in agreement with Ref. 8. As discussed in Ref. 10, an increase in the mesh size from 12 (which is the mesh size of the present investigation) to 25 M does not alter the dominant St number.
In the present case, the comparison is extended to force coefficients at higher angles of attack ( Fig. 6), surface flow visualization, and wake characteristics (vortex formation length) at a ¼ 0 ( Table I). The discrepancies observed at angles of attack beyond maximum lift are linked with stall cell formation and their tendency to provide bifurcating flows in either experiments 41 or simulations. 52 The interested reader is directed to Refs. 10 and 11 and the comparisons are not repeated here in the interest of brevity.

Measurements
The CFD simulations are compared the wind tunnel data from Ref. 5. The experiments concerned a 30% thick FB airfoil with a 10.6% thick TE (LI30-FB10). The airfoil is designed by adding thickness to the mean camber line of a low-induction rotor profile. 53 In the experiments, the wing model had a chord of c ¼ 0:5 m and an aspect ratio of AR ¼ 2:0, while the chord Reynolds number was 1.5 Â 10 6 . The lift was integrated from pressure tap measurements around the airfoil chord at midspan, while drag was measured with a wake rake for attached flow conditions and it was assumed equal to pressure drag for cases with the separated flow. Stereo PIV and hot wire anemometry measurements were taken in the locations shown in Fig. 1. The 95% confidence interval for the lift and drag values is 1% and 4%, respectively. All measurements refer to free transition cases, but it is noted that fixing transition had minimal effect on the performance of the specific profile. 53 Extensive details on the wind tunnel test setup can be found in Ref. 5.

AR study: How much is enough?
In this section, the results from simulations with different AR are examined, from AR ¼ 0:2 to AR ¼ 2:0. Two different 3D phenomena are involved in the flow in question. On one hand, there is the 3D bluff body shedding downstream of the blunt TE that begins at prestall AoA and on the other hand, we have the 3D stall cells that develop poststall. It is vital to assess what the minimum suitable computational domain is to study this flow. Both IDDES and URANS data are reported in this section.
Bluff body shedding Figure 2 shows the mean lift and drag coefficient values at a ¼ 0 for the different AR values. Results show that while the lift mean value remains largely unaffected, the standard deviation is reduced with AR increase. The drag mean value and standard deviation are reduced as AR increases, too, and the results approach the experimental value. Going from AR ¼ 1 to AR ¼ 2 does not alter either the mean values or the standard deviation of either lift or drag.
For completeness, URANS results are also given in this plot, where the AR ¼ 0 data correspond to the 2D simulations. The latter significantly overpredict drag value and standard deviation, but lift and drag predictions are closer to IDDES values for AR ¼ 1.0, in agreement with Ref. 10. These observations suggest that the inherent 3D character of the bluff body flow is artificially constrained for the low AR cases (AR < 1).

Physics of Fluids
The spectral content of the drag coefficient time series is given in, Fig. 3 left. The Strouhal number for this case is based on the TE height (St ¼ fhTE U1 ). The spectra from all cases are very similar, all having a dominant frequency of St % 0:21, which corresponds to bluff body shedding. This indicates that even the lowest AR simulations are able to capture the dominant phenomenon (von K arm an-B enard vortex street). A second peak is observed at the second harmonic, which is because drag is affected similarly by both vortices shed from the TE. The lift coefficient spectrum (not shown here for brevity) provides only the main frequency peak. Interestingly, the peak amplitude at the dominant frequency is reduced as AR increases up to AR ¼ 1:0. With increasing AR, the 3D disturbances allow for the redistribution of the energy content of the spectrum among the "near-peak" frequencies. Again, it appears that for AR > 1:0; there is no significant difference. In agreement with the previous observation, that low AR simulations artificially constrain the 3D character of the flow, it appears that the lower the AR the more dominant the 2D von K arm an-B enard shedding appears to be in the flow.
In order to identify the spanwise correlation length in the flow, the Pearson correlation coefficient (r) of the pressure signal at different spanwise locations is presented in Fig. 4. The correlation coefficient is calculated with respect to the pressure at a point located on the TE, at midheight and at midspan [i.e., at ðx=c; y=c; z=bÞ ¼ ð1; 0; 0Þ]. The flow appears strongly correlated in the spanwise direction for the low AR cases (r % 1 for AR ¼ 0:2 and 0:5). On the contrary, for AR ! 1:0; the maximum spanwise correlation length appears to be k % 0:3c, where c is the chord length. In agreement with previous observations, this suggests that the 3D character of the flow requires a wide enough computational domain (AR > 0:5) to fully develop.
These correlation lengths are linked to the behavior of the von K arm an-B enard vortices shed from the airfoil's blunt TE. In practice, these spanwise vortices experience oblique shedding, vortex dislocations and breakdown. Looking at the concentrated vortices (visualized by isosurfaces of the Q criterion 54,55 ) in Fig. 5, it appears that for such phenomena to be modeled in a CFD simulation the AR needs to be greater than at least 0:5, or if expressed in terms of the TE height (h TE Þ, the computational domain width should be greater than $5h TE .  Moving to even higher AoA, we see in Fig. 6, left, that a C l ;max depends on the computational domain AR: the higher the AR the lower the predicted a C l ;max is. For AR ! 1:0 maximum lift appears at a C l ;max ¼ 13 , practically equal to the experimental value of a C l ;max ¼ 12:6 . For lower AR values a C l ;max is overpredicted by 2 .

Physics of Fluids
The reason for this discrepancy is that low AR simulations do not allow for SC formation. As indicated in Refs. 15, 17, and 41, for the case of extruded airfoils, AR values of AR ! 1:0 are required for SC structures to form. Indeed, the time-averaged surface streamlines shown in Fig. 7 reveal the 3D structure of the flow on the airfoil suction surface for the AR ! 1:0 cases, at a ¼ 15 . It is also observed that the three-dimensionality of the flow on the wing surface significantly affects the flow in the wake as well. The separation line vortex and the trailing edge vortex are three-dimensional, growing at the center of the SC, as originally described in Ref. 17.
It is noted that for AR ¼ 1:0 only half a SC is formed, while for AR

Flow analysis
In this section the flow in question is examined in detail in terms of force coefficient mean values and time correlations, as well as frequency domain characteristics and spanwise correlations. Finally, the wake structure is examined in detail. Based on the findings of "Stall

Force coefficients
Time domain. The lift and drag polars for the two different modeling approaches are shown in Fig. 11. It is observed that the agreement at low AoA (0 a 6 ) is good. Then, at a ¼ 7 the lift gradient for the IDDES simulation drops slightly until a ¼ 13 . The lift standard deviation in this AoA range, as shown in Fig. 12, left, is smaller than in the low AoA range. Then at a ¼ 14 (a C l ;max j IDDES Þ the lift coefficient increases slightly, while its variation increases significantly. For the URANS simulation, the lift gradient remains constant up to a ¼ 10 , while it is reduced to a lower constant value in the range 11 a 14 . At a ¼ 14 (a C l ;max j RANS ) the maximum lift value is observed, as in the IDDES simulations. In the reduced lift gradient range (11 a 14 ), the lift standard deviation is reduced, as the error bars indicate and again as is in the IDDES simulations. The poststall lift level is higher in the URANS predictions compared to the IDDES results.
Looking at the drag variation with AoA (Fig. 11, right), it is observed that there is a low-drag pocket for both modeling approaches, albeit at different incidence ranges (IDDES: 7 a 13 , URANS: 11 a 14 ). In that range, the drag coefficient variation is also reduced, as shown in Fig. 12, left for the IDDES

ARTICLE
scitation.org/journal/phf simulations. In the experiments, the low-drag range is 7:5 < a < 13:5 , see Fig. 6, right, practically equal to the range predicted by IDDES. For both IDDES and URANS predictions, this low-drag/drag variation pocket coincides with the reduced lift gradient range, which indicates a more generic change in the flow structure.
The correlation between the IDDES lift and drag coefficient signals is given in Fig. 12, left. It is observed that the correlation increases as the AoA increases until a ¼ 14 , where there is a sudden drop. In the low-drag range (7 a 13 ), the two signals are strongly correlated (r cor > 0:9).
The results (mean force coefficient gradient, standard deviation and correlation) suggest that the prestall regime can be split into two characteristically different regimes, namely, the low AoA and what is named here as the low-drag pocket regime.  Looking closer at the frequency spectra, it is observed that at a ¼ 12 (low-drag pocket), the frequency peak about St % 0:21 is more concentrated than in the case of a ¼ 0 (low AoA range), where more broadband peak appears. The drag spectra also show a secondary peak at the second harmonic. This secondary peak is not as pronounced in the lift spectra. In the poststall region (a ¼ 16 ), the dominant frequency is very low, at St ¼ 0:005 and the bluff body shedding frequency peak is no longer observed.

Physics of Fluids
The St variation with AoA is given in Fig. 14, for all AoA. Only data from the lift spectra are presented here, in the interest of brevity, and it is worth mentioning that the findings would not change if the drag spectra were used. It is found that the bluff body shedding frequency is dominant up to a ¼ 14 , or the a C l ;max . Beyond this angle, very low dominant St numbers are observed.
It is noted at this point that very low frequencies have been observed in 3D separated flow over airfoils and more specifically when stall cells appear. [59][60][61][62] This low frequency has been attributed to a wholesale expansion/contraction of the stall cells 61 and, when normalized with the projection of the airfoil chord (St c;proj ¼ fcsin a ð Þ U1 ), corresponds to Strouhal number of St c;proj $ Oð10 À2 Þ. This nondimensional quantity is also plotted in Fig. 14 and, indeed, values close to the SC-related frequency are observed.
Spanwise correlation. The Pearson correlation coefficient (r) of the drag time series with respect to the midsection (placed at z=b ¼ 0) for the IDDES simulations is shown in Fig. 15, for a ¼ 0 (low AoA), a ¼ 12 (low-drag pocket), and a ¼ 16 (post-stall range). For completeness, the correlation coefficient of the lift, drag, and pressure time series for the IDDES simulations and for all AoA is given in the Appendix, in Fig. 21.
According to Fig. 15, the flow spanwise correlation differs in the three identified flow regimes. At a ¼ 0 , the flow at the ends of the computational domain is not correlated with the flow at the center of the wing. On the contrary, at a ¼ 12 the flow becomes strongly correlated across the wing span. In the post-stall regime (a ¼ 16 ) there is strong negative correlation between the flow at the ends of the computational domain. This is conceivably due to the stall cell that forms on the wing surface, as discussed later.

Wake structure
Bluff body wake flows at sufficiently high Reynolds numbers are inherently three-dimensional, dominated by two types of coherent structures, namely the von K arm an-B enard spanwise vortices and the streamwise vortices, known as rolls and ribs, respectively. Rolls, which are the primary flow instability, experience small and large scale secondary instabilities. 63 Ribs are referred to as the small-scale, while roll dislocations are large scale secondary instabilities. 63 The structures identified in the wake of the wing at various AoA can help explain the observed differences between the flow regimes. Figure 16 shows instantaneous isosurfaces of the Q criterion, at the top, and instantaneous contours of spanwise vorticity, at the bottom, for a ¼ 0 (low AoA), a ¼ 12 (low-drag pocket), and a ¼ 16 (poststall range).
At low AoA, the wake is dominated by the presence of the von K arm an-B enard vortices shed from both edges of the blunt TE. These vortices are three-dimensional with dislocations and oblique shedding, as expected. They also give rise to significant streamwise ribs, which resemble the secondary instabilities described in Ref. 28 although in that case the bluff body Re number was much smaller. At a ¼ 12 , the number of ribs per unit span is considerably reduced and the wake appears significantly more organized than at a ¼ 0 . It is noted that despite the stark difference in the wake structures, the flow over the airfoil is two-dimensional and very similar between the two cases (see surface streamlines in Fig. 16). In the post-stall case (a ¼ 16 Þ; the surface streamlines reveal the trace of half a stall cell, while the wake structures are significantly more complex and three-dimensional. The wake for the low AoA cases (a ¼ 0 , a ¼ 12 ) is examined further in terms of the vortex formation length, L f , and the

Physics of Fluids
ARTICLE scitation.org/journal/phf wavelengths involved in the flow. The formation length is defined as the location of the maximum wake velocity fluctuation at the wake centerline 64 or, alternatively, as the position of the downstream edge of the recirculation zone along the wake centerline, determined based on the sectional streamlines of the mean flow field. 25 Here, for consistency with Refs. 25, 25, and 65, the latter definition is used and the values are given in Table I. The prediction at 0 , L f % 0:75h TE , is in excellent agreement with the wind tunnel measurements 5 and in very good agreement with previous bluff body wake studies, 25,28,65 where L f % 0:82. It is noted that the latter studies concerned flows with lower bluff body Reynolds numbers, Re d ð Þ 50 Â 10 3 , where d is the bluff body thickness. In the present case, the Reynolds number based on the TE thickness is Re TE ¼ qUhTE l ¼ 159 Â 10 3 . Formation length is strongly correlated with base pressure and base drag in bluff body flows 5,66 and, indeed, the formation length increase at a ¼ 12 corresponds to the observed profile drag decrease.
The spatial structure of the primary instability, the von K arm an-B enard vortex street, is quantified by its wavelength, k x , defined as the streamwise distance between two consecutive vortices shed from the same edge of the blunt TE. Again, for consistency with previous studies, k x can be determined according to the equation k x ¼ Uc fs ; where f s is the shedding frequency and U c is the convective velocity. The latter is defined as the mean velocity on the wake centerline at a distance 4h TE from the wing TE. In the present case, this is at x ¼ 1:424c. The values are given in Table I and, again, the 0 value, k x % 3:6h TE ; is in very good agreement with the bluff body studies at lower Reynolds numbers, 25,28,65 where k x % 3:7h TE . In the low-drag regime (a ¼ 12 ), the primary instability wavelength is reduced to k x % 3:3h TE : The small-scale secondary instability (streamwise vortex ribs) is quantified by the spanwise wavelength, k z , which is a function of flow geometry, bluff body length to thickness ratio, Reynolds number, and inflow turbulence characteristics. 67-69 The interaction of ribs and rolls leads to velocity undulations of streamwise velocity on the XZ plane, which can be used to estimate the spanwise wavelength of the secondary instability. 68,70 Figures 17 and 18 show instantaneous Q criterion isosurfaces and the velocity undulations on a XZ plane at y=c ¼ 0:04 for a ¼ 0 and a ¼ 12 , respectively. The spanwise wavelength is found to be k z % 1:0h TE at a ¼ 0 and k z % 1:4h TE at a ¼ 12 . The vortex structure is similar to the model proposed in Refs. 25 and 28, especially at a ¼ 12 . Instantaneous velocity data from Ref. 25 also suggest similar spanwise wavelength values (k z % 1:25d, from Figs. 8 and 14,in Ref. 25). It is noted that in Refs. 25 and 28, a POD based method was used to estimate the spanwise wavelength and the value k z % 2:4d was found. Later, however, it was shown that the implemented analysis was inadequate and it was suggested that a more complete proper orthogonal decomposition analysis would lead to lower spanwise wavelength values, similar to the ones reported here. 71

ARTICLE
scitation.org/journal/phf In an experimental study of an airfoil with relatively thin finite trailing edge (h TE ¼ 1:8% chord), 26 it was found that the spanwise distribution of ribs does not change with incidence. However, that investigation was performed at significantly lower Reynolds numbers (Re ¼ 0:7 Â 10 6 ; Re TE ¼ 1:26 Â 10 4 ) and the examined incidence range did not exceed a ¼ 6 . As discussed below, the increase in AoA at higher Re is an important factor, as the flow remains attached on the suction side of the airfoil, under conditions of increased adverse pressure gradient. Figure 19 shows Q isosurfaces and spanwise vorticity contours for the time-averaged IDDES case at a ¼ 0 (low AoA), a ¼ 12 (low-drag pocket). The Q isosurface for the top vortex is much larger in the first case and, at the same time, the spanwise vorticity contours are much shorter, indicating a thinner boundary layer (BL) compared to the second case. At the same time, the lower vortex appears larger and stronger at a ¼ 12 , while the top vortex is significantly weaker and smaller. This happens because the BL on the suction side of the airfoil must overcome a significant adverse pressure gradient as the AoA increases. This decelerates the flow and the BL that feeds the vortex shed from the top edge of the TE (see also Refs. 10,72,and 73).
These observations suggest that when we consider the unsteady vortex shedding from the wing TE, the vortices shed from the two edges are of comparable strength at a ¼ 0 , but this is not the case at a ¼ 12 , where the vortices shed from the lower edge will be much stronger. It is conceivable that this difference leads to different vortex interaction in the wake and, as a result, to the different wake structure in the low-drag pocket regime (a ¼ 12 ) compared to the low AoA regime (a ¼ 0 ).
The discussion of the complex vortex shedding behavior at high AoA is limited to qualitative terms, given the discrepancies between simulations and experiments highlighted in Fig. 6. Figure 20 shows where spanwise vorticity contours are plotted on three planes normal to the wing span. In this case, which is identical to the case shown in Fig. 10, top right, half a SC appears on the wing. The end of the SC is at the left-hand side of the wing while the center of the SC is at the righthand side, where the separation line is at its most upstream location.
Looking at the vorticity contours in Fig. 20, we observe that the separation shear layer (blue in Fig. 20) is shed from different chordwise locations on the wing suction side, due to the presence of the SC and the three-dimensional separation line. The separation shear layer is hence markedly three-dimensional and so is its roll up into discreet vortices. The separation line is also not fixed and the most upstream separation location varies by 0.2c for the images shown in Fig. 20.
The lower edge shear layer (red in Fig. 20) on the other hand, has a fixed shedding point and is always shed from the lower edge of the blunt TE. As such it always leaves the wing in relatively twodimensional form. The roll-up of the lower shear layer is affected by the three-dimensionality of the top shear layer.
At the left end of the wing, i.e., at the side of the stall cell, the shedding strongly resembles a von K arm an-B enard vortex street. Both shear layers roll up into discrete spanwise vortices in a highly organized manner. At the other end of the wing, i.e., at the center of

DISCUSSION AND CONCLUSIONS
In the present study, the flow past a flatback wind turbine airfoil has been investigated for a range of AoA. Both IDDES and URANS simulations were performed and validated against experimental data. The flow phenomena involved in this case, such as bluff body vortex

ARTICLE
scitation.org/journal/phf shedding and stall cell formation, are inherently unsteady and threedimensional and a challenge for the CFD community. Results show that the selection of CFD modeling approach (URANS or IDDES) is not as important as the spanwise extent of the computational domain. In fact, when applied on the same mesh of AR ¼ 1:0, there is good agreement between IDDES and URANS predictions in (a) force coefficient values at low AoA (b) a C l ;max prediction and (c) the fact that a low-drag pocket regime is predicted in both cases.
On the other hand, the spanwise extent of the computational domain has been found to be crucial. More specifically, it is found that for both 3D phenomena involved in this flow (bluff body vortex shedding and stall cell formation) to develop unrestrained, the AR needs to be greater than 0.5. At lower angles of attack, AR values lower or equal to 0:5 artificially constrain the 3D character of the bluff body wake. Limiting the computational domain width to less than less than 5h TE and enforcing symmetry conditions at the lateral boundaries significantly alters the vorticity evolution in the wake. Still, all simulations predict a dominant frequency of St % 0:21, which suggests that even the lowest AR simulations are able to capture the dynamics of the primary instability (bluff body spanwise vortex shedding).
With regard to higher AoA flow, it is found that when the spanwise extent of the computational domain is limited (AR 0:5), the formation of stall cells is prevented and, as a result, a C l ;max is overpredicted. Stall cells appear slightly earlier for AR ¼ 2:0 than for AR ¼ 1:0, but the a C l ;max remains the same for the two cases, practically equal to that of the experiments. The number of stall cells per unit span is the same for AR ¼ 1:0 and AR ¼ 2:0 simulations for a ! 15 .
The flow past the airfoil was further investigated based on the AR ¼ 1:0 IDDES simulations. To the best of the authors' knowledge, a low-drag pocket regime at AoA just before a C l ;max has been identified and analyzed for the first time. This pocket is observed in the experiments, IDDES and URANS simulations. Previous studies did report lower drag values in the pre-stall region, but the phenomenon was not investigated further. The existence of a separate, low-drag regime is confirmed by lift and drag coefficient mean values and signal behavior along with wake structures, frequency content, and spanwise correlations. Furthermore, the vortex formation length and the wavelength of the secondary instability increase, while the wavelength of the primary instability decreases. Conceivably, this structure change is linked to the reduction and increase in strength of the top and lower vortex, respectively, in the bluff body wake. The vortices are no longer comparable in strength and their interaction in the wake is no longer as strong and three-dimensional as it was at low AoA.
At AoA beyond a C l ;max , the wake is highly complex and threedimensional. Both shear layers, the separation and lower edge, roll-up in discreet vortices. The roll-up is more structured at the side of the SC. On the other hand, at the center of the SC, the roll-up is more complex and less canonical. The presence of the stall cells leads to bent and highly 3D vortical structures.
FIG. 20. Instantaneous contours of spanwise vorticity on three planes parallel to the flow at z/b ¼ À0.5, z/b ¼ 0 and z/b ¼ 0.5, for AR ¼ 1.0 and Re ¼ 1.5 Â 10 6 , a ¼ 15 . Contour levels: red: À10 to À9 and blue: 9 to 10. The last two figures are the same, but on the lower right one, Q criterion instantaneous isosurfaces of Q ¼ 100 colored by spanwise vorticity are also shown to highlight the complexity of the flow. Physics of Fluids ARTICLE scitation.org/journal/phf