Wall-resolved LES modeling of aWind turbine airfoil at different angles of attack

Journal of Marine Science and Engineering Article Wall-Resolved LES Modeling of a Wind Turbine Airfoil at Different Angles of Attack Irene Solís-Gallego, Katia María Argüelles Díaz, Jesús Manuel Fernández Oro and Sandra Velarde-Suárez * Fluid Mechanics Area, Department of Energy, University of Oviedo, 33204 Gijón, Spain; [email protected] (I.S.-G.); [email protected] (K.M.A.D.); [email protected] (J.M.F.O.) *Correspondence: [email protected]; Tel.: +34-985-18-2101 Received: 18 February 2020; Accepted: 15 March 2020; Published: 19 March 2020   Abstract: Noise has arisen as one of the main restrictions for the deployment of wind turbines in urban environments or in sensitive ecosystems like oceans for offshore and coastal applications. An LES model, adequately planned and resolved, is useful to describe the noise generation mechanisms in wind turbine airfoils. In this work, a wall-resolved LES model of the turbulent flow around a typical wind turbine airfoil is presented and described in detail. The numerical results obtained have been validated with hot wire measurements in a wind tunnel. The description of the boundary layer over the airfoil provides an insight into the main noise generation mechanism, which is known to be the scattering of the vortical disturbances in the boundary layer into acoustic waves at the airfoil trailing edge. In the present case, 2D wave instabilities are observed in both suction and pressure sides, but these perturbations are diffused into a turbulent boundary layer prior to the airfoil trailing edge, so tonal noise components are not expected in the far-field noise propagation. The results obtained can be used as input data for the prediction of noise propagation to the far-field using a hybrid aeroacoustic model. Keywords: wall-resolved LES model; wind turbine airfoil; wind tunnel measurements 1. Introduction The growing awareness for the use of renewable energy has led to the development of ambitious projects to meet the increasing energy demands in a less aggressive way with the environment. To that end, wind energy technology has been growing considerably over the last few years. In this area, noise is one of the main restrictions for the deployment of wind turbines in urban environments or in sensitive ecosystems like oceans for offshore and coastal applications. Therefore, aeroacoustics in wind turbine profiles is a technological field in development that is increasingly demanding advances that allow us to reduce the environmental impact generated. Thus, noise, considered until recently as a by-product in wind turbines, tends to become an essential and indispensable objective in the design stages of the new prototypes. Most of the noise generated by wind turbines is aerodynamic because the mechanical noise has been considerably reduced. Therefore, only through an exhaustive knowledge of the fluid dynamics around wind turbine airfoils, the basic mechanisms of noise generation can be determined and relevant actions with the aim of reducing the consequent noise emissions could be performed. Frequently, small wind turbines operate at low-to-moderate Re numbers. Arcondoulis et al. [ 1 ] exposed a classification of the noise generated by airfoils at this range of Re numbers. Among all types of noise, trailing edge (TE) noise is considered one of the major noise generation mechanisms for rotor blades of wind turbines [ 2 , 3 ]. This trailing edge noise limits the use of wind turbines in urban areas and offshore applications due to acoustic impact. TE noise is due to the scattering of the vortical J. Mar. Sci. Eng. 2020,8, 212; doi:10.3390/jmse8030212 www.mdpi.com/journal/jmse J. Mar. Sci. Eng. 2020,8, 212 2 of 18 disturbances in the boundary layer into acoustic waves at the airfoil TE. It is an unavoidable noise source, being the most significant component of the broadband noise in a frequency range from 750 Hz up to 2500 Hz [ 4 ]. Therefore, the understanding of the flow around the airfoil is strongly necessary for aerodynamic and acoustic design purposes. For the comprehension of the flow phenomena involved in the generation of this type of noise, numerical simulation tools are desirable, because they allow a reduction in the number of laboratory tests and optimize the design of new prototypes. Reynolds Averaged Navier–Stokes equation (U-RANS) models, which involve time-averaging to the Navier–Stokes equations to model the turbulent part of the flow, offer the most economic approach for computing complex turbulent industrial flows. They are suitable for many engineering applications and typically provide the required level of accuracy, but they are not suitable to describe 3-D unsteady turbulent flows with the level of precision required in aeroacoustics applications. An alternative to these models are Scale-Resolving Simulation (SRS) models, which resolve at least a portion of the turbulence for at least a portion of the fluid domain (typically, larger and more problematic scales), leaving the turbulence model to account for just the effects of more universal and smaller isotropic scales [ 5 ]. This family includes Large Eddy Simulations (LES), which solve the largest flow scales. A LES model, adequately planned and resolved, is useful to describe the airfoil noise generation mechanisms. In addition, the results obtained can be used as a starting point for the application of hybrid aeroacoustic models, in which the far-field acoustic pressure is predicted from the LES source terms using methods based on Lighthill’s acoustic analogy [ 6 ]. Sol í s-Gallego et al. [ 7 ] have applied Curle’s surface approach [ 8 ] and Ffowcs–Williams and Hall’s volumetric analogy (FW-Hall) [ 9 ] to predict the far-field trailing edge noise radiated by a high-lift wind turbine airfoil. LES modeling techniques require extraordinarily fine meshes and time steps small enough to capture the fluctuations of variables in the scales to be resolved. A correct selection of parameters is essential for obtaining results that faithfully reproduce the physical flow phenomena involved. In particular, the mesh size should be established in such a way that a significant percentage of the turbulent kinetic energy can be resolved so that some experimental information on the turbulence integral length scale should be provided for this purpose [10]. According to the statements exposed above, the aim of this work is the development, application and validation of an LES model of the turbulent flow around a typical wind turbine airfoil. The numerical procedures are presented and described in detail, in order to serve as a guide to another LES modeling works of similar features. The numerical results obtained have been validated with hot wire measurements in a wind tunnel. Based on the analysis of the results obtained, the features of the turbulent flow and the boundary layer developed on the profile are explained. More specifically, the determination of the boundary layer characteristics and its interaction over the airfoil trailing edge allow the main noise generation mechanisms to be identified. These results can also be used as input data for the prediction of noise propagation to the far-field using a hybrid aeroacoustic model. Firstly, the paper describes the numerical methodology, LES simulations and experimental validation. Then, the results obtained are exposed, analyzed and discussed. Finally, a conclusion section outlines the main findings of the work and proposes future developments. 2. Numerical Methodology and LES Computations In this section, the numerical methodology and the LES modeling used for the simulations are described in detail. Firstly, a brief description of the global characteristics presented by the airfoil is given. 2.1. FX 63-137 Airfoil This airfoil, introduced by F.X. Wortmann [ 11 ] in 1972 for the Liver Puffin human-powered aircraft [ 12 ], has been typically used for many low-Re-number applications. It presents a remarkable J. Mar. Sci. Eng. 2020,8, 212 3 of 18 high-lift, soft-stall characteristics and an overall good performance. In the case of small wind turbine facilities, it has been used by several companies, like Aeromag or Southwest Windpower, for the development of the blades of different wind turbines, like the Lakota Unit or the H-40 and H-80 Models [13]. The FX 63-137 airfoil presents a maximum thickness of 13.7% of the chord, located at 30.9% of the chord length. The maximum camber, equivalent to 6% of the chord, is placed at 53.3% of the chord length. For the present study, a reference chord of 0.305 m has been considered for the numerical model (see Figure 1, left). This value has been selected to match the geometrical dimensions of a previous physical prototype made in aluminum (with a 1.1 m span, see figure right). J. Mar. Sci. Eng. 2020, 8, x FOR PEER REVIEW 3 of 19 development of the blades of different wind turbines, like the Lakota Unit or the H-40 and H-80 Models [13]. The FX 63-137 airfoil presents a maximum thickness of 13.7% of the chord, located at 30.9% of the chord length. The maximum camber, equivalent to 6% of the chord, is placed at 53.3% of the chord length. For the present study, a reference chord of 0.305 m has been considered for the numerical model (see Figure 1, left). This value has been selected to match the geometrical dimensions of a previous physical prototype made in aluminum (with a 1.1 m span, see figure right). Figure 1. Airfoil FX 63-137. (a) 2D section. (b) 3D numerical blade. (c) Experimental model. Intensive experimental measurements of the aerodynamic performance of this airfoil have been performed by Selig and McGranaham [13,14] at NREL (US Energy Department). It has been determined that the FX 63-137 produces a maximum lift coefficient (CL,max) of approximately 1.7 for a wide range of Re numbers (between 100,000 and 500,000). It is also generally accepted that for Re = 100,000, the airfoil suffers the consequences of a large laminar separation bubble, especially at the lower angles of attack, with a severe fall of the lift coefficient and a quite drag force high. The situation improves for angles of attack higher than 4°, as the bubble begins to attach to the airfoil, so the lift increases, and the drag is correspondingly reduced. Another significant characteristics of the airfoil are that (1) the FX 63-137 is susceptible for suffering a reduction of the maximum lift performance if simulated roughness is added (it is estimated to be in a drop of 0.2 in the CL,max) and that (2) the airfoil exhibited a soft stall with little unsteadiness. 2.2. Numerical Scheme The commercial CFD software FLUENT® was used to solve the Navier–Stokes set of equations in an incompressible fashion, introducing an unsteady 3D viscous scheme for the finite volume method, with second-order accuracy for the temporal discretization. A bounded second-order upwind formulation has been used for the convection terms, providing a spatial accuracy with a reduced numerical diffusion, especially for complex three-dimensional flows. The diffusion terms are central-differenced and second-order accurate. On the other hand, a pressure-based solver with the SIMPLE algorithm and a Green-Gauss cell-based discretization scheme for the gradient computation has demonstrated an accurate compromise between stability and CPU time. In addition, for the turbulence closure, an LES scheme was used to solve directly the large scales of the flow, modeling the effect of eddies smaller than the grid cell size [15]. For the subgrid-scale model, the Smagorinsky–Lilly model [16] was employed after filtering the incompressible Navier– Stokes equations (in the following, a hat denotes a subgrid average): The continuity equation is given by:    =0. (1) Figure 1. Airfoil FX 63-137. (a) 2D section. (b) 3D numerical blade. (c) Experimental model. Intensive experimental measurements of the aerodynamic performance of this airfoil have been performed by Selig and McGranaham [ 13 , 14 ] at NREL (US Energy Department). It has been determined that the FX 63-137 produces a maximum lift coefficient (C L,max ) of approximately 1.7 for a wide range of Re numbers (between 100,000 and 500,000). It is also generally accepted that for Re =100,000, the airfoil suffers the consequences of a large laminar separation bubble, especially at the lower angles of attack, with a severe fall of the lift coefficient and a quite drag force high. The situation improves for angles of attack higher than 4 ◦ , as the bubble begins to attach to the airfoil, so the lift increases, and the drag is correspondingly reduced. Another significant characteristics of the airfoil are that (1) the FX 63-137 is susceptible for suffering a reduction of the maximum lift performance if simulated roughness is added (it is estimated to be in a drop of 0.2 in the CL,max) and that (2) the airfoil exhibited a soft stall with little unsteadiness. 2.2. Numerical Scheme The commercial CFD software FLUENT ® was used to solve the Navier–Stokes set of equations in an incompressible fashion, introducing an unsteady 3D viscous scheme for the finite volume method, with second-order accuracy for the temporal discretization. A bounded second-order upwind formulation has been used for the convection terms, providing a spatial accuracy with a reduced numerical diffusion, especially for complex three-dimensional flows. The diffusion terms are central-differenced and second-order accurate. On the other hand, a pressure-based solver with the SIMPLE algorithm and a Green-Gauss cell-based discretization scheme for the gradient computation has demonstrated an accurate compromise between stability and CPU time. In addition, for the turbulence closure, an LES scheme was used to solve directly the large scales of the flow, modeling the effect of eddies smaller than the grid cell size [ 15 ]. For the subgrid-scale model, the Smagorinsky–Lilly model [ 16 ] was employed after filtering the incompressible Navier–Stokes equations (in the following, a hat denotes a subgrid average): J. Mar. Sci. Eng. 2020,8, 212 4 of 18 The continuity equation is given by: ∂ˆ ui ∂xi =0. (1) The momentum equation is given by: ρ∂ˆ ui ∂t+ρ∂ˆ uiˆ uj ∂xj =−∂ˆ p ∂xi +µ∇2ˆ ui+∂τij ∂xj , (2) where τij =−ρˆ uiˆ uj−ˆ uiˆ uj is the subgrid-scale stress tensor, modeled with the Smagorinsky–Lilly closure model [16]: τij =2νTˆ Sij (3) νT=L2 Sˆ S=CSˆ ∆2ˆ S. (4) In Equation (4), L S is the mixing length for subgrid scales, νT is the subgrid eddy viscosity, C S is the Smagorinsky constant (0.18 in our case), ˆ ∆ is the local grid size, Sij is the resolved scale strain rate tensor and ˆ S=q2ˆ Sij ˆ Sij. The Smagorinsky model is the simplest and most robust option for the subgrid-scale modeling in LES computations [ 17 ]. Despite the well-known dissipative issues of this model, especially for shock flow situations, a wide number of researchers are still relying on this SGS model for their LES computations, even in complex bladed geometries for turbomachinery (rotor/stator stages), because of its simplicity and versatility [ 18 – 20 ]. Moreover, it is recognized to still provide a good prediction of the important flow patterns and also an accurate reproduction of secondary flow features in complex flows [21]. 2.3. Numerical Domain and Boundary Conditions (BCs) Domain sizes found in the literature in airfoil chord units cfor airfoil simulation dictate typical values around 10 chord lengths for inlet distance and 20 chord lengths for outlet boundaries (see literature survey in [ 22 ]). However, it is not unusual to find reduced domains in the order of 5c upstream and just 10cdownstream for LES schemes [ 23 , 24 ]. For the present investigation, preliminary studies [ 25 ] were performed on a 2D-RANS basis concluding the marginal effect of the boundary conditions on the flow instabilities of the wake flow. As a consequence, a reduced domain comprising 4.92cupstream from the leading edge and 8.84cdownstream from the TE was finally adopted as an optimal selection for CPU costs. Because LES simulations preclude the use of full-3D domains to resolve the three-dimensional anisotropy of the largest flow scales, it is necessary to model the spanwise direction of the airfoil with accurate precision. Since the largest scales in a boundary layer are in the order of δ , and these scales are probably also apparent in the spanwise direction [ 26 , 27 ], the ratio δ /L z should at least be less than one [ 28 ], L z being the spanwise extent of the domain. Taking into account that the boundary layer thickness obtained from the experimental measurements is about 2% of the chord [ 7 ], the spanwise dimension of the domain was set as L z /c=0.164, eight times larger than the maximum length scale expected to be found in the flow. Figure 2shows a sketch with the final dimensions of the numerical domain. Regarding the boundary conditions, a pressure outlet condition was applied in the far boundary, while the velocity-inlet condition was specified for the rest of the far-field boundaries. A moderate-to low Re number of 350,000, based on the chord length, was used at the inlet with a characteristic 0.7% turbulence intensity and an integral length scale of 0.075 m (obtained in [ 29 ] by hot wire experimental measurements). The no-slip condition was used for the airfoil walls and the symmetry condition was applied in the spanwise boundaries (top and bottom). J. Mar. Sci. Eng. 2020,8, 212 5 of 18 J. Mar. Sci. Eng. 2020, 8, x FOR PEER REVIEW 5 of 19 Figure 2. Numerical domain. (a) Upstream and downstream extension. (b) Boundary conditions. Regarding the boundary conditions, a pressure outlet condition was applied in the far boundary, while the velocity-inlet condition was specified for the rest of the far-field boundaries. A moderateto low Re number of 350,000, based on the chord length, was used at the inlet with a characteristic 0.7% turbulence intensity and an integral length scale of 0.075 m (obtained in [29] by hot wire experimental measurements). The no-slip condition was used for the airfoil walls and the symmetry condition was applied in the spanwise boundaries (top and bottom). 2.4. Computational Mesh A block-structured C-mesh of 19M elements, refined at the boundary layer and TE regions, was used for the numerical simulations (Figure 3). Concerning the in-plane mesh requirements [27], recommendations are Δ x+ = 50 to 150 and Δ y+ = 1 for wall-resolved LES resolution. To ensure a good mesh resolution, stream and normal directions near the wall contours were discretized imposing Δ x+ = 45 and Δ y+ = 0.8. These requirements result in typical cell sizes of 0.7 mm and 0.0132 mm in the xdirection and in the y-direction, respectively. These values are also aligned with the value proposed by Pope [17], Δ x+ ~ δ /10 = 0.61 mm. Regarding the number of nodes in the spanwise direction, the guidelines proposed in references [5,17,30] were followed. Once again, in order to obtain a wall-resolved LES resolution, a dimensionless wall distance Δ z+ = 10 to 40 is recommended. Taken 30 as a reasonable value, cell size was fixed to 0.4 mm, which resulted in 15 cells inside the boundary layer, quite close to the 20-cells recommendation inside the boundary layer according to Sagaut [31]. Figure 2. Numerical domain. (a) Upstream and downstream extension. (b) Boundary conditions. 2.4. Computational Mesh A block-structured C-mesh of 19M elements, refined at the boundary layer and TE regions, was used for the numerical simulations (Figure 3). Concerning the in-plane mesh requirements [ 27 ], recommendations are ∆ x + =50 to 150 and ∆ y+ = 1 for wall-resolved LES resolution. To ensure a good mesh resolution, stream and normal directions near the wall contours were discretized imposing ∆ x+ =45 and ∆ y+ = 0.8. These requirements result in typical cell sizes of 0.7 mm and 0.0132 mm in the x-direction and in the y-direction, respectively. These values are also aligned with the value proposed by Pope [17], ∆x+~δ/10 =0.61 mm. J. Mar. Sci. Eng. 2020, 8, x FOR PEER REVIEW 6 of 19 Figure 3. (a) Computational mesh. (b) Minimum cell sizes in the different coordinates. 2.5. LES Modeling and Resolved Scales Additionally, the selection of an appropriate temporal resolution to resolve the turn-out time of the resolved eddies in the LES scheme is also critical. For wall-resolved Large Eddy Simulations, the time step may be calculated as the ratio between the smallest cell size and the fluctuating velocity u’ in the airfoil walls. Assuming Δ y+ ∼ 1, the required time step is ∆𝑡 ≅∆ 󰆓=.· . =2.6·10𝑠, where u’ has been considered as low as 3% of the bulk velocity, according to the experimental measurements for low angles-of-attack (2.5°) in the vicinity of the airfoil wake (𝑢′𝑈 ⁄=0.001). This selection is also compatible with the satisfaction of the Courant number, given the freestream velocity and the size of the cell in x-direction, corresponding to the time step finally used in these simulations, 4.05 × 10−5 s. Not only the wall-resolved scales are important in WR-LES modeling. The typical cell size outside of the airfoil boundary layer must also ensure that at least 80% of the turbulent kinetic energy of the flow is resolved in the free-stream regions. According to [26], this leads to a typical cut-off wavenumber of the LES filter around 𝜅𝐿≅38, or ℓ/𝐿≅0.16, since 𝜅~2𝜋 ℓ ⁄. Consequently, assuming that the integral length scale must be at least one order of magnitude lower than the chord length (L ~ 0.03 m), this leads to a requirement for the cell size around 4.8 mm. Effectively, for the 19M cells discretization over the considered volume domain (approx. 13.7 × 13.7 × 0.16·c3), the average cell size can be estimated in Δcell ∼ (ΔVcell)1/3 = 3.5 mm, well-in-range with the required sizes. Moreover, the time step required to track these scales can be computed as ∆𝑡 ≅  = .  =7.0 × 10𝑠, following a similar rationale than URANS computations where 25 intermediate instants per cycle are usually adopted to capture those large-scale fluctuations. Note that this restriction falls within the current time-step of the modeling. 2.6. Solution Procedure, Convergence and Post-Processing Convergence was guaranteed by monitoring the residual history of the solution, which must be dropped below 10−4 for all the resolved variables. The simulations were run for approximately 44 flowthrough times, based on the freestream velocity and the airfoil chord length, until reaching a statistically steady state. This corresponds to approximately 0.8 s of throughflow time over the airfoil. Seven Figure 3. (a) Computational mesh. (b) Minimum cell sizes in the different coordinates. Regarding the number of nodes in the spanwise direction, the guidelines proposed in references [5,17,30] were followed. Once again, in order to obtain a wall-resolved LES resolution, a dimensionless wall distance ∆ z+ = 10 to 40 is recommended. Taken 30 as a reasonable value, cell size J. Mar. Sci. Eng. 2020,8, 212 6 of 18 was fixed to 0.4 mm, which resulted in 15 cells inside the boundary layer, quite close to the 20-cells recommendation inside the boundary layer according to Sagaut [31]. 2.5. LES Modeling and Resolved Scales Additionally, the selection of an appropriate temporal resolution to resolve the turn-out time of the resolved eddies in the LES scheme is also critical. For wall-resolved Large Eddy Simulations, the time step may be calculated as the ratio between the smallest cell size and the fluctuating velocity u’ in the airfoil walls. Assuming ∆ y+~ 1, the required time step is ∆tLES  ∆y u0=0.0132·10−3 0.51 = 2.6 × 10 −5 s, where u’ has been considered as low as 3% of the bulk velocity, according to the experimental measurements for low angles-of-attack (2.5 ◦ ) in the vicinity of the airfoil wake ( u02/U2= 0.001 ) . This selection is also compatible with the satisfaction of the Courant number, given the freestream velocity and the size of the cell in x-direction, corresponding to the time step finally used in these simulations, 4.05 ×10−5s. Not only the wall-resolved scales are important in WR-LES modeling. The typical cell size outside of the airfoil boundary layer must also ensure that at least 80% of the turbulent kinetic energy of the flow is resolved in the free-stream regions. According to [ 26 ], this leads to a typical cut-offwavenumber of the LES filter around κcL 38, or `c/L 0.16, since κc∼ 2 π/`c . Consequently, assuming that the integral length scale must be at least one order of magnitude lower than the chord length (L~ 0.03 m), this leads to a requirement for the cell size around 4.8 mm. Effectively, for the 19M cells discretization over the considered volume domain (approx. 13.7 × 13.7 × 0.16 · c 3 ), the average cell size can be estimated in ∆cell ~ ( ∆ V cell ) 1/3 =3.5 mm, well-in-range with the required sizes. Moreover, the time step required to track these scales can be computed as ∆tLES 1 25 L U=1 25 0.03 17 = 7.0 × 10 −5s , following a similar rationale than URANS computations where 25 intermediate instants per cycle are usually adopted to capture those large-scale fluctuations. Note that this restriction falls within the current time-step of the modeling. 2.6. Solution Procedure, Convergence and Post-Processing Convergence was guaranteed by monitoring the residual history of the solution, which must be dropped below 10 −4 for all the resolved variables. The simulations were run for approximately 44 flow-through times, based on the freestream velocity and the airfoil chord length, until reaching a statistically steady state. This corresponds to approximately 0.8 s of throughflow time over the airfoil. Seven computers with 4-core i5 processors at 2.67 GHz and 4 Gb DDR3 RAM memory were used for the simulations, which took approximately one month to be completed and to collect the required data for post-processing. For the identification of coherent structures comprising the movement of larger scales in the flow, the Q-criterion has been used to describe vortex-interaction phenomena. This formulation takes advantage of the nature of the big vortices arising in turbulent flows, which despite their chaotic and random nature, can be identified as “fluid regions that maintain some of their properties for a relatively large spatial and/or temporal extent” (also known as “coherent patterns”). This method [ 32 ] defines a vortex as a spatial region where the Euclidean norm of the vorticity tensor dominates that of the rate of strain: Q=1 2ΩijΩij 2−SijSij 2> 0. Using this criterion, vortices in any LES computation can be detected and visualized by rendering iso-surfaces of a given Q threshold. Figure 4shows the baseline iso-surface for Q =5000 s−2in the case of 2.5◦of the angle of attack. J. Mar. Sci. Eng. 2020,8, 212 7 of 18 J. Mar. Sci. Eng. 2020, 8, x FOR PEER REVIEW 7 of 19 computers with 4-core i5 processors at 2.67 GHz and 4 Gb DDR3 RAM memory were used for the simulations, which took approximately one month to be completed and to collect the required data for post-processing. For the identification of coherent structures comprising the movement of larger scales in the flow, the Q-criterion has been used to describe vortex-interaction phenomena. This formulation takes advantage of the nature of the big vortices arising in turbulent flows, which despite their chaotic and random nature, can be identified as “fluid regions that maintain some of their properties for a relatively large spatial and/or temporal extent” (also known as “coherent patterns”). This method [32] defines a vortex as a spatial region where the Euclidean norm of the vorticity tensor dominates that of the rate of strain: 𝑄= 󰇡ΩΩ−SS󰇢>0. Using this criterion, vortices in any LES computation can be detected and visualized by rendering iso-surfaces of a given Q threshold. Figure 4 shows the baseline iso-surface for Q = 5000 s−2 in the case of 2.5° of the angle of attack. Figure 4. Example of vortex identification using the Q-criterion for a 2.5° angle of attack. 2.7. Validation For validation purposes, an FX 63-137 airfoil model with a span of 1.1 m and a chord length of 0.305 m was aerodynamically tested in a closed-loop wind tunnel via hot-wire measurements of the flow field. Figure 5a shows a picture of the test chamber (a cross-sectional area equal to 1.0x1.2 m2) where the FX airfoil has been placed in a vertical arrangement. The different experimental equipment employed to complete the measurements is listed in the figure and includes (1) dual HW probe; (2) BNC connectors; (3) Inclined manometer; (4) CTA anemometer IFA-100; (5) Acquiring card; (6) PC; (7) Anechoic chamber and (8) FX 63-137 airfoil. More details about the measuring devices can be found in [7,25]. The measurements were performed at 10 kHz over 25 s for every position. A convergence study of the mean velocity was made by increasing the number of samples for each measurement point until the results no longer differed. Measurements were made by sweeping two rakes at different streamwise locations: L1, in the airfoil wake (x = 1.108c) and L2, at around 75% of the airfoil chord (x = 0.764c). This last location was chosen because it is a representative position of the high airfoil curvature. Both positions may be seen in Figure 5-right. Wake (L1) and Airfoil (L2) will be used to refer to these positions. In each rake, different positions were measured and for every position, a sample of 250,000 point values was obtained with hot-wire anemometry (see Figure 5). The velocity profile for every angle was built through the time average of the point values in every position of the rake. The airfoil was placed at four different incidence angles: -2.5°, 2.5°, 7.5° and 12.5°, none close to stall, and at a 350,000 Re number. Figure 4. Example of vortex identification using the Q-criterion for a 2.5◦angle of attack. 2.7. Validation For validation purposes, an FX 63-137 airfoil model with a span of 1.1 m and a chord length of 0.305 m was aerodynamically tested in a closed-loop wind tunnel via hot-wire measurements of the flow field. Figure 5a shows a picture of the test chamber (a cross-sectional area equal to 1.0 × 1.2 m 2 ) where the FX airfoil has been placed in a vertical arrangement. The different experimental equipment employed to complete the measurements is listed in the figure and includes (1) dual HW probe; (2) BNC connectors; (3) Inclined manometer; (4) CTA anemometer IFA-100; (5) Acquiring card; (6) PC; (7) Anechoic chamber and (8) FX 63-137 airfoil. More details about the measuring devices can be found in [7,25]. J. Mar. Sci. Eng. 2020, 8, x FOR PEER REVIEW 8 of 19 Figure 5. (a) Experimental setup for hot-wire measurements. (b) Sketch of the measuring database. From these measurements, instantaneous values of in-place velocity, velocity angle, turbulence intensities and integral length scales can be obtained for different angles of attack. Hence, all these experimental results will be later compared to the numerical results, so the computations via LES modeling may be positively validated. In addition, the lift and drag coefficients obtained with the present LES simulation have been compared with the experimental results provided by Selig and McGranaham [14] for a smooth FX 63-137 airfoil, as shown in Figure 6. The solid black lines provide the experimental results measured for Re = 350,000 in the UIUC low-speed subsonic wind tunnel (NREL), at free-stream turbulence levels below 0.1%. The white dots correspond to the present CFD results after time-averaging, obtained when both aerodynamic coefficients are stabilized (typically, after 25–30 flow-through times). The lift coefficient is perfectly matched in the unstalled region by the computations, with maximum deviations in the range of just 2.5%. Only at 12.5°, the LES computations and the experimental results show a slight discrepancy. In the case of the drag coefficient, the differences are remarkable for all the range of angles of attack simulated, being the CD always higher in the simulations. This can be attributed to the difference in the free-stream turbulence level between the NREL wind tunnel (roughly 0.1%) and the turbulence intensity imposed in the LES modeling (0.7%, in resemblance to the wind tunnel shown in Figure 5). In particular, the effect of the free-stream turbulence on the aerodynamic performance of airfoils can be drastic when they are operated at low turbulence levels. Huang and Lee [33] have reported severe drag increments with the increase in freestream turbulence intensity, especially if the value is below 0.45%, which is in correspondence to the present database. The artificial overestimation of the drag coefficient in the computations can thus be perfectly associated with the high free-stream turbulence level employed in the model. Another source of uncertainty can be identified in the momentum method employed by Selig to estimate the drag force over the airfoils indirectly. Depending on the position of the traverse hot-wire that it is measuring the velocity profiles at the wake sections, the two-dimensional theory may lead to significant errors in the determination of the drag coefficient. Note that for the CL, Selig and McGranaham do employ a beam balance to measure the lift force directly. Figure 5. (a) Experimental setup for hot-wire measurements. (b) Sketch of the measuring database. The measurements were performed at 10 kHz over 25 s for every position. A convergence study of the mean velocity was made by increasing the number of samples for each measurement point until the results no longer differed. Measurements were made by sweeping two rakes at different streamwise locations: L1, in the airfoil wake (x=1.108c) and L2, at around 75% of the airfoil chord (x=0.764c). This last location was chosen because it is a representative position of the high airfoil curvature. Both positions may be seen in Figure 5-right. Wake (L1) and Airfoil (L2) will be used to refer to these positions. In each rake, different positions were measured and for every position, a sample of 250,000 point values was obtained with hot-wire anemometry (see Figure 5). The velocity profile for every angle was built through the time average of J. Mar. Sci. Eng. 2020,8, 212 8 of 18 the point values in every position of the rake. The airfoil was placed at four different incidence angles: -2.5◦, 2.5◦, 7.5◦and 12.5◦, none close to stall, and at a 350,000 Re number. From these measurements, instantaneous values of in-place velocity, velocity angle, turbulence intensities and integral length scales can be obtained for different angles of attack. Hence, all these experimental results will be later compared to the numerical results, so the computations via LES modeling may be positively validated. In addition, the lift and drag coefficients obtained with the present LES simulation have been compared with the experimental results provided by Selig and McGranaham [ 14 ] for a smooth FX 63-137 airfoil, as shown in Figure 6. The solid black lines provide the experimental results measured for Re =350,000 in the UIUC low-speed subsonic wind tunnel (NREL), at free-stream turbulence levels below 0.1%. The white dots correspond to the present CFD results after time-averaging, obtained when both aerodynamic coefficients are stabilized (typically, after 25–30 flow-through times). The lift coefficient is perfectly matched in the unstalled region by the computations, with maximum deviations in the range of just 2.5%. Only at 12.5 ◦ , the LES computations and the experimental results show a slight discrepancy. In the case of the drag coefficient, the differences are remarkable for all the range of angles of attack simulated, being the C D always higher in the simulations. This can be attributed to the difference in the free-stream turbulence level between the NREL wind tunnel (roughly 0.1%) and the turbulence intensity imposed in the LES modeling (0.7%, in resemblance to the wind tunnel shown in Figure 5). In particular, the effect of the free-stream turbulence on the aerodynamic performance of airfoils can be drastic when they are operated at low turbulence levels. Huang and Lee [ 33 ] have reported severe drag increments with the increase in freestream turbulence intensity, especially if the value is below 0.45%, which is in correspondence to the present database. The artificial overestimation of the drag coefficient in the computations can thus be perfectly associated with the high free-stream turbulence level employed in the model. Another source of uncertainty can be identified in the momentum method employed by Selig to estimate the drag force over the airfoils indirectly. Depending on the position of the traverse hot-wire that it is measuring the velocity profiles at the wake sections, the two-dimensional theory may lead to significant errors in the determination of the drag coefficient. Note that for the C L , Selig and McGranaham do employ a beam balance to measure the lift force directly. J. Mar. Sci. Eng. 2020, 8, x FOR PEER REVIEW 9 of 19 Figure 6. Comparison of experimental (the data from Selig and McGranaham [13]) and numerical aerodynamic coefficients at Re = 350,000. 3. Results and Discussion This section, presenting the most remarkable results of the simulation, has been divided in five different subsections. Sections 3.1 and 3.2 include results from the velocity fields, that are compared with HW measurements for validation purposes. Following, the other subsections are focused on a detailed description of the boundary layer characteristics over the airfoil. 3.1. Velocity Components and Reynolds Stresses The normal-to-airfoil distribution of the streamwise velocity (x-coordinate) is shown at L1, for the angles of attack tested, in Figure 7. Both experimental and numerical values, normalized by the maximum velocity, are compared for validation purposes, showing a good agreement. Particularly, the wake width is perfectly predicted by the LES modeling, though its deficit is slightly overestimated. This is due to a geometrical penalty of the experimental airfoil in the trailing edge, associated with mechanical deviations. As a consequence, the flow is longer attached to the airfoil in the numerics, resulting in a modeled wake deeper than the real one. The largest discrepancies between experimental and numerical results are found at 12.5°, when the flow is more detached and the three-dimensional effects of the trailing edge are more evident. Figure 6. Comparison of experimental (the data from Selig and McGranaham [ 13 ]) and numerical aerodynamic coefficients at Re =350,000. J. Mar. Sci. Eng. 2020,8, 212 9 of 18 3. Results and Discussion This section, presenting the most remarkable results of the simulation, has been divided in five different subsections. Sections 3.1 and 3.2 include results from the velocity fields, that are compared with HW measurements for validation purposes. Following, the other subsections are focused on a detailed description of the boundary layer characteristics over the airfoil. 3.1. Velocity Components and Reynolds Stresses The normal-to-airfoil distribution of the streamwise velocity (x-coordinate) is shown at L1, for the angles of attack tested, in Figure 7. Both experimental and numerical values, normalized by the maximum velocity, are compared for validation purposes, showing a good agreement. Particularly, the wake width is perfectly predicted by the LES modeling, though its deficit is slightly overestimated. This is due to a geometrical penalty of the experimental airfoil in the trailing edge, associated with mechanical deviations. As a consequence, the flow is longer attached to the airfoil in the numerics, resulting in a modeled wake deeper than the real one. The largest discrepancies between experimental and numerical results are found at 12.5 ◦ , when the flow is more detached and the three-dimensional effects of the trailing edge are more evident. J. Mar. Sci. Eng. 2020, 8, x FOR PEER REVIEW 10 of 19 Figure 7. Comparison between numerical and experimental results for the transversal distributions of streamwise velocity at different angles of attack. Another interesting comparison comes from the different distributions of the Reynolds stresses in the wake region, obtained from the time-averaging of the velocity fluctuations. Figure 8a compares numerical and experimental results including longitudinal (𝑢′𝑈 ⁄), transversal (𝑣′𝑈 ⁄) and crossed (𝑢′𝑣′ 𝑈 ⁄) components of the Reynolds stresses for all the tested angles of attack. In general, simulation results present a similar trend when compared to experimental measurements, exhibiting a double-peak pattern in the longitudinal component with a local minimum that corresponds to the minimum wake velocity. On the other hand, the magnitude of the numerical results for this component shows significant differences, especially at higher angles of attack. In the case of the transversal component, a one single peak characteristic for reduced angles of attack is observed, which is progressively evolving into a double-peak distribution as the angle of attack increases. Finally, it is significant for the crossed component how the sign is switched in the wake center at the minimum velocity, being more evident for positive angles of attack with a positive peak in the pressure side and a negative peak in the suction side. Additionally, it is remarkable that the suction side contributes more than the suction side to the Reynolds stresses in the wake fluid. This is due to the reinforcement of the instabilities in the shear layer of the suction side, associated with the detached flow conditions as the angle of attack increases. Complementarily, Figure 8b reveals the flow streamlines along the airfoil, computed by the numerical simulations for the different angles of attack of the incoming flow. Note that the streamlines have been colored by the local intensity of the turbulent fluctuations of the streamwise velocity. As expected, in the case of a low angle of attack, the flow is smoothly attached to the airfoil geometry (see for instance the suction side at 2.5°). However, for higher incidence angles, the streamlines become unstable, following irregular trajectories like in the case of the suction side at 12.5°. Besides, the longitudinal fluctuations tend to be more concentrated towards the leading edge in the suction side as the angle of attack is enlarged. The opposite trend is observed in the pressure side, with the largest area of instabilities in the case of −2.5°. (the negative angles of attack). This behavior, directly related to the boundary layer, is explained in more detail in the following sections. Figure 7. Comparison between numerical and experimental results for the transversal distributions of streamwise velocity at different angles of attack. Another interesting comparison comes from the different distributions of the Reynolds stresses in the wake region, obtained from the time-averaging of the velocity fluctuations. Figure 8a compares numerical and experimental results including longitudinal u02/U2 , transversal v02/U2 and crossed u0v0/U2 components of the Reynolds stresses for all the tested angles of attack. In general, simulation results present a similar trend when compared to experimental measurements, exhibiting a double-peak pattern in the longitudinal component with a local minimum that corresponds to the minimum wake velocity. On the other hand, the magnitude of the numerical results for this component shows significant differences, especially at higher angles of attack. In the case of the transversal component, a one single peak characteristic for reduced angles of attack is observed, which is progressively evolving into a double-peak distribution as the angle of attack increases. Finally, it is significant for the crossed component how the sign is switched in the wake center at the minimum velocity, being more evident for positive angles of attack with a positive peak in the pressure side and a negative peak in the suction side. Additionally, it is remarkable that the suction side contributes more than the suction side to the J. Mar. Sci. Eng. 2020,8, 212 16 of 18 J. Mar. Sci. Eng. 2020, 8, x FOR PEER REVIEW 17 of 19 Figure 12. (a) Pressure coefficient for different angles of attack. (b) Pressure vectors on the airfoil. 4. Conclusions A wall-resolved LES model of the flow around a typical wind turbine airfoil has been developed and resolved, in order to describe the main features of the turbulent structures and the boundary layer. The numerical results obtained have been validated with hot wire measurements in a wind tunnel. Despite some differences observed and explained for the higher angle of attack, the results obtained from the LES simulations for the best angles of attack have shown overall trends and magnitudes in agreement with the experimental data. The detailed description of the development of the boundary layer over the airfoil provides an insight into the main noise generation mechanism in wind turbines, which is known to be the scattering of the vortical disturbances in the boundary layer into acoustic waves at the airfoil trailing edge. In this case, 2D wave instabilities are observed in both suction and pressure sides, but these perturbations are diffused into a turbulent boundary layer prior to the airfoil trailing edge, so tonal noise components are not expected in the far field noise propagation. The detailed description of this experimentally validated wall-resolved LES model can be useful as a guide to another modeling works of similar features. The results obtained can also be used as input data for the prediction of noise propagation to the far-field using a hybrid aeroacoustic model. Figure 12. (a) Pressure coefficient for different angles of attack. (b) Pressure vectors on the airfoil. 4. Conclusions A wall-resolved LES model of the flow around a typical wind turbine airfoil has been developed and resolved, in order to describe the main features of the turbulent structures and the boundary layer. The numerical results obtained have been validated with hot wire measurements in a wind tunnel. Despite some differences observed and explained for the higher angle of attack, the results obtained from the LES simulations for the best angles of attack have shown overall trends and magnitudes in agreement with the experimental data. The detailed description of the development of the boundary layer over the airfoil provides an insight into the main noise generation mechanism in wind turbines, which is known to be the scattering of the vortical disturbances in the boundary layer into acoustic waves at the airfoil trailing edge. In this case, 2D wave instabilities are observed in both suction and pressure sides, but these perturbations are diffused into a turbulent boundary layer prior to the airfoil trailing edge, so tonal noise components are not expected in the far field noise propagation. The detailed description of this experimentally validated wall-resolved LES model can be useful as a guide to another modeling works of similar features. The results obtained can also be used as input data for the prediction of noise propagation to the far-field using a hybrid aeroacoustic model. J. Mar. Sci. Eng. 2020,8, 212 17 of 18 Author Contributions: Conceptualization, K.M.A.D., J.M.F.O., S.V.-S.; methodology, K.M.A.D., S.V.-S.; software, I.S.-G., J.M.F.O.; validation I.S.-G., K.M.A.D., J.M.F.O.; investigation, I.S.-G.; data curation, I.S.-G.; writing—original draft preparation, I.S.-G., K.M.A.D., J.M.F.O., S.V.-S.; writing—review and editing, K.M.A.D., J.M.F.O., S.V.-S.; visualization, I.S.-G., J.M.F.O.; supervision, project administration and funding acquisition, K.M.A.D., S.V.-S. All authors have read and agreed to the published version of the manuscript. Funding: This work was supported by Projects (1) “Caracterizaci ó n y predicci ó n de la generaci ó n aerodin á mica de ruido en perfiles de turbinas e ó licas”, DPI2011-25419, provided by the Spanish Ministry of Economy and Competitiveness; (2) “Desarrollo y construcci ó n de turbinas e ó licas de eje vertical para entornos urbanos”, ENE2017-89965-P, from the Spanish Ministry of Economy and Business, as well as by both “Severo Ochoa” predoctoral research scholarship provided by the Principality of Asturias. Conflicts of Interest: The authors declare no conflict of interest. References 1. Arcondoulis, A.; Doolan, C.; Zander, A.; Brooks, L. A review of trailing edge noise generated by airfoils at low to moderate Reynolds number. Acoust. Aust. 2010,38, 129–133. 2. Oerlemans, S.; Fisher, M.; Maeder, T.; Kögler, K. Reduction of wind turbine noise using optimized airfoils and trailing-edge serrations. AIAA J. 2009,47, 1470–1481. [CrossRef] 3. Blake, W. Mechanics of Flow Induced Sound and Vibration; Academic Press: New York, NY, USA, 1986. 4. Wagner, S.; Bareiss, R.; Guidati, G. Wind Turbine Noise; Springer: Berlin, Germany, 1996. 5. Tucker, P. Unsteady Computational Fluid Dynamics in Aeronautics; Springer: Berlin, Germany, 2014. 6. Lighthill, M. On sound generated aerodynamically. Part I: General theory. Proc. R. Soc. Lond. 1952 ,211, 564–587. 7. Sol í s-Gallego, I.; Meana-Fern á ndez, A.; Fern á ndez Oro, J.M.; Argüelles D í az, K.M.; Velarde-Suarez, S. LES-based numerical prediction of the trailing edge noise in a small wind turbine airfoil at different angles of attack. Renew. Energy 2018,120, 241–254. [CrossRef] 8. Curle, N. The influence of solid boundaries upon aerodynamic sound. Proc. R. Soc. Lond. 1955 ,231, 505–514. 9. Ffowcs Williams, J.; Hall, L. Aerodynamic sound generation by turbulent flow in the vicinity of a scattering half plane. J. Fluid Mech. 1970,40, 657–670. [CrossRef] 10. Roache, P. Verification and Validation in Computational Science and Engineering; Hermosa Publishers: Alburquerque, NM, USA, 1998. 11. Althaus, D.; Wortmann, F. Stuttgarter Profilkatalog I; Institut für Aerodynamik, Friedr. Vieweg & Sohn: Braunschweig/Wiesbaden, Germany, 1981. 12. Man-powered planes get a new lift. In Popular Science; Popular Science Publishing Co.: New York, NY, USA, 1972. 13. Selig, M.S.; McGranahan, B.D. Wind tunnel aerodynamic tests of six airfoils for use on small wind turbines. In Technical Report, NREL/SR-500-34515; National Renewable Energy Laboratory, U.S. Department of Energy: Golden, CO, USA, 2003. 14. Selig, M.S.; McGranahan, B.D. Wind tunnel aerodynamic tests of six airfoils for use on small wind turbines. J. Sol. Energy Eng. 2004,126, 986–1001. [CrossRef] 15. Piomelli, U. Large Eddy Simulation and Related Techniques; Lecture Series 2006-04; Von Karman Institute for Fluid Dynamics: Brussels, Belgium, 2006. 16. Smagorinsky, J. General circulation experiments with the primitive equations. I. The basic experiment. Mon. Weather Rev. 1963,91, 99–164. [CrossRef] 17. Pope, S.B. Turbulent Flows; Cambridge University Press: London, UK, 2000. 18. Wang, G.; Duchaine, F.; Papadogiannis, D.; Duran, I.; Moreau, S.; Gicquel, L. An overset grid method for large eddy simulation of turbomachinery stages. J. Comput. Phys. 2014,274, 333–355. [CrossRef] 19. Zauner, M.; Sandham, N.; Wheeler, A.; Sandberg, R.D. Linear stability prediction of vortex structures on high pressure turbine blades. Int. J. Turbomach. Propuls. Power. 2017,2, 8. [CrossRef] 20. McMullan, W.; Page, G. Towards large eddy simulation of gas turbine compressors. Prog. Aerosp. Sci. 2012 , 52, 30–47. [CrossRef] 21. Papadogiannis, D.; Duchaine, F.; Gicquel, L.; Wang, G.; Moreau, S. Effects of subgrid scale modeling on the deterministic and stochastic turbulent energetic distribution in large-eddy simulations of a high-pressure turbine stage. J. Turbomach. 2016,138, 091005. [CrossRef] J. Mar. Sci. Eng. 2020,8, 212 18 of 18 22. Meana-Fern á ndez, A.; Fern á ndez Oro, J.M.; Argüelles D í az, K.M.; Velarde-Su á rez, S. Turbulence-model comparison for aerodynamic-performance prediction of a typical vertical-axis wind-turbine airfoil. Energies 2019,12, 488. [CrossRef] 23. Cao, H. Aerodynamics Analysis of Small Horizontal Axis Wind Turbine Blades by Using 2D and 3D CFD Modelling. Master’s Thesis, University of Central Lancashire, Preston, UK, 2011. 24. Mendez, B.; Muñoz, A.; Munduate, X. Study of distributed roughness effect over wind turbine airfoils performance using CFD. In Proceedings of the 33rd Wind Energy Symposium, Kissimmee, FL, USA, 5–9 January 2015. 25. Solis-Gallego, I. Caracterizaci ó n Del Comportamiento Aeroac ú stico De Perfiles De Turbinas E ó licas En Flujo Turbulento. Master’s Thesis, University of Oviedo, Asturias, Spain, 2017. (In Spanish). 26. Tucker, P. Computation of unsteady turbomachinery flows: Part 2—LES and hybrids. Prog. Aerosp. Sci. 2011 , 47, 546–569. [CrossRef] 27. Dahlström, S.; Davidson, L. Large eddy simulation of the flow around an aerospatiale A-aerofoil. In Proceedings of the European Congress on Computational Methods in Applied Sciences and Engineering ECCOMAS 2000, Barcelona, Spain, 1–3 September 2000; pp. 1–20. 28. Davidson, L.; Dahlström, S. Hybrid LES-RANS: An approach to make LES applicable at high Reynolds number. Int. J. Comput. Fluid Dyn. 2005,19, 415–427. [CrossRef] 29. Lastra, M.; Fern á ndez Oro, J.M.; Galdo Vega, M.; Blanco Marigorta, E.; Santolaria Morros, C. Novel design and experimental validation of a contraction nozzle for aerodynamic measurements in a subsonic wind tunnel. J. Wind Eng. Ind. Aerodyn. 2013,118, 35–43. [CrossRef] 30. Chapman, D. Computational aerodynamics development and outlook. AIAA J. 1979 ,17, 1293–1313. [CrossRef] 31. Sagaut, P. Large Eddy Simulation for Incompressible Flows; Springer: Berlin, Germany, 2006. 32. Haller, G. An objective definition of a vortex. J. Fluid Mech. 2005,525, 1–26. [CrossRef] 33. Huang, R.F.; Lee, H.W. Effects of freestream turbulence on wing-surface flow and aerodynamic performance. J. Aircr. 1999,36, 965–972. [CrossRef] 34. Comte-Bellot, G. Hot-wire and hot-film anemometers. In Measurement of Unsteady Fluid Dynamic Phenomena; Richards, B.E., Ed.; Hemisphere: Washington, DC, USA, 1977; pp. 123–162. 35. Drela, M. XFOIL: An analysis and design system for low reynolds number airfoils. In Lecture Notes in Engineering; Mueller, T.J., Ed.; Springer: New York, NY, USA, 1989; Volume 54. 36. Cebeci, T.; Mosinskis, G.; Smith, M. Calculation of separation points in incompressible turbulent flows. J. Aircr. 1972,9, 618–624. [CrossRef] 37. Mayle, R. The 1991 IGTI Scholar Lecture: The role of laminar-turbulent transition in gas turbine engines. J. Turbomach. 1991,113, 509–537. [CrossRef] 38. Winkler, J.; Carolus, T.; Moreau, S. Airfoil trailing edge blowing: Broadband noise prediction from large eddy simulation. In Proceedings of the 15th AIAA/CEAS Aeroacoustics Conference (30th AIAA Aeroacoustics Conference), Miami, FL, USA, 11–13 May 2009. © 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).