On the prospects of improving the numerical analysis of symmetric airfoils at low Reynolds numbers and high angles of attack by a LES approach: the case of NACA0021 profile

. The working conditions of airfoils along modern wind turbine blades are putting new focus on the importance of properly characterizing the aerodynamic performance of different airfoil families also at high angles of attack (AoAs) beyond stall and at Reynolds numbers much lower (from few thousands to one million) than those commonly analyzed before. Several test cases are showing that even higher-order computational methods (like RANS/URANS CFD) are unable to properly capture the complex flow physics taking place past the blades, when deep stall occurs or when the AoA changes so rapidly to provoke the onset of dynamic stall. To fill this gap, the use of high-fidelity methods, like the Large Eddy Simulation (LES) is proposed, even though it implies a massive increase of the calculation cost. In order to analyze the prospects of using LES in comparison to RANS for low Reynolds, high AoAs, this work presents an in-depth study of the NACA 0021 aerodynamics at the Reynolds number of 80,000, by means of both traditional RANS approaches and high-fidelity (LES) simulations using the OpenFOAM suite. The selected airfoil has been showing in fact several issues in the correct characterization of its performance in similar conditions in many recent wind energy applications. The LES approach showed the ability to overcome the limitations of traditional RANS simulations, improving the accuracy of the results and reducing their dispersion thanks to the fact that the flow structures in the separated-flow regions are properly captured. Overall, this work underlines that accurate investigations of the aerodynamic performance of the NACA 0021 at low Reynolds require multiple sensitivity studies when RANS approaches are used, and suggests the use of LES simulations in order to increase the accuracy of estimations, especially when studying the stalled-flow operating conditions of the airfoil.


Introduction
The numerical simulations presented in this work are aimed at providing a reliable evaluation of the NACA 0021 airfoil performance at low Reynolds number. Different numerical setups are proposed to assess the most suitable to solve this kind of flows, focusing on the comparison between RANS and LES approaches. .

RANS Investigation
RANS (Reynolds Averaged Navier-Stokes) simulations are a consolidated tool in many industrial applications. Their application to airfoil aerodynamics for wind energy is more recent, and their accuracy at high AoAs is not completely assessed yet. However, RANS simulations have been considered in the study as the correct starting point for benchmarking the potential of LES in similar applications. In the following paragraphs, the numerical setup used for RANS simulations is described, and their results are presented and discussed.

Numerical Setup
The RANS computations carried out for the study are two-dimensional, as the inflow is confined in the airfoil cross sectional plane. Simulations were run at Reynolds 80,000 with undisturbed flow at atmospheric pressure. Moreover, the flow was considered as incompressible, since its Mach number was extremely low. The gradient terms were discretized by means of a second order central-differencing scheme to improve the solution accuracy. The same scheme was used to discretize the velocity divergence terms, while the other divergence terms are solved with first-second order blended schemes. Simulations were run through an incompressible, bi-dimensional, steady code, which implements the SIMPLE algorithm to solve the pressure-velocity coupling. In external aerodynamics, the boundaries of the numerical domain are loosely defined, but they can affect the quality of the results if they are not placed far enough from the airfoil. For this reason, a domain sensitivity study was performed running simulations at three angles of attack (13°, 15°, 22°) on four different domains, with simulation boundaries located 10, 20, 40, and 50 chords far from the airfoil respectively (these chord numbers will be used in this work to label each domain). A similar wall resolution is maintained in the four domains discretization, changing properly the number of cells and grading along the direction perpendicular to the airfoil symmetry axis. The characteristics of the tested domains are listed in Table 1. The grid refinement effect was evaluated running simulations at different angles of attack with two different setups: the first one is related to the 20-chord domain, while the second one refers to the 40-chord domain using a finer grid. Both the increase in the domain size and the grid refinement improved the accuracy. The finer grid shows an average y + of 0.408, at 13° of incidence, with respect to the 0.670 obtained with the coarser grid. On the other hand, the number of cells increased from 2.1x10 5 to 6.3x10 5 , with a consequent increase in computational time.
Two different turbulence closures were used in the simulations to compare their performance in solving the flow under investigation. More specifically, the k-ω SST model [15,16] was selected as representative of the state-of-the-art of two-equation models, additionally, the transitional behaviour of the flow was investigated coupling this turbulence model with the γ-Reθ transition model [17].

Results
In this section, the main outcomes of the RANS simulations are reported. First, the domain sensitivity analysis is discussed, then the grid sensitivity results are reported, and finally the results obtained using the two turbulence closures are compared.

Domain sensitivity analysis
When the outlet boundaries are far enough from the airfoil, so that the outflow can be supposed uniform, the boundary location does not affect the results. The flow was simulated with an AoA of 13°, 15°, 22° with the four domains using the k-ω SST model. Figure 2 and Figure 3 show the lift coefficient and the drag coefficient as a function of the incidence angle, respectively. A dispersion in the lift curve is apparent only for the simulation with the highest angle of attack, while at low incidence the lift is almost unaffected by the domain size. This discrepancy at high incidence is not remarkable, being lower than 10% for each domain. The drag curve is less affected by the domain size, with the smaller one already able to provide satisfactory results. Both quantities show the expected trend as a function of the angle of attack. Based on the above, the influence of the boundaries location was considered negligible, and the smaller domain sizes were chosen to reduce the computational cost of the simulation.

Mesh sensitivity analisys
In addition to the domain size sensitivity, a mesh sensitivity analysis was performed with the k-ω SST turbulence model. The lift curves obtained running simulations at several angles of attack on two grids are reported in Figure 4. Slight differences can be noted as long as the angle of attack is low, but around the peak, the higher resolution of the finer grid brings about an increase in the predicted lift values. Furthermore, in the finer grid, stall is more prominent at high angles of attack as the predicted loss in performance increases. The lift-to-drag ratio is defined as the efficiency of the airfoil: the higher lift predicted with the finer grid leads to a higher efficiency, but after stall the efficiency decreases for the finer grid as a consequence of the rise in drag. The lift under-prediction near the peak with the coarser grid is a consequence of an overestimation of the separation bubble length due to the anticipated separation onset as underlined by skin friction curves shown in Figure 5.

Turbulence closure sensitivity
The simulations previously discussed have been performed using the k-ω SST turbulence model. In this section, the reliability of this model is compared with the transitional model. For the latter, simulations were run with both the 20-chord and the 10-chord domain size to further investigate the boundary location effects on the calculated flow field.
Upon examination of Figure 6, which shows the lift curves obtained with the two models, it can be seen that the k-ω SST predicts a lower value of the lift at low angle of attack and a higher value at high angle of attack. Moreover, no evidence of stall phenomena can be detected by this model. On the other hand, the transitional model captures a sudden lift drop between 16° and 17°, indicating the occurrence of the aerodynamic stall in this range of angles of attack. A finer angular resolution would be necessary to exactly define the stall angle. The same trend is observed in the drag curves shown in Figure 7. Skin friction curves provide a clear comparison of the separation bubble predicted by the two models. At 13° of AoA (see Figure 8), the γ-Reθ predicts a laminar separation, as shown by the negative friction coefficient, in the front part of the suction side as consequence of the flow acceleration around the leading edge. The flow reattaches as turbulent until a second major bubble occurs just after x/c=0.5. Conversely, the SST model predicts just one separation bubble that occurs upstream at x/c=0.35. The total length of the separated zone is similar for the two cases as, for the transitional flow, the length of the first bubble almost balances the length reduction of the second bubble. Nevertheless, this first recirculation core has a lower impact on the airfoil performance as it is characterized by a lower thickness. Since the separated-flow losses are mainly due to the bubble thickness rather than to its length, it is clear that the γ-Reθ simulations consistently estimate a higher lift.
The differences between the two domain sizes using the transitional model are maximum in the proximity of the stall as the 20-chord domain seems to slightly anticipate the phenomenon. These discrepancies are much less marked than the ones deriving from the turbulence closure, confirming the scarce impact on the results of the boundary location provided it is at least 10 chords far away from the airfoil. At higher angle of attack the airfoil operates with stalled flow conditions: the γ-Reθ model detects an earlier separation point and the two bubbles previously detected merge in a large separated area, as shown in Figure 9. The recirculation core moves towards the front part of the suction side showing no dependency on the domain size. The trailing edge pressure values are similar, but the pressure at the trailing edge, as well as the mean pressure on the pressure side, predicted by the SST model are higher, resulting in a positive contribution to the airfoil load. This one is added to the positive term associated with the delayed separation and a consequently higher acceleration on the suction side. Results obtained with the γ-Reθ model show that with 19° of angle of attack the flow is not able to remain attached: the NACA 0021 behaves like a bluff body and the load capacity is compromised. Figure 10 and Figure 11 emphasize the effect of increasing the angle of attack on the separation bubble captured by the transitional model. When the flow is attached, two distinct separated cores are detected, as the angle of attack increases the leading-edge bubble separation and reattachment points move towards the leading edge and the bubble in the aft part of the airfoil increases its length and thickens. When the stall occurs the two separated cores merge in a single bubble.

LES Investigation
In this section, the results obtained with a Large Eddy Simulation (LES) approach are discussed. This high-fidelity model was able to provide a reliable comparison with RANS results and to achieve a detailed study of separated flow structures and unsteady phenomena. Since high fidelity simulations can be particularly demanding in terms of computational resources requirements, it has been also investigated whether a LES approach is worth being applied in order to get accuracy improvements. The LES simulation was run at three angles of attack (13°, 16° and 19°) to reproduce attached flow, near-stall and stalled-flow operating points.

Numerical setup
The LES approach is intrinsically three-dimensional, for this reason, a layer composed by 67 cells has been added to the 10-chord mesh in the span-wise direction within a thickness of 0.25c, leading to a Δz + of about 15 and a total of 9,849,000 cells. The cube root volume of the cell is selected as the filter for the system of equations and the WALE model [18] is chosen to model the sub-grid turbulent scales because of the good balance between accuracy and computational cost, since it is an algebraic model, no additional transport equation must be solved. A cyclic boundary condition applies for the external faces in span-wise direction, setting the periodicity pattern. To minimize the reflections at the outlet that could negatively influence the solution, an advective velocity boundary condition is specified at the outlet boundary. The high temporal accuracy was obtained with a second order scheme. Gradients and divergence terms were discretized with a second order central-differencing scheme, as in a LES approach, numerical diffusion should be minimized to mitigate its impact on the natural flow oscillations arising from turbulence. In order to speed-up the numerical simulation and to improve its robustness, the RANS solution was used to initialize the flow field at the same incidence angle. The solution was advanced in time with a time step of 4μs for a total simulated time of 0.4s, corresponding to 21 through-flow times along the airfoil. Laminar flow is imposed at the inlet. The solver implements the PIMPLE loop and the stability is ensured even if the maximum Courant number slightly exceeds the unity. Both the RANS and the LES simulations have been run using 64 processors on a Linux cluster of an Intel CPU E5-2680-2.8GHz: each RANS analysis lasted around half an hour, while 4 days were needed for the LES computation.

Results
At 19° of AoA, the airfoil works in stalled conditions. The LES approach allows to study properly the vortex shedding phenomenon and the temporal evolution of the separation bubble. Figure 12, Figure 13, Figure 14 and Figure 15 show the temporal evolution of the zones with negative axial velocity highlighting the presence of separated flow regions. The reported frames cover one vortex shedding period of approximately 0.03s. At the beginning of the period, the recirculation core covers most of the suction side, then it thickens in the aft part of the airfoil until a vortex is shed at almost half a period. As consequence of the vortex shedding, the separation bubble length decreases, and it does not cover the trailing edge area. Later on, the bubble grows again towards the trailing edge and the process is repeated. As the vortex is convected downstream, it is dissipated by the fluid viscosity and numerical diffusion due to the grid coarsening far from the airfoil.
The vortex shedding phenomenon takes place also at the leading edge and it is characterized by smaller spatial scales than the trailing edge ones. The oscillations of the stalled flow affect the temporal evolution of performance parameters such as the lift and the drag shown in Figure 16 and Figure 17. A numerical transient can be appreciated in the convergence history of the monitored quantities up to almost 0.3s, but once a pseudoconvergence is reached, they start oscillating with an almost regular amplitude. The period of these oscillations is the same observed for the vortex shedding: it is worth noting that the lift is minimum when the vortex has just been shed as it impinges on the suction side with a consequent pressure increase. It is harder to clearly identify any regular oscillations in the drag curve, but it seems that its trend is somehow phase-shifted compared to the lift one. Finally, LES results were compared with RANS and experimental ones. From an experimental point of view, the results of Jacobs [13] are taken as a benchmark, since they appeared to be the most reliable ones. Figure 18 shows the lift curves resulting from the different approaches. Numerical results are affected by a high dispersion: the SST model does not provide reliable results, since it underestimates the lift for each angle of attack and provides a nearly flat trend without an evidence of stall phenomena. With a transitional model the experimental trend is well captured even if the airfoil performance are still underestimated. The results obtained with LES are time-averaged in the interval 0.3s-0.4s to remove the numerical transient, the predicted lift coefficient is found to be in better agreement with the experimental measurements. To further outline the separation bubble on the suction side of the airfoil, the Q-criterion contours reported in Figure 19 clearly highlight the turbulent structures captured by the simulation. Figure 20 and Figure 21 show the skin friction curves and the pressure distributions along the airfoil obtained from the three numerical simulations for the same AoA, respectively. Assuming LES results as the most accurate ones, it is clear that the SST model underestimates the length of the separation, moving the separation point downstream at x/c=0.1, while the RANS transitional model accurately catches the separation onset. On the other hand, both RANS models are not able to capture correctly the core of the recirculation, leading to a higher trailing edge pressure than the one predicted by the LES simulation, as apparent upon comparison of the pressure distributions shown in Figure 21.

Conclusions
The aerodynamic performance of the NACA 0021 airfoil available from the open literature have been found to be affected by a huge dispersion. For this reason, the present work was devoted to the numerical investigation of the airfoil performance at low Reynolds number in a wide range of angles of attack, extending also beyond stall.
Since the RANS approach is still widely used to select airfoil geometries for many aerodynamic applications, the ability of this approach to capture the flow field around and isolated airfoil has been assessed. To this end, sensitivity studies have been carried out, focusing on the evaluation of the most important parameters for the RANS numerical setup. From a perusal of the results, one could notice that the size of the computational domain does not affect the results, provided that the boundaries are located more than 10 chords far away from the airfoil. The wall resolution plays an important role as well, but further improvements to bring y + below unity lead to a reduced gain in accuracy with respect to the computational cost increase. To fill this gap, an average y + of around 0.65 has been considered the best trade-off. Numerical results highlighted the key role of the turbulence closure to obtain reliable results. Due to the low value of the Reynolds number (80,000), the SST turbulence model was not able to provide satisfactory results, and the coupling with a transitional model was required to properly capture the flow evolution at different angles of attack. Despite an appreciable improvement, discrepancies were found with the available experimental data in terms of lift polar chart, with an underestimation of the airfoil performance. Adopting a highfidelity approach is thus required to further improve the accuracy of the numerical predictions. LES simulations were run at three angles of attack and the results were focused on the 19° simulation, i.e. a stalled-flow operating point. The separation bubble was found to be correctly captured by LES, as well as the unsteady temporal evolution of the airfoil aerodynamics, strongly affected by vortex shedding phenomena. Although the use of highfidelity techniques is still limited by the need for large computational resources, it was found that only the use of a LES approach allowed one to study accurately low Reynolds stalled flows, since this method only is able to solve both the recirculation cores and the temporal evolution of the flow structures.
Future works could be focused on further improving the accuracy of LES simulations extending the domain in the span-wise direction, refining the grid or extending the total simulated time to pursue a better statistical convergence of the results. The LES approach could be used to study the flow in a wider range of angles of attack to clearly identify the stall position and gain information on the deep stall behavior. The effect on performance of the turbulence intensity should be investigated as well, as it could delay stall phenomena.