Numerical Framework for Aerodynamic Characterization of Wind Turbine Airfoils: Application to Miniature Wind Turbine WiRE-01

: A numerical framework for the aerodynamic characterization of wind turbine airfoils is developed and applied to the miniature wind turbine WiRE-01. The framework is based on a coupling between wall-resolved large eddy simulation (LES) and application of the blade element momentum theory (BEM). It provides not only results for the airfoil aerodynamics but also for the wind turbine, and allows to cover a large range of turbine operating conditions with a minimized computational cost. In order to provide the accuracy and the ﬂexibility needed, the unstructured ﬁnite volume method (FVM) and the wall-adapting local eddy viscosity (WALE) model are used within the OpenFOAM toolbox. With the purpose of representing the turbulence experienced by the blade sections of the turbine, a practical turbulent inﬂow is proposed and the effect of the inﬂow turbulence on the airfoil aerodynamic performance is studied. It is found that the consideration of the inﬂow turbulence has a strong effect on the airfoil aerodynamic performance. Through the application of the framework to WiRE-01 miniature wind turbine, a comprehensive characterization of the airfoil used in this turbine is provided, simplifying future studies. In the same time, the numerical results for the turbine are validated with experimental results and good consistency is found. Overall, the airfoil and turbine designs are found to be well optimized, even if the effective angle of attack of the blades should be reduced close to the hub. study is performed on two representative cases with Re = 10 4 and AoA = 4 ◦ and 8 ◦ . The default mesh ( ≈ 6 × 10 6 cells) is compared to a mesh about twice as coarse ( ≈ 3 × 10 6 cells). The coarser mesh follows the same mesh topology as the default one, and due to the time-stepping method used, the convergence study is not only performed in space but also in time.


Introduction
Nowadays, many of the biggest technological challenges stem certainly from the need for renewable, clean energy and wind energy is clearly a part of the solution. This calls for a better understanding of the complex fluid flow phenomena involved in wind energy, such as the wind turbine wakes interaction with the surrounding turbulent atmospheric boundary layer (ABL), including thermal and topographic effects, and the interaction between multiple wind turbine wakes in wind farms [1,2].
The starting point of most wind turbines aerodynamics is the airfoil aerodynamics. Therefore, when designing or studying wind turbines, an accurate aerodynamic characterization of the airfoils used in the turbine design is desirable. For example, when performing BEM-based simulations, airfoil aerodynamic data are among the required inputs [3]. Recently, numerous studies have pointed out the importance of this topic and investigated it [4][5][6][7][8][9][10][11]. In particular, Sarlark et al. [4] investigated the aerodynamic performance of the NREL S826 airfoil using a wind tunnel experiment and LES.
By comparing the two results, they concluded that LES is able to provide reasonable accuracy. Frère et al. [5] examined the transitional physics of the flow over the same airfoil. They used LES with two different codes, one based on high-order discontinuous Galerkin method (DGM) and the other based on low-order FVM. Results from both codes gave similar results. The numerical results were compared with experimental results and some differences were found. However, they showed that the experimental results can be very sensitive to the details of the experimental configuration, e.g., wind tunnel lateral walls. Due to the consistency of the LES results between the two codes, they concluded that the LES approach can be trusted. Sarlark [8] investigated the aerodynamic performance of the SD7003 airfoil using LES. By a comparison with experimental results, he showed again that LES can produce reasonably accurate results. In addition, the role of the sub-grid scale (SGS) modeling on the results was studied. It was found that the SGS model has a small influence for small Reynolds number but can have a large influence for a large Reynolds number.
Based on the operating conditions of wind turbines, wind turbine airfoils are expected to experience a wide range of Reynolds numbers (Re) and angles of attack (AoA), including large angles of attack. Complex phenomena such as separation bubbles (SBs) are also expected to be critical for their aerodynamic behavior. These elements call for the use of wall-resolved LES in order to provide an accurate simulation of the flow. Numerous LES investigations of airfoil flows have been conducted by the aerodynamic research community. The collaborative study from the LESFOIL project [12] is a pioneering work in this field, providing a general overview of the most important aspects to be considered in LES of airfoil flows. However, this study focuses on a specific case at high Reynolds number and faces the issue of the limited computational resources available at that time. Multiple recent wall-resolved LES studies have been conducted using two main approaches. The first is based on high-order methods (with mainly structured grids), e.g., [5,[13][14][15][16], and the second is based on low-order methods (with mainly unstructured grids), e.g., [4,5,8,17,18]. Both approaches seem to be equally successful. In this study, a low-order method with unstructured grids is used in order to provide the flexibility needed to cover the large span of operating conditions of the turbine blade sections. In particular, the open-source computational fluid dynamics (CFD) toolbox OpenFOAM [19,20] is used. OpenFOAM-based LES has been validated in several studies, e.g., [21][22][23][24][25], and has been already used with success in high Reynolds number airfoil flows; e.g., Roohi et al. [26] investigated the cavitation around the Clark-Y hydrofoil and Jawahar et al. [27] studied the self-noise of the NACA 6512-63 airfoil.
Since the flow around wind turbine airfoils represents the flow experienced by the blade sections of the turbine, two important elements should be taken into account. First, a turbulent inflow should be considered since the blade sections of the turbine operate in the ABL. However, LES studies of airfoil flows are mostly conducted with laminar inflows, e.g., [4,8,[13][14][15][16][17][18]26,27] and very few examples of LES studies of airfoil flows using turbulent inflows are found in the literature. One of the few examples is the work of Gilling et al. [28], in which they used detached eddy simulation (DES) to study the flow around the NACA 0015 airfoil at high Reynolds numbers. In particular, they showed that taking into account the inflow turbulence can increase the accuracy of the results. Here, a practical turbulent inflow definition is proposed in order to represent the turbulence experienced by the blade sections of the turbine, and the effect of the inflow turbulence on the airfoil performance is studied. Secondly, the airfoil simulations should be coupled with the turbine flow. In this context, a procedure based on the BEM is developed. This leads to a coupled LES-BEM framework which provides a characterization of both the airfoil and the turbine aerodynamic coefficients, for prescribed turbine operating conditions. The wind turbine power and thrust coefficients obtained are validated with experimental results [29,30]. To the best of the authors' knowledge, it is the first time that such a framework has been proposed. The important characteristics of this framework are: (1) the lift and drag coefficients of the turbine airfoil are computed from wall-resolved LES and cover a prescribed range of turbine operating conditions; (2) a practical way to consider the turbulence experienced by the blade sections is proposed/used; (3) an LES-BEM coupling procedure is proposed/used; (4) the turbine thrust and power coefficients are also computed; and (5) the computational cost is minimized.
One of the ways to study wind turbine flows is to perform wind tunnel experiments using miniature wind turbines. Compared to their field-scale counterparts, miniature wind turbines flows are often characterized by much lower Reynolds numbers, at which a simple down-scaling of the turbine geometry leads to unrealistic turbine performance. This, in turn, has a negative impact on the wind turbine flows generated in the wind tunnel experiments, since such flows are highly dependent on the performance of the turbine. To address this issue, a new miniature wind turbine has been developed in [29,30] with similar power and thrust coefficients to those of field-scale wind turbines. Another way to study wind turbine flows is to perform numerical simulations. Those need to be validated, mostly against wind tunnel experiments, and can then be used for field-scale problems. Following this approach, wind tunnel experiments and numerical simulations can be used in a very complementary way. The new miniature wind turbine developed in [29,30] has a diameter of D = 0.15 m, a hub height of z H = 0.125 m and a cord length c ≈ 0.01 − 0.02 m. In addition, the miniature wind turbine is based on a design angle of attack of AoA = 4 • and blade sections designed with a thin airfoil with 5% thickness, 5% centered camber and sharp leading and trailing edges, c.f. Figure 1. It has been shown that this kind of airfoil is suitable for low Reynolds numbers, e.g., [31][32][33][34][35]. Since no aerodynamic study of this new airfoil design is available at the moment, the present study aims to provide such an investigation. In the same time, the study allows for a characterization of the overall turbine aerodynamic design and offers insights for its further optimization. In addition, it establishes the aerodynamic database required in numerical studies of wind turbines through advanced modeling methods such as the actuator disk model with rotation (ADM-R) [36], simplifying possible validation studies. The rest of the paper is structured as follows. In Section 2, the numerical method is introduced, including the large eddy simulation model, the turbulent inflow definition, a comprehensive sensitivity study and the LES-BEM coupling procedure. In Section 3, the aerodynamic behaviors of the airfoil and turbine are analysed. In this context, the results for the turbine are compared to the experimental results from [29,30] and a database of the aerodynamic performance of the airfoil is presented. Finally, the conclusions are given in Section 4.

Numerical Method
The governing equations corresponding to the filtered incompressible Navier-Stokes equations can be expressed as: whereũ is the filtered velocity,p =p + 2 3 k SGS ,p the filtered kinematic pressure, k SGS the SGS turbulence kinetic energy, ν e f f = ν + ν SGS , ν the molecular kinematic viscosity and ν SGS the SGS kinematic viscosity. In addition, the equations are non-dimensionalized with the airfoil free-stream velocity U and the airfoil chord length c. The equations are solved based on the cell-centered unstructured FVM within the OpenFOAM toolbox [19]. The filtering is taken as the finite volume based implicit filtering or equivalently the top-hat filtering with the filter width taken as ∆ = V 1/3 , where V is the volume of the cell [37,38]. The SGS contributions are modeled through the wall-adapting local eddy-viscosity (WALE) model [37,39]. The WALE model is chosen for three main reasons. First, it is known to capture the correct near-wall behavior [39][40][41] and is thus particularly well adapted for the present study, where complex near-wall flows are expected. Secondly, it is computationally relatively inexpensive, which is a significant advantage in this study, where the SGS contributions are small since a very fine resolution is used. Finally, airfoil studies with very similar numerical setup and the WALE model have been successfully performed with OpenFOAM, e.g., [23].
The equations are solved based on a collocated grid arrangement and the Rhie-Chow interpolation [42][43][44]. They are coupled through the segregated PISO algorithm with two corrector steps [43][44][45]. Centered linear interpolation is used in the finite volume spatial discretization of all terms except for the convective term where the linear upwind stabilized transport interpolation (LUST) is used [43,44,46]. The LUST interpolation scheme is a blend with 75% centered linear interpolation and 25% linear upwind interpolation and is known to be well optimized in terms of accuracy and stability for the kind of problem presently considered [23,24,46]. The implicit second order backward scheme is used in the temporal discretization with an adaptive time stepping [43,44]. Considering all these elements, the overall method is expected to be 2nd order accurate [43][44][45].
The coordinate system is defined with x direction parallel to the free-stream, z direction parallel to the span and y direction orthogonal to x and z directions. The airfoil in centered in the origin, i.e., in (x/c = y/c = z/c = 0). The domain is taken as a rectangular parallelepiped with −4.5 ≤ x/c ≤ 25.5, −7.5 ≤ y/c ≤ 7.5, −0.075 ≤ z/c ≤ 0.075. The boundary conditions are defined with a velocity inlet and a pressure outlet in x direction, symmetry boundary conditions in y direction, periodic boundary conditions in z direction and a no-slip on the airfoil. The airfoil geometry definition follows the airfoil geometry shown in Figure 1, with the sharp leading and trailing edges slightly rounded with a diameter of d/c ≈ 0.004, corresponding to the turbine manufacturing resolution [29].
The meshing process is performed with OpenFOAM automatic mesh generator snappyHexMesh [20] which generates unstructured meshes with mainly hexahedral cells and prescribed mesh quality. This mesh generator enables flexible local refinement capabilities and a fully automatic meshing process, which is essential for the LES-BEM coupling. The computational mesh is composed mainly of two different regions corresponding to the far-field and the near-wall regions, both of which are structured and aligned with the local flow. In addition, the mesh is highly refined around the airfoil and the wake regions. The near-wall region is composed of only hexahedral cells and extends over 14 cells. In this region, the mesh is stretched in the wall normal direction with a cell-to-cell expansion ratio of R = 1.1. The spatial resolution at the wall is defined such as ∆n + ≈ 0.2, ∆q + = ∆s + ≈ 6 · ∆n + , where ∆n, ∆q, ∆s are the size of the wall cells in the wall normal, streamwise and spanwise directions, ∆ + = u * · ∆/ν, u * is the friction velocity. It should be noted that the resolution is much finer than the common practice for wall-resolved LES [12]. This is especially important in the present low-order LES approach for which the spatial resolution should be very fine in order to minimize the numerical dissipation [24]. Meanwhile, the full 3D unstructured approach allows a relatively small number of mesh cells despite the very fine resolution with a total number of cells around 6 × 10 6 . An example of the computational mesh for AoA = 10 • is shown in Figure 2. Finally, the ratio of the modeled turbulence kinetic energy to the total turbulence kinetic energy is used as mesh resolution criterion [47] and the value is kept under 20% in the region of interest and for all cases considered.
The adaptive time stepping is defined such that the cell Courant number Co < 2. This can be interpreted as a unstructured time stepping and ensures small time steps when needed and larger time steps when possible. Since the spatial resolution is very fine, the selected Courant number threshold leads already to very fine time steps with ∆t ≈ 0.001 · c/U. In particular, the time step is similar to the one used in [8]. Finally, the statistics are computed from t 1 = 15 · c/U until t 2 = 45 · c/U, leading to an averaging time of t 2 − t 1 = 30 · c/U. It is verified for all cases that this definition is adapted based on the mean flow properties and the time evolution of the aerodynamic coefficients.

Turbulent Inflow
In this study, a turbulent inflow is considered in order to reproduce as closely as possible the operating conditions of the wind turbine blade sections. In this context, the flow experienced by the whole turbine and the one experienced by the blade sections should be differentiated. In the following, the "turbine inflow" refers to the incoming turbulent flow experienced by the whole turbine and the "airfoil inflow" refers to the incoming turbulent flow experienced by the blade sections of the turbine. The following refers only to the definition of the airfoil inflow turbulence. By taking the ABL nature of the turbine inflow into account, the airfoil inflow turbulence is different for every operating condition of the turbine, every radial position on the turbine and every time of the rotation period. Here, it is chosen to define a unique and representative airfoil inflow turbulence for obvious practical reasons. The following assumptions are used in the definition of the airfoil inflow turbulence: (1) the effect of turbine dynamics is neglected; (2) the turbine and airfoil operating conditions are fixed at characteristic conditions with the possible operating conditions evolving around them; (3) the turbulence is assumed isotropic since the blade sections are experiencing many different relative directions of motion. The airfoil inflow is then imposed using the divergence free synthetic method (DFSM) [48]. In particular, following the previous assumptions and the experimental conditions [30], the airfoil inflow is defined with its turbulence intensity I and non-dimensional turbulent length scale l computed as following: where I t is the characteristic turbulence intensity of the turbine inflow, U t /U is the characteristic ratio of the turbine inflow wind speed to the airfoil inflow wind speed and δ t /c is the characteristic ratio of the thickness of the ABL to the chord of the blade sections. It should be noted that the DFSM is not expected to generate a perfect turbulent inflow. In fact, it is designed to accelerate the process of obtaining the statistically converged turbulent flow when the turbulence generation is part of the simulation, e.g., in turbulent boundary layer simulations [48]. This is not considered as an issue since the present turbulent inflow is artificial and only considered in order to reproduce as closely as possible the wind turbine blade sections' operating conditions. Overall, the results show that the turbulent inflow experienced by the airfoil in the simulations is reasonable enough for the purpose of this study.

Sensitivity Study
In order to optimize the numerical configuration, a comprehensive sensitivity study is performed. This is especially important since many numerical simulations must be performed and no validation data for the airfoil are available.
First, a general sensitivity study is performed on a test case with Re = 10 4 and AoA = 4 • , which can be considered as the design case [29]. The default xy domain (−4.5 ≤ x/c ≤ 25.5, −7.5 ≤ y/c ≤ 7.5) is compared with a xy domain twice as large (−9 ≤ x/c ≤ 51, −15 ≤ y/c ≤ 15). The domain length in the x direction is very important in terms of computational cost since the turbulent flow imposed at the inlet takes a longer time to pass through the domain for larger domains. In particular, the domain length in the x direction is identical to the one used in [28], where a turbulent inflow is considered. The comparison of the aerodynamic coefficients obtained for the two domains and presented in Table 1 shows that the difference in those coefficients is small and that the default xy domain is large enough. In addition, the default z domain (−0.075 ≤ z/c ≤ 0.075), corresponding to a span s/c = 0.15, is compared with a z domain twice as large (−0.15 ≤ z/c ≤ 0.15), corresponding to a span s/c = 0.3. The domain span is essential to allow for the 3D turbulent structures to develop. In particular, the default span is very near the value of s/c = 0.2 prescribed in [49]. The comparison of the aerodynamic coefficients obtained for the two domains and presented in Table 2 shows that the difference in those coefficients is small and that the default z domain is large enough. Furthermore, the default Courant number threshold (Co < 2) is compared with a Courant number threshold twice as small (Co < 1). Obviously, the Courant number threshold is very important in terms of computational cost. The comparison of the aerodynamic coefficients obtained for the two Courant number thresholds and presented in Table 3 shows that the difference in those coefficients is small and that the default Courant number threshold is small enough.   z Domain  Secondly, a convergence study is performed on two representative cases with Re = 10 4 and AoA = 4 • and 8 • . The default mesh (≈6 × 10 6 cells) is compared to a mesh about twice as coarse (≈3 × 10 6 cells). The coarser mesh follows the same mesh topology as the default one, and due to the time-stepping method used, the convergence study is not only performed in space but also in time.
The comparison of the aerodynamic coefficients obtained with the two resolutions and presented in Tables 4 and 5 shows that the difference in those coefficients is small. It can be concluded that the coarse resolution is already fine enough. However, many different cases are considered, with some cases potentially requiring finer resolution. Therefore, knowing that the computational cost of the simulations with the fine resolution is reasonable, as a precaution, the fine mesh configuration is chosen for all cases. Table 4. Sensitivity study on mesh resolution for Re = 10 4 and AoA = 4 • . Default correspond to a mesh with ≈6 × 10 6 cells, coarse corresponds to a mesh with ≈3 × 10 6 cells. δ is the relative difference with the default mesh resolution case.  Table 5. Sensitivity study on mesh resolution for Re = 10 4 and AoA = 8 • . Default correspond to a mesh with ≈6 × 10 6 cells, coarse corresponds to a mesh with ≈3 × 10 6 cells. δ is the relative difference with the default mesh resolution case. Finally, a sensitivity study on the inflow turbulence is performed on two representative cases with Re = 10 4 and AoA = 4 • and 8 • , respectively. The results of the sensitivity study are shown in Tables 6 and 7. For AoA = 4 • , the results show a small influence of the inflow turbulence, with a difference of around 1% between the cases. For AoA = 8 • , the results show a large influence of the inflow turbulence, with a difference of around 40% between the cases with and without inflow turbulence and a difference of around 4% between the cases with nonzero inflow turbulence. In addition, it seems that the results are much more sensitive to the presence or not of inflow turbulence than to the non-zero intensity of the turbulence. The difference between the sensitivity to the inflow turbulence for AoA = 4 • and AoA = 8 • can be explained by the different flow behaviors: for AoA = 4 • the flow does not show any SB, while for AoA = 8 • the flow does show one SB, c.f. Section 3. Since SBs are expected to be very sensitive to the details of the inflow, it is consistent to find a small influence of the inflow turbulence for AoA = 4 • and a large one for AoA = 8 • . As a consequence, it is found that the inflow turbulence should be considered and that the exact non-zero level of the turbulence is not of major importance. This last point justifies the idea of representing the multiple different turbulent inflows experienced by blade sections by one unique representative turbulent inflow. Table 6. Sensitivity study on the inflow turbulence intensity I for Re = 10 4 and AoA = 4 • . δ is the relative difference with the turbulence intensity I = 0% case.  Table 7. Sensitivity study on the inflow turbulence intensity I for Re = 10 4 and AoA = 8 • . δ is the relative difference with the turbulence intensity I = 0% case.

LES-BEM Coupling
In order to couple the airfoil simulations with the wind turbine flow, an automatic and iterative procedure based on BEM is developed. The procedure computes the airfoil aerodynamic coefficients (C L , C D ) with LES and the airfoil operating conditions (Re, AoA) with BEM, for a prescribed range of turbine operating condition. In this case, a prescribed range of tip-speed ratios (λ) and a fixed turbine inflow wind speed U t are defined following [30]. In the same time, the BEM also computes the turbine performance parameters (C P , C T ). In this context, the BEM implementation follows [3]. The procedure starts with a guessed and coarse resolution of airfoil operating conditions (Re, AoA), and iteratively refines them, up to a predefined resolution. This allows to increase and refine the database of the airfoil aerodynamic data while computing only the interesting/necessary LES cases. The automation of the procedure is mainly made possible by the use of an automatic mesh generation procedure. In this context, the refinement of the wake region is automatically adapted to the angle of attack.

Results
The prediction of the turbine performance is shown in Figure 3. The predicted power coefficient shows a negative quadratic behavior, with the maximum value C P = 0.36 for λ = 3.45. The predicted thrust coefficient shows a consistent increase from C T = 0.71 for λ = 2.85 to C T = 0.93 for λ = 5.64. It can be seen that the numerical results are consistent with the experimental ones. Especially, the numerical and experimental values of the maximum power coefficients and the optimal tip-speed ratios are very similar. Still, the power coefficient is mainly underestimated and the thrust coefficient is mainly overestimated by the simulations. This could be explained by an overestimation of the lift to drag ratio near the hub and an underestimation of the lift to drag ratio near the tip, which would be consistent with a delayed stall of the LES flow. However, it is difficult to discern how much of the difference is due to the errors in the LES, the BEM or even the experiment. Overall, the difference between the numerical results and the experimental ones is in the order of 5-15%. Considering the expected errors of the coupled LES-BEM simulation and the experiment, this result is considered satisfactory.
The prediction of the airfoil operating conditions is shown in Figure 4. The black dots emphasize the LES cases considered. The prescribed final resolution is ∆Re = 0.2 × 10 4 and ∆AoA = 2 • and the final number of LES cases is N = 46. It is found that the operating conditions are in the range 0.4 × 10 4 < Re < 1.6 × 10 4 and 0 • < AoA < 18 • . For small tip-speed ratios (λ < λ opt ), the operating conditions of the airfoil show mainly small Reynolds numbers (Re < 10 4 ) and a large span of angles of attack, with mainly large angles of attack (AoA > 8 • ). For large tip-speed ratios (λ > λ opt ), the operating conditions of the airfoil show mainly relatively large Reynolds numbers (Re > 10 4 ) and a small span of angles of attack, with mainly small angles of attack (AoA < 8 • ). As a consequence, the operating conditions show a lower triangular distribution, meaning that the airfoil does not operate at large Reynolds numbers and large angles of attack at the same time. In addition, it is interesting to note that, because of the sparse distribution of the airfoil operating conditions, the coupling procedure saves up to about 50% of LES cases w.r.t. the computation of all cases within 0.4 × 10 4 < Re < 1.6 × 10 4 and 0 • < AoA < 18 • . The airfoil performance at the operating conditions presented is shown in Tables 8 and 9. It is found that the lift/drag coefficient generally increase/decrease with the Reynolds number and both generally increase with the angle of attack. Overall, the lift and drag coefficients are consistent with the expected behavior for thin cambered airfoils at low Reynolds number [31][32][33][34].
It can be seen in Figures 4 that the overall operating conditions are well centered on Re = 10 4 with a deviation of about ±2 × 10 4 , and that the operating conditions around the optimal tip-speed ratio are centered on Re = 10 4 and AoA = 4-8 • . In addition, it can be seen in Tables 8 and 9 that the aerodynamic coefficients are much more sensitive to the angle of attack than to the Reynolds number. For all these reasons, the effect of the angle of attack is now studied in detail at the fixed Reynolds number Re = 10 4 .   Table 9. Drag coefficient C D for the Reynolds numbers Re and the angles of attack AoA experienced by the blade sections for the tip-speed ratios around the optimal one. The aerodynamic coefficients of the airfoil at Re = 10 4 is shown in Figure 5, using an extended range of angles of attack. The lift coefficient shows a nearly linear increase for −4 • ≤ AoA ≤ 10 • and a nearly constant behavior for 10 • ≤ AoA ≤ 24 • , with the maximum value C L = 1.53 for AoA = 24 • . The airfoil is therefore found to be very robust to stall. The drag coefficient shows a positive quadratic behavior with the minimum value C D = 0.058 for AoA = 2 • . As a consequence, the lift to drag ratio shows a strong increase for −4 • ≤ AoA ≤ 8 • and a strong decrease for 8 • ≤ AoA ≤ 24 • with the maximum value C L /C D = 13.6 for AoA = 8 • . The lift coefficient corresponding to the maximum lift to drag ratio is C L = 1.27. Again, the performance of the airfoil follows the expected performance for thin cambered airfoils at low Reynolds number [31][32][33][34].
By a comparison of Figures 3-5, the aerodynamic performance of the turbine can be related to the one of the blade sections. For example, optimum turbine tip-speed ratios (λ ≈ λ opt ) correspond directly to high-performance angles of attack (4 • ≤ AoA ≤ 8 • ). Alternatively, suboptimal tip-speed ratios (λ < λ opt and λ > λ opt ) correspond directly to low-performance angles of attack (8 • ≤ AoA ≤ 20 • and 0 • ≤ AoA ≤ 4 • , respectively). It is interesting to notice that the high performance configuration is consistent with the design of the turbine, especially the design angle of attack. Moreover, the estimated values of the lift coefficient and lift to drag ratio used in the design process correspond very well to the values found here [29]. However, it is found that the turbine design could be further improved by reducing the effective angles of attack near the hub. With the operating conditions of the blade sections being well centered on Re = 10 4 and the performance of the airfoil being not very sensitive to the Reynolds number, the fixed Reynolds number Re = 10 4 is expected to be sufficient to describe the airfoil flow. In order to test this hypothesis, a BEM simulation using only the LES aerodynamic data for Re = 10 4 was performed and very similar results were found (not shown here). This confirms that the Reynolds number Re = 10 4 is sufficient in order to characterize the performance of the blade sections and the focus will be now on this specific Reynolds number.
The airfoil flow for Re = 10 4 and AoA = 0 • is presented in Figure 6. In addition, velocity profiles are presented in Figure A1a. This flow is representative for small angles of attack. The stagnation point is located near the sharp leading edge on the suction side. This leads to one region of flow acceleration on the pressure side of the sharp leading edge and another on the rounded leading edge of the suction side. On the suction side, the flow is almost fully attached, except for the small trailing edge section where a small separation can be seen. On the pressure side, the flow shows a separation from the leading edge to approximately the middle of the airfoil, which reattaches due to the geometry of the airfoil. The wake evolves from the suction side small separation to a turbulent vortex street. It is interesting to notice the shortness of the mean wake due to the thin geometry of the airfoil and the inflow turbulence. As a general consequence, the flow shows a suction area in the pressure side which reduces the airfoil aerodynamic performance.
The airfoil flows for Re = 10 4 and AoA = 4 • and 8 • are presented in Figures 7 and 8, respectively. In addition, velocity profiles are presented in Figure A1b,c. These flows are representative for moderate angles of attack. With the increase of the angle of attack from AoA = 0 • to AoA = 4 • , the stagnation point moves in the direction of the pressure side. This reduces the flow acceleration on the pressure side of the sharp leading edge and increases the one on the rounded leading edge of the suction side. On the pressure side, the flow evolves to a fully attached flow. On the suction side and in the wake, the flow does not evolve much. As a consequence, the flow is almost fully attached in both sides of the airfoil. With the further increase of the angle of attack to AoA = 8 • , the stagnation point moves even further to the pressure side and is now located clearly on the pressure side. The region of flow acceleration on the pressure side disappears and the one on the suction side grows with two peaks, one near the sharp leading edge and the other on the rounded leading edge of the suction side. On the suction side, the flow shows a SB which starts just after the rounded leading edge and ends on the third quarter of the airfoil. This is especially clear by looking at the turbulence kinetic energy and the vortex structures. It is also clear that the SB contributes to the good aerodynamic performance of the airfoil at this angle of attack. On the pressure side, the flow does not evolve much. Finally, the mean wake grows together with the turbulent vortices. The airfoil flow for Re = 10 4 and AoA = 12 • is presented in Figure 9. In addition, velocity profiles are presented in Figure A1d. This flow is representative for large angles of attack. With the increase of the angle of attack to AoA = 12 • , the stagnation point moves further to the pressure side. The flow acceleration on the suction side of the sharp leading edge increases, while the one on the rounded leading edge decreases. On the suction side, the flow shows again a separation which starts just after the rounded leading edge, but without reattachment. The separation point moves in the direction of the leading edge and the turbulence gets stronger. This is especially clear by looking to the turbulence kinetic energy and the vortex structures. On the pressure side, the flow is similar to the one at AoA = 4 • and 8 • with a fully attached flow. In addition, the mean wake grows together with the turbulent vortices. With the further increase of the angle of attack, the previous behavior is confirmed: the stagnation point moves further to the pressure side, the suction side separation point moves further in the direction of the leading edge and the turbulence gets stronger in the wake region (not shown here).
The findings regarding the flow behavior explain clearly the small lift, small drag and small lift to drag ratio at small angles of attack, the large lift, small drag and large lift to drag ratio at moderate angles of attack and the large lift, large drag and small lift to drag ratio at large angles of attack. In addition, the analysis of the flow behavior shows that the high performance at AoA = 8 • is closely related to the SB behavior, which is expected to be sensitive to the details of the inflow. For example, it was found in Section 2 that the performance of the airfoil at AoA = 8 • is sensitive to the airfoil inflow characteristics. Considering the variable inflow experienced by the blade sections in the turbine, it is very likely that the airfoil performance for such angles of attack behaves irregularly in the wind turbine. In contrast, the performance at the design angle of attack (AoA = 4 • ) is not expected to be very sensitive to the details of the inflow, while being very near to the optimal performance. Therefore, the design angle of attack seems to be a good compromise. Finally, the analysis of the flow behavior shows that the rounded leading edge is responsible for a large part of the lift for 0 • ≤ AoA ≤ 8 • , which is interesting considering that this feature was originally added solely due to manufacturing limitations [29].

Conclusions
A numerical framework for the aerodynamic characterization of wind turbine airfoils is proposed and applied to the miniature wind turbine WiRE-01. Thanks to an LES-BEM coupling procedure, the developed framework provides not only results for the airfoil aerodynamics, but also for the entire turbine, and minimizes the number of LES cases required. In order to take into account the turbulent nature of the flow experienced by the blade sections of the turbine, a practical turbulent inflow is proposed. In particular, a representative turbulent inflow is implemented through the DFSM. Based on this, the effect of the inflow turbulence is investigated. It is found that the consideration of the inflow turbulence can have a strong impact on the airfoil performance, especially when the flow shows separation bubbles (AoA ≈ 8 • ).
The predicted turbine power coefficients show a positive quadratic behavior with the maximum values of C P = 0.36 for λ = 0.45. The predicted turbine thrust coefficients show a consistent increase from C T = 0.71 for λ = 2.85 to C T = 0.93 for λ = 5.64. These results are compared with experimental results [29,30] and good consistency is found, with a difference between the simulation and experimental results in the range of 5 − 15%. Especially, the general trends, the maximum power coefficients and the optimal tip-speed ratio are very consistent. The predicted blade sections operating conditions show values in the range of 4 000 < Re < 16 000 and 0 • < AoA < 18 • , mainly centered on Re = 10 4 . The operating conditions corresponding to the optimal tip-speed ratios are mainly centered on Re = 10 4 and 4 • < AoA < 8 • . It is found that the blade sections and turbine performance are much more sensitive to the operating angle of attack than to the operating Reynolds number, and that a Reynolds number of Re = 10 4 is sufficient to describe the blade sections flows.
The airfoil flow behavior can be summarized as follows. At small angles of attack (−4 • < AoA < 4 • ), the flow shows a fully attached flow on the suction side and a separation with a geometrical reattachment on the pressure side. The stagnation point is located in the suction side and the flow produces small lift, small drag and small lift to drag ratio. At moderate angles of attack (4 • < AoA < 8 • ), the flow shows whether a fully attached flow or a short SB on the pressure side and a fully attached flow on the suction side. In any case, a strong acceleration is seen in the rounded leading edge and the flow produces large lift, small drag and large lift to drag ratio. At large angles of attack (8 • < AoA), the flow shows a separation without reattachment on the suction side, and a fully attached flow on the suction side. The flow produces then large lift, large drag and small lift to drag ratio. Overall, the optimal angles of attack are found to be AoA ≈ 4-8 • with C L /C D ≈ 12-14 and C L ≈ 0.8-1.4.
From the results, it is clear that the blade sections operating angles of attack are directly responsible for the turbine performance. The design angle of attack is found to be the best compromise in terms of performance and expected behaviour in the wind turbine. The turbine design is found to be well optimized, except near the hub, where the effective angle of attack should be reduced. The design of the airfoil is also found to be well optimized, showing very good performance and strong robustness at large angles of attack [31][32][33][34][35]. Especially, the geometry of the leading edge seems to be very interesting in terms of performance.
When combining the current airfoil-oriented LES framework with a turbine-oriented LES framework, e.g., [36], a self-contained LES approach for the simulation of wind turbines is possible. This is remarkable considering the range of turbulent scales involved in the problem and made possible by the application of BEM in both the airfoil-oriented and turbine-oriented LES frameworks. In addition, application of the proposed framework to field-scale wind turbines seems quite straightforward, but it has not been demonstrated in this work. Finally, it should be noted that both the proposed framework and the airfoil design could be used for other applications such as micro air vehicles (MAVs).  Acknowledgments: The authors would like to thank Majid Bastankhah for the work related to WiRE-01 turbine. We also would like to thank the three anonymous reviewers whose comments helped improve this work. Computing resources were provided by EPFL through the use of the facilities of its Scientific IT and Application Support Center.

Conflicts of Interest:
The authors declare no conflict of interest.