Abstract
This work presents the analysis, application, and comparison of thirteen fluid flow models in the prediction of two-dimensional airfoil aerodynamics, considering laminar and turbulent subsonic inflow conditions. Diverse sensitivity analyses of different free parameters (e.g., the domain topology and its discretization, the flow model, and the solution method together with its convergence mechanisms) revealed important effects on the simulations’ outcomes. The NACA 4412 airfoil was considered throughout the work and the computational predictions were compared with experiments conducted under a wide range of Reynolds numbers () and angles-of-attack (). Improvements both in modeling accuracy and processing time were achieved by considering the RS LP-S and the Transition SST turbulence models, and by considering finite volume-based solution methods with preconditioned systems, respectively. The RS LP-S model provided the best lift force predictions due to the adequate modeling of the micro and macro anisotropic turbulence at the airfoil’s surface and at the nearby flow field, which in turn allowed the adequate prediction of stall conditions. The Transition-SST model provided the best drag force predictions due to adequate modeling of the laminar-to-turbulent flow transition and the surface shear stresses. Conclusions, recommendations, and a comprehensive research agenda are presented based on validated computational results.
1. Introduction
The measurement and prediction of aerodynamic forces on two-dimensional airfoils is a problem that has been widely investigated since the early 1930s and its development has produced important improvements in the aerospace, automotive, and wind-based sciences, among others [1–4]. Prior to the experimental assessment of aerodynamic forces, the state-of-the-art procedures [1, 2] impose major prerequisites such as the detailed manufacture of the tested airfoil [5], the setup of expensive wind tunnel facilities [6], and the use of special sensing equipment to characterize both the aerodynamic behavior of the airfoil and the disturbances it produces on the free stream (e.g., streamlines, flow attachment/detachment, flow compression dynamics, and wake aerodynamics). In addition, correction factors [1, 4, 6–10] are often applied to account for nonideal inflow conditions (e.g., buoyancy, solid blockage, wake blockage, or streamline curvature corrections). These prerequisites and issues, together with the overall propagation of uncertainty, turn the experimentation procedures into daunting tasks. A useful, inexpensive, and faster alternative to perform aerodynamic characterizations involves the implementation of computational methods for the theoretical estimation of aerodynamic forces, which are predicted through the numerical solution of the governing equations of fluid mechanics. This approach is formally known as Computational Fluid Dynamics (CFD).
It is generally acknowledged that there is no universal model/method that ultimately describes the complete characteristics of a fluid flow and its interactions with objects with reasonable accuracy while employing a reasonable amount of computational resources. This modeling problem becomes more complex as more physical phenomena are considered (e.g., if turbulent, compressible, and multiphase flows are considered, among other relevant conditions). Therefore, depending on the case study conditions and the assumptions made, different CFD-based approaches with different levels of sophistication can be employed. Some of the most important fluid flow modeling techniques are briefly presented next: (1) the potential flow theory [11], considered the coarsest modeling approach, does not account for turbulence or vorticity effects in its basic formulations. Nonetheless, recent advances in potential flow theory and Boundary Layer (BL) modeling have led to the development of the vortex modeling approach [12], in which viscous and vorticity effects have been successfully integrated into the fluid flow modeling, resulting in improved aerodynamic predictions. (2) The turbulence modeling approach, considered as the industry standard approach for design purposes, is much more complex and computationally demanding since both small-scale and large-scale turbulence effects are modeled by solving either the Reynolds-Averaged Navier-Stokes (RANS) or the Favre-Averaged Navier-Stokes (FANS) equations complemented with turbulence models [13, 14]. (3) Advanced techniques that solve large-scale turbulence effects and model only the small-scale turbulence effects are based on large eddy simulations (LES) complemented with subgrid-scale models [15]. Finally, (4) more advanced techniques based on Direct Numerical Simulations (DNS) [16], which are typically implemented for theoretical research purposes, have the ability to solve the whole range of spatial and temporal scales of the turbulence and predict all the effects and interactions between fluids and solids at the cost of an extraordinary large amount of computational resources.
Currently, there is no consensus regarding the minimum level of flow modeling required to accurately predict airfoil aerodynamics, while considering different flow regimes/states (e.g., laminar, transitional, and turbulent flows) and different inflow conditions (e.g., Reynolds numbers and angles of attack). Despite this issue, in works mainly focused on airfoil shape development and optimization [3, 17] it is often considered sufficient to solve the compressible Euler equations [18] or solve a set of potential flow equations that are coupled with integral BL formulations, some of which are implemented in popular public domain and commercial codes such as XFOIL [19, 20] or VisualFoil [21], respectively, to estimate aerodynamic forces under subsonic, transonic, or supersonic flows. The selection of these approaches, however, was mostly based on the convenient amount of computational resources they require rather than their performance for predicting aerodynamic forces.
Only a limited number of research works have attempted to determine the accuracy of different fluid flow modeling techniques for predicting the aerodynamic behavior of two-dimensional airfoils undergoing different inflow conditions. Wolfe and Ochs [22] presented a laminar/turbulent flow analysis, considering an airflow at a Reynolds number (Re) of and the range of angles of attack (, AOA) , to determine the asymmetric S809 airfoil aerodynamics. They employed the commercial code CFD-ACE, which solves the FANS equations coupled with the Standard - turbulence model. Wolfe and Ochs contrasted the computed pressure coefficient distributions and the computed aerodynamic coefficients with experimental measurements obtained under laminar inflow conditions and observed a drag force overprediction when fully turbulent computations were considered. To address the issue of simulating transitional flows, they developed a mixed laminar/turbulent calculation method, in which the computational domain was split into one laminar and one turbulent region at a guessed transition point, which in turn improved the drag force predictions. They concluded that more research on both the determination of the laminar-to-turbulent flow transition point [23] and the accurate modeling of turbulent effects under stall conditions was necessary to reduce observed discrepancies.
Some of the discrepancies observed by Wolfe and Ochs are related to BL modeling issues. In their work, the modeled dimensionless wall distance ( [24], which is a parameter typically used to determine what sublayers of the BL are solved) was of the order of , thus limiting the probed sample volume of the BL to the logarithmic and outer layers. Therefore, the modeling of the near-wall flow dynamics and the calculation of wall shear stresses depended on the use of wall functions. The standard wall functions [24, 25], such as the ones used by Wolfe and Ochs, have proven to be inaccurate while modeling BLs subject to large adverse pressure gradients (like the ones encountered on airfoils undergoing inflow conditions at large AOA), which in turn induce flow detachment conditions. The appropriate description of the complex BL, from which aerodynamic forces are calculated, requires the accurate modeling of the viscous, turbulent, and rotational properties of the flow found within the airfoil’s vicinity. Therefore, the better the airfoil’s BL is modeled, with special emphasis on the viscous sublayer (or laminar sublayer, which is located at the inner part of the BL), the better the agreement between the computed and the measured aerodynamic forces and the observed flow dynamics. In order to model the viscous sublayer, a solved is required over the entire airfoil surface [24, 26]. The interested reader is directed to [1, 2, 27] for additional and comprehensive descriptions of the physics of airfoil aerodynamics.
Eleni et al. [28] presented a work focused on determining which turbulence model, among the Spalart-Allmaras, the Realizable -, and the SST - turbulence models, was the best performer for predicting the symmetric four-digit National Advisory Committee for Aeronautics (NACA) 0012 airfoil aerodynamics, while considering an airflow at a and for the range of AOA . Similar to the study performed by Wolfe and Ochs, Eleni et al. found drag overpredictions when comparing fully turbulent computations with experimental measurements that considered laminar inflow conditions. To improve the drag force predictions, they conducted mixed laminar/turbulent simulations, similar to the procedure proposed by Wolfe and Ochs. However, in order to determine the laminar-to-turbulent flow transition point, Eleni et al. developed an iterative method that depends on already measured data. In addition, they conducted simulations, considering five different Reynolds numbers (), at a zero AOA by assuming both fully turbulent and mixed laminar/turbulent conditions. In such simulations, the laminar-to-turbulent flow transition took place at the same axial point at both the top and the bottom surfaces of the airfoil due to its symmetry. The resultant drag coefficients were compared satisfactorily with the experimental measurements. They concluded that for both fully turbulent and mixed laminar/turbulent flows the best performer turbulence model was the SST -.
Kumar et al. [29] presented a work in which the asymmetric NACA 4412 airfoil was simulated in a turbulent airflow, at a and for the range of AOA , while considering the Spalart-Allmaras and the Standard - turbulence models. The computed predictions were contrasted with the experimental data (at a ) provided by Abbott and von Doenhoff [1]. Kumar et al. reached the same conclusions as Wolfe and Ochs about the importance of the accurate determination of the laminar-to-turbulent flow transition points at both the upper and the lower surfaces of the airfoil, since a notable overprediction of the drag coefficient was observed by simulating a fully turbulent flow over the airfoil’s vicinity. They concluded that the Standard - model was the best performer turbulence model.
Villalpando et al. [30] presented an assessment of the ability of different turbulence models to predict the asymmetric NACA 63-415 airfoil aerodynamics, while considering an airflow at a , the range of AOA , and an inlet turbulence intensity (TI) of 1%. The contrasted turbulence models were the Spalart-Allmaras, the RNG -, the SST -, and the Reynolds Stress Low-Re S-. The computational results were compared with experimental measurements performed under laminar inflow conditions, which were provided by the Risø National Laboratory for Sustainable Energy. Unlike the observations of the previously described works, the contrasted turbulence models accurately predicted the experimental drag coefficient for . Only at large AOA (i.e., under stall conditions) all the tested models overpredicted both the lift and drag coefficients. Villalpando et al. did not report the numerical convergence criteria of their finite-volume-based solution method but reported oscillatory convergence, an observation pointing to numerical instabilities, for moderately large AOA (e.g., ). In such cases, an averaging procedure was performed to estimate the aerodynamic forces. The approach, however, is believed to be inadequate since such convergence issues are often related to the quality of the domain discretization and/or to the numerical approach used to solve the problem. They concluded that the SST - model was the best performer model, with the Reynolds Stress Low-Re S- being the worst performer model.
Other studies related to the prediction of airfoil aerodynamics (e.g., [31–34]), in which laminar/transitional flows are involved, may be affected by the already noted overprediction of the drag coefficient, which is an issue directly related to the lacking ability of full-turbulent models to simulate the laminar-to-turbulent flow transition. In order to deal with this issue, improved turbulence models, named Transition-based models [24, 26, 35], have been developed to accurately estimate the laminar-to-turbulent flow transition zones and solve the flow as turbulent at downstream locations. So far, improved drag predictions have been obtained as described by Yuhong and Congming [36], Yao et al. [37], Aranake et al. [38], and Khayatzadeh and Nadarajah [26].
Yuhong and Congming [36] presented a work in which the asymmetric S814 airfoil aerodynamics was predicted using the Transition SST model, considering an airflow at a and for the range of AOA . They concluded that the Transition SST model predicts the lift and drag coefficients more accurately than full-turbulent models only in prestall conditions. Under stall conditions, both tested turbulence models (the Transition SST and the SST -) failed to predict the aerodynamic behavior of the tested airfoil.
Yao et al. [37] performed a similar work by predicting the symmetric NACA 0018 airfoil aerodynamics for an airflow at a , a relative Mach number of 0.023, and for the range of AOA . The contrasted turbulence models were the Standard -, the RNG -, the four-equation Transition SST model, and a five-equation Reynolds Stress model. For all the considered inflow conditions, all turbulence models overpredicted the drag coefficient, as concluded in previous works. The magnitude of the lift coefficient was systematically underpredicted for both large negative AOA and large positive AOA. The Reynolds Stress model was the top performer in that study.
Aranake et al. [38] presented an evaluation of the RANS-based Transition model, which was coupled with the Spalart-Allmaras turbulence model (named the Transition model), for predicting the asymmetric S827 and the S809 airfoil aerodynamics. For the S827 airfoil aerodynamics prediction, an airflow at a , a Mach number of 0.1, an inlet , and the range of AOA were considered. The same conditions, but considering a , were assumed for the evaluation of the S809 airfoil aerodynamics. After comparing fully turbulent computations performed with the Spalart-Allmaras model, laminar/transitional computations performed with the Transition model, and experimental measurements conducted under laminar inflow conditions, Aranake et al. found improved lift, drag, and pressure coefficients predictions under both prestall and stall conditions. Nevertheless, significant discrepancies were found at intermediary AOA, where a characteristic double-stall condition arises for both tested airfoils. They concluded that in situations for which the flow remains completely attached or massively separated, the Transition model qualitatively and quantitatively exhibits the same behavior as the fully turbulent SA model. But when moderate flow separation occurs, the coupled model significantly improves the quality of the predictions.
Khayatzadeh and Nadarajah [26] presented a comprehensive evaluation of the RANS-based Transition model, which was coupled with the SST - turbulence model, for predicting the asymmetric NLF(1)-0416 and the S809 airfoil aerodynamics. They provided modifications for both the SST - and the Transition models in order to allow an appropriate gradual activation of the SST - model along the BL upon the onset of transition from laminar-to-turbulent flow conditions. The proposed modifications improved the aerodynamic predictions of the original models, which were compared with experimental measurements performed under laminar inflow conditions. The improved model predicted the laminar-to-turbulent transition locations, the skin friction coefficient distribution, the pressure coefficient distribution, and the lift, drag, and moment coefficients with reasonable accuracy for different inflow conditions, which included the ranges and . Khayatzadeh and Nadarajah highlighted the importance of using an adequate domain discretization (considering average values below 1) and the importance of developing robust and well-calibrated relations for the prediction of the laminar-to-turbulent flow transition, which were shown to be very sensitive to the considered inflow condition.
The evidence presented in the above-described works indicates that by employing Transition-based models, or by performing mixed laminar/turbulent procedures, improved predictions can be obtained when comparing computational results with experiments conducted under laminar/transitional inflow conditions. However, for the mixed laminar/turbulent procedures proposed in [22, 28], three unsolved issues have been identified: (1) no insights of the boundary conditions used between the split regions were provided; thus, the procedure may have resulted in an unphysical representation of the flow field in the vicinity of that boundary, (2) the determination of the transition point was guessed or computed from already measured data, and (3) the laminar/turbulent domain segmentation must be changed for different scenarios, resulting in a complex task when considering asymmetric airfoils undergoing nonhomogeneous inflow conditions. Moreover, grid-independence tests must be performed for each considered scenario. If a mixed laminar/turbulent procedure is to be implemented, instead of guessing the transition location, the authors of the present work recommend computing the airfoil’s top transition point () and bottom transition point () by first using specialized software such as XFOIL that incorporates the method for transition prediction [39], which has proven to be accurate [40], thus avoiding possible metastability issues. It should be noted that one of the key advantages of employing Transition-based turbulence models is that they implicitly solve the three above-described issues, at the expense of increased computational requirements and the incorporation of complex formulations for the accurate prediction of the laminar-to-turbulent flow transition [24, 26, 35].
As it becomes apparent from this literature review, most of the research works have focused on validating the effectiveness of a limited number of fluid flow models for predicting two-dimensional airfoil aerodynamics, while considering a limited range of inflow conditions. None of them provided a justification for the selection of specific full-turbulence models (e.g., while coupling with Transition models [26, 38]) during the validation process and their scope was limited mainly due to the availability of experimental measurements. Moreover, only a few works have investigated the sensitivity of the simulations’ outcomes to the different free parameters, and the majority of the works have not provided detailed descriptions of the numerical solution process and its convergence mechanisms. As consequence, only limited conclusions can be drawn about the performance of the tested models due to the impossibility of extrapolating the findings to situations different from the original case studies. It should be noted that a wide range of conflicting findings arise in the literature with regard to which models, or combination of models, are the most effective in terms of solution quality and computational efficiency. Thus, in order to overcome these issues, this work reports on an extended assessment of the accuracy of thirteen state-of-the-art fluid flow models applied to the problem of quantifying aerodynamic forces on two-dimensional airfoils for a wide range of inflow conditions. The outcomes of the assessment allowed to (1) identify the best performing fluid flow models, (2) understand the modeling pitfalls for different conditions (e.g., stall conditions), (3) determine which free parameters are the most important during the computational evaluation, (4) identify strategies for improved numerical convergence, and (5) identify key research needs and provide a comprehensive research agenda.
A full comparison of two different CFD-based methodologies was performed; the first one is based on the potential flow modeling approach [11, 41] complemented with an transition model and a set of integral BL formulations, which are implemented in the XFOIL 6.96 software. The second methodology is based on the turbulence modeling approach, where twelve different turbulence models were tested. Transition-based turbulence models were adopted when laminar/transitional inflow conditions were considered. For the remaining turbulence models, the free-stream flow was considered to be turbulent (i.e., no attempts were made to separate the laminar and the turbulent regions within the computational domain) and, therefore, an overprediction of the drag coefficient was expected when comparing the computational outputs with experimental measurements that were typically conducted under laminar inflow conditions. The asymmetric four-digit NACA 4412 airfoil was considered as a test case. The rationale for selecting this airfoil is twofold: on the one hand, abundant information of the NACA 4412 airfoil can be found in the literature [1, 42–45], containing experimental measurements for up to 23 different Reynolds numbers ranging from . On the other hand, the asymmetric features of the NACA 4412 airfoil pose a reasonable challenge to the flow modeling techniques.
The remainder of this work is structured as follows: Section 2 presents the different setups of the computational study. Section 3 presents the computational results, their validation, and the corresponding physical interpretations. Finally, Section 4 presents the overall conclusions and outlines future research.
2. Computational Study
The NACA 4412 is an airfoil that has a maximum camber of 4%, which is located at 40% from the leading edge, and has a maximum thickness of 12%, all percentages measured with respect to the airfoil’s chord length. The NACA 4412 ordinates were obtained from [20, 45]. The airfoil’s trailing edge was smoothed, with the aid of the XFOIL 6.96 software, on the last 5% of the chord to produce a sharp closed profile since the original ordinates had an open section at the tip of the trailing edge (i.e., blunt shape). This is a common issue for analytically developed airfoils and affects all the NACA four-digit airfoils.
The present work reproduced the experimental tests reported in [1, 10, 42, 43, 45] in which transonic and supersonic flows were avoided since the experiments were performed in the Langley two-dimensional low-turbulence pressure tunnel [6], which compressed the airflow up to an absolute pressure of 4 atmospheres in order to increase the air density. Therefore, the simulated airflow was considered to be incompressible and standard sea level air properties ( (kg/m3), (kg/m-s), (K), and (kPa)) were considered in all case studies. Furthermore, the earth’s gravitational force was neglected.
The computer used to perform the computational assessment was a customized machine containing a water-cooled Intel Core i7-2600k processor operating at 5.2 GHz, 16 GB of memory DDR3 @1600 MHz, a Corsair Force GT solid state drive operating at 6 Gb/s, and an ASUS Maximus IV GENE-Z motherboard based on the Intel Z68 chipset. Double-precision parallelized simulations were performed for all case studies using an Intel 64 Message Passing Interface (MPI-2).
2.1. Computation of the NACA 4412 Airfoil Aerodynamics Using XFOIL 6.96
XFOIL [19, 20] is Fortran-based software, created by Drela in 1986 at the Massachusetts Institute of Technology (MIT) during the MIT Daedalus project, for the analysis of the subsonic aerodynamics of isolated airfoils. The XFOIL computations are based on the panel method [11, 19, 41, 46], which is combined with an laminar-to-turbulent transition method [35, 39] and a set of integral boundary layer formulations. Given an initial inflow condition, the flow velocity distribution around the airfoil is computed from the panel method while accounting for viscous forces and the induced vorticity from the airfoil surface. The resultant boundary layer and wake are interacted with a surface transpiration model. The resultant flow field is incorporated into the fluid mechanics viscous equations, yielding a nonlinear elliptic system of equations which is solved by a Newton-Raphson algorithm, resulting in both a complete pressure and velocity distributions in the airfoil vicinity. The lift force coefficient () is calculated by direct surface pressure integration, as viscous contributions to the lift force are often neglected, and the pressure coefficient () is calculated using the Karman-Tsien compressibility correction. The drag force coefficient () is determined from the wake momentum thickness at a location far downstream of the airfoil and calculated with use of the Squire-Young formulation. The methods, corrections, and the boundary layer formulation used in XFOIL are extensively described in [19, 46].
For panel-based methods, the first set of free parameters is related to the discretization of the airfoil’s geometry. In all the XFOIL-based simulations, a constant number of panel nodes (160) were considered. The panel nodes were concentrated towards both the leading and the trailing edges of the airfoil, as shown in Figure 1, with the aim of increasing the density of nodes in these sensitive zones. The trailing edge to the leading edge panel density ratio was of 0.15. The panel density ratio at the leading edge was 0.2 and the maximum panel angle was 7.87°. The second set of free parameters is related to the definition of the specific inflow conditions to be studied. In XFOIL, the laminar-to-turbulent flow transition begins when one of the following two scenarios occur: (1) a free transition occurs when the criterion is met or (2) a forced transition occurs when a trip or the airfoil’s trailing edge is encountered. For this Transition-based model, a free parameter named the critical amplification factor (), which affects the laminar-to-turbulent flow transition location, must be defined. A suitable value of this parameter depends on the ambient disturbance, or Turbulence Intensity (TI), in which the airfoil operates and mimics the effect of such disturbances on the flow state transition. A value of = 2.6232 was set to simulate a free-stream TI of 1%, which was the same flow condition of the experimental measurements obtained from the literature. Finally, it was observed that the XFOIL computations were usually very fast, in the order of milliseconds for one simulated case, and consume almost negligible computational resources.
2.2. Computation of the NACA 4412 Airfoil Aerodynamics through Turbulence Modeling
Two-dimensional, incompressible, steady-state, turbulence modeling-based simulations were performed for predicting the NACA 4412 airfoil aerodynamics using the specialized commercial package ANSYS Fluent 13.0. On the preprocessing stage, the geometry, the boundary topology, and the finite volume domain discretization were set/performed as described in next subsections.
2.2.1. Boundary Topology
The first set of free parameters on finite volume-based simulations refers to the definition of the boundary topology, which is strongly dependent on the geometric complexity of the tested object/system. The main concern is related to the requirement of allowing enough free space between the tested object/system and the simulation’s boundary so that no blockage effects occur. Blockage effects may produce significant side effects on the modeled physics, even if numerical convergent simulations are obtained. Different boundary topologies have been considered in the literature within the context of two-dimensional airfoil aerodynamics prediction. Most of them are based on the two-dimensional immersed boundary method [47], which is an approach used for solving external fluid flow problems (i.e., problems in which the fluid flows through the external surface of the tested object). Therefore, the geometry of interest is immersed in a uniform or nonuniform nonboundary-conforming Cartesian grid; thus the flow and its interactions with the tested object are modeled.
In the present work, three different boundary topologies based on the C-type shape, as shown in Figure 2, were studied with the aim of improving the quality of the numerical solutions and reducing the resources spent during the computational assessment. Type 1 topology consisted in an upwind semicircular boundary (denoted by the letters I and J) of radius equal to 10 times the chord length of the airfoil and a downwind rectangle of total height and axial length of 20 and 10 times the chord length of the airfoil, respectively. Type 2 and 3 topologies consisted of upwind semielliptical boundaries with semimajor axis of 15 and 20 times the chord length of the airfoil and a semiminor axis of 10 times the chord length of the airfoil. The boundaries were extended at least 10 times the chord length of the airfoil in all directions in order to avoid flow blockage effects, thus properly encompassing all the relevant zones of the simulation (e.g., the wake zone and the zones where the pressure, velocity, and main turbulence variables distributions still affected the computation of aerodynamic forces). Elliptical-based topologies enabled discretizing the domain using structured quadrilateral elements (e.g., see Section 2.2.2) without significantly affecting the final grid quality and its use can be advantageous when highly asymmetric and/or sharper airfoils are considered. However, elliptical topologies often require more quadrilateral elements and require finer grid refinements towards the airfoil surface in order to obtain acceptable values of dimensionless wall distance (). In the boundary topology analysis performed for the NACA 4412 airfoil, it was found that the use of any of the proposed topologies is suitable as it was easy to get enough grid quality to get quality-converged results, considering any of the studied inflow conditions. However, type 1 topology was used in the rest of this work as it requires fewer discretized elements in the final refined grid. In cases where the grid quality is low due to the geometric features of the tested airfoil, the use of cubic splines or the airfoil shape itself can provide a guideline for the construction of the upwind boundaries (i.e., edges I and J in Figure 2).
2.2.2. Spatial Discretization
The finite volume-based domain discretization consisted of structured quadrilateral elements that were refined towards airfoil’s surface. The use of high-quality structured quadrilateral grids allows numerical algorithms to converge faster, thus requiring less computational resources. The grid quality is commonly measured in terms of the cell aspect ratio, cell Jacobian ratio, cell parallel deviation, cell maximum corner angle, cell skewness, and cell orthogonal quality, among other measures. Each edge of the boundary topology (labeled with a specific letter, as shown in Figure 2) contained a different number of nodal points and different nodal refinements (e.g., arrows in Figure 2 define the refinement direction in which the average quality of the grid is enhanced). In previous studies [28, 48], 80,000 quadrilateral elements formed the total grid and were considered sufficient since successful grid-independence tests were performed, and, therefore, the same grid was used for all the computational experiments. In the present work, it was found that the use of a unique grid was not sufficient to properly test the different fluid flow models and thus several grid-independence tests were performed for each model. This fact can be understood in terms of BL modeling; since each model predicts different results in the converged solution, as described in Section 3, different grids are required to obtain suitable values of the dimensionless wall coordinate (), hence adequately modeling the complete airfoil’s BL. The grids were constructed in such a way that the converged varied along the airfoil’s surface with an approximated maximum value of 1 for cases considering a Reynolds number of and below 1 for cases considering lower Reynolds numbers. The first nodal point was typically located at a distance of m from the airfoil’s surface. The best nodal distributions, determined through grid-independence tests (as discussed in Section 3.1), are summarized for each turbulence model in Table 1. The bias factor, or refinement factor, is defined as the ratio of the largest distance between two adjacent nodes in the discretized edge to the smallest distance between two adjacent nodes in the discretized edge, so the distance between nodes increased linearly towards the nonrefined section of the edge.
2.2.3. Turbulence Models, Numerical Schemes, Boundary Conditions, and Convergence Criteria
In ANSYS Fluent 13.0, the set of RANS equations [13, 14, 24], complemented with turbulence models, are solved. A total of 12 different turbulence models were tested and are summarized in Table 2. A complete description of each turbulence model can be found in [24].
The models’ characteristic constants were maintained at their default values. Enhanced wall functions were considered instead of standard wall functions, although no meaningful differences in flow modeling were observed while considering standard wall functions [24, 49] for cases in which the grid was sufficiently refined towards the airfoil surface, so a value of 1 was obtained. In all cases, a density-based solver with an absolute velocity formulation was considered. The solution method considered an implicit formulation and a Roe flux-difference splitting (FDS) convective flux type scheme. The gradients and derivatives were computed using the least squares cell-based method and the flow variables (i.e., the and velocities, the pressure, and all the turbulence variables) were solved and interpolated using second-order upwind discretization schemes. A full description of such methods and schemes can be found in [24, 25]. The edges A, I, J, and D in Figure 2 were prescribed with an inlet-velocity boundary condition in which the magnitude and direction of the free-stream airflow varied with the Reynolds number and the AOA. The values of the inlet turbulence variables were calculated as a function of the free-stream TI (1%), which was reported in the experimental data found in the literature, and the eddy length scale, which was defined as 1% of the length of the airfoil (0.01 m) [24] since the experiments were conducted in a low turbulence wind tunnel where larger free-stream eddies were hardly found. The edges B and C were prescribed with a pressure-outlet condition, where the gauge pressure was set to 0 Pa and the backflow turbulence intensity and backflow turbulence length scale were set to 1% and 0.01 m, respectively. The edges K and L (the airfoil’s top and bottom surfaces) were prescribed as stationary walls with a no-slip shear condition. The airfoil was considered to be solid and made of aluminum with a density of = 2,719 kg/m3. The wall roughness height and the roughness constant were defined as 0 m and 0.5, respectively, as no other information was available from the literature results. A zero wall roughness height (smooth wall) means that no roughness effects [27, 50] were included during the work development.
The numerical convergence definition is based on the computation of double-precision residuals. In ANSYS Fluent 13.0, the residual of a variable being solved, considering a density-based solver, is defined as the time rate of change of the conserved variable (). Equation (1) shows the unscaled root-mean-square residual definition, which is used for all the variables being solved:where is the number of nodal points in the spatial discretization. In order to properly judge the numerical convergence, the scaled residual definition presented in (2) was considered:where is the unscaled residual of the th iteration and is the largest absolute residual value found within the first five iterations. In the present work, an absolute convergence criterion of ≤ was considered for all the solved variables. This value was considered adequate as no significant improvements in the solution’s precision (e.g., up to 5 significant digits when computing aerodynamic coefficients) were obtained when considering lower convergence criteria values.
In ANSYS Fluent 13.0, the coupled set of governing equations are discretized in time, for both steady-state and unsteady-state computations, and are solved with the use of an explicit or an implicit time-marching algorithm. The present work adopted a Euler-type implicit discretization in time of the governing equations, which was combined with a Newton-type linearization of the flow fluxes, to produce the following preconditioned linearized system in delta form [24]:where the center and off-diagonal coefficient matrices and are given bywhere the residual vector and the time step , computed from the Courant number (), are defined aswhere is the cell volume with surface area A, is the cell face area, is an intermediate solution, is a spatial coordinate in the simulation domain, is a preconditioned matrix, whose aim is to help in reducing numerical divergence issues at the beginning of the simulations [24], is the maximum of the local eigenvalues of the preconditioned system, is a vector containing the convective variables to solve, is a vector that contains the viscous stress tensor terms, and contains source terms such as body forces and energy sources:where , , , , and are the density, velocity, total energy per unit mass, pressure, and heat flux, respectively. In the present work, the energy equation is not solved, since the flow is considered to be subsonic and incompressible, so the last term of and is not considered.
During computations, the variation of the values of the explicit underrelaxation factors (URFs) was crucial to get quality-converged results and to speed up convergence [24, 51]. When oscillatory convergence on the monitored residuals was observed, the URFs values (which typically range around 0 ≤ URFs ≤ 1) were decreased to almost a fifth of its default values and later, after numerical oscillations were reduced, they were increased to almost the unity. Nonetheless, the variation of the CFL number provided the best speed-up improvements. By increasing the CFL number an increased time step is computed; on the one hand, the increased time step develops the steady-state solution faster through the time-marching algorithm and on the other hand makes the first term of (4) less important and therefore more dependent on the computed linear strain rate tensor , thus reducing the effects of the preconditioned system on the Euler-type implicit scheme. From the linear stability theory, it is well known that the implicit formulations are unconditionally stable (e.g., the Courant-Friedrichs-Lewy condition [25]) in contrast to explicit formulations. This fact provides the ability to propose any desired time step to the problem being solved. However, nonlinearities in the governing equations will often limit the numerical stability at the beginning of the simulations (e.g., first 300–500 iterations) and thus the CFL number was maintained low (e.g., 0.1 ≤ CFL ≤ 5) to avoid numerical divergences. Later, the CFL number was slowly increased (e.g., up to a value of 200) until the convergence criterion was met. This simple procedure accelerated the numerical convergence at least three times in contrast with constant CFL number simulations. However, numerical convergence is often difficult unless a good initial solution is provided. In ANSYS Fluent 13.0 three different initialization procedures are available. The first type of initialization (or standard initialization) assigns the same value, which is typically determined from the inlet condition, to all the variables being solved within the discretized domain. The second type of initialization, named the “hybrid initialization” [24, 51], solves the Laplace equation (i.e., is based on a potential flow procedure) in a preprocessing stage to produce a velocity field and a pressure field which smoothly connects high and low pressure values in the computational domain. The third type of initialization consists in incorporating a solution from a previous simulation by using solution interpolations (multigridding approach). The present work considered a multigridding technique in which solutions computed using coarser grids (i.e., having less quadrilateral elements and where the converged values were larger than 30) were used as initial solutions to solve refined cases, which had grids such as the ones reported in Table 1. Nonetheless, the hybrid initialization procedure was considered while conducting simulations using coarser grids. Even after these considerations, in some cases numerical convergence could not be obtained if a proper multigrid initialization and a proper variation of the URFs and the CFL number are not provided/performed.
The computation of aerodynamic forces employing turbulence modeling approaches is very slow in contrast with the XFOIL computations, requiring execution times in the scale of tens of minutes per simulated case and consuming significant computational resources as a function of grid size and considered turbulence model.
2.2.4. Test Cases
In order to properly determine the performance of the tested models while predicting two-dimensional airfoil aerodynamics, more than 1,200 CFD-based simulations, considering the range of AOA and five different Reynolds numbers (, , , , and ), were performed and the outputs were contrasted with experimental measurements reported in the literature. Some of the literature results were not reported in a tabular fashion; thus the data was digitized using the commercial software GetData Graph Digitizer [52]. In addition, most of the drag-based characteristics are commonly reported as a function of the lift coefficient. In order to express the drag coefficient as function of the AOA, the MATLAB 2012a Piecewise Cubic Hermite Interpolating Polynomial (PCHIP) function was used. Finally, in order to compare the CFD-based computations with the experimental measurements provided in the literature, bicubic interpolations were performed using the MATLAB R2012a griddata function, which is based on the Boris Delaunay triangulation method, which in turn fits a surface of the form = and = to the set of experimental data.
3. Results and Discussion
3.1. Grid-Independence Tests
The case chosen for the determination of the minimum amount of quadrilateral elements required to compute grid-independent values of the aerodynamic coefficients was an airflow having the properties described in Section 2, a Reynolds number of , and the AOA 12.22°. The selection of this inflow condition posed a reasonable challenge during the grid-independence characterization due to the required flow modeling complexity, as opposed to previous works [26, 30, 47] in which a prestall condition was considered. The experimental values of the lift and drag coefficients for this test case are and [1, 45]. For this case, XFOIL predicts a and . Figure 3 summarizes the grid-independence test results for all the considered models. From the figure three important findings can be highlighted: (1) it can be noted that the different models differ quite substantially in their predictions of both the lift and the drag coefficients values. (2) At least 120,000 quadrilateral elements are required (with their corresponding refinements towards the airfoil’s surface) to obtain stable values of the computed aerodynamic coefficients and to reach a maximum value of approximately 1 along the airfoil surface. Note, however, that the values are highly dependent on the refinement factor described in Section 2.2.2. (3) In all cases, a drag coefficient overprediction is observed when comparing with the experimental measurement.
3.2. Boundary Layer Modeling Effects
After performing multigridding techniques, as described in Section 2.2.3, insights into how the different solved layers of the BL affect the computation of aerodynamic coefficients were obtained. Figure 4 shows the difference between the absolute values of the aerodynamic coefficients computed using grids with converged maximum values of along the airfoil’s surface and grids with converged maximum values of (i.e., ) for different AOA, Reynolds numbers, and five turbulence models (for the sake of brevity). It can be observed that there is an important effect in both the lift and the drag coefficients predictions as stall conditions are reached. The magnitude of the effect is partially attributed to the change of the wall shear stresses magnitudes, which substantially increase in near-stall and stall conditions. These changes were captured differently in grids that do not solve the inner parts of the BL and depended on the considered wall models [24, 49]. Moreover, the effect varied in magnitude depending on the considered turbulence model; it can be observed that the effect is lower for the Standard - and the Standard - models (note that, for all the - turbulence models, enhanced wall functions were considered for both cases showing converged and ). However, the predictions of other turbulence models (e.g., the Transition-based models and all the Reynolds Stress models) were substantially affected. In particular, the RS LP-S model showed a maximum difference in the lift coefficient magnitude of 0.3275 and a maximum difference in drag coefficient magnitude of 0.0326.
The ability of the Transition-based models to predict the laminar-to-turbulent flow transition is strongly influenced by the converged value [24, 26], and therefore an erratic behavior was obtained when characterizing the effects of not explicitly modeling all the layers of the BL. The side effect of not meeting a value of (e.g., ) is that the transition onset location moves upstream with increasing , resulting in higher drag and lower lift forces predictions. However, it must be noted that this effect can be confused with modeling performance since the comparison between the computed lift coefficients (e.g., considering ) and the experimental measurements performed under laminar inflow conditions tends to improve. It has been observed (not shown) that lift predictions typically improve as larger values are considered, but drag predictions worsen even more.
Figure 5 shows the evolution of the converged maximum along the airfoil’s surface, for both the coarse (e.g., ) and the refined (e.g., ) grids, as a function of the Reynolds number and AOA while considering the RS LP-S model. The converged maximum value increased as both the Reynolds number and the magnitude of the AOA increased. It can be observed that the magnitude of the slope decreases in the range as the maximum wall shear stresses are reached (stall condition). It is expected, as shown in the coarse grid results, that at larger AOA this slope will tend to zero and/or become negative, as no larger wall shear stresses would result in deep stall conditions, consequence of the massive flow separation observed in such conditions.
Not unexpectedly, the variation of the converged maximum value as a function of Reynolds number, for a given AOA, is linear, as shown in Figure 6. It should be noted that the obtained linear regressions are a function of the grid, the considered turbulence model, and the airfoil’s shape. A similar behavior was observed for the results of the remaining turbulence models. The graph provided in Figure 6 is a very useful tool to determine up to which Reynolds number the grid will be able to simulate the viscous sublayer. When the maximum value, evaluated at the first grid layer along the airfoil’s surface, is much larger than 1, the grid will not be able to simulate the viscous sublayer and the resultant airfoil aerodynamics could be questionable, even if the simulations were numerically convergent.
3.3. Validation of the Computational Predictions
It has been experimentally observed [1, 27, 53] that as the Reynolds number increases the lift coefficient will increase but the drag coefficient will decrease for a given AOA. Nonetheless, the net effects of varying the Reynolds number are higher at large AOA. Through the simulation of flows with different free-stream velocities, the ability of the different turbulence models to predict prestall and stall conditions can be quantified. Figure 7 illustrates the Reynolds number effects predictions of each turbulence model on the computed lift-to-drag ratio (direct comparisons of the individual predicted coefficients against experimental measurements performed under laminar inflow conditions [1, 42, 45] can be found in Figures 11, 12, 13, 14, 15, 16, and 17). It can be observed that both Spalart-Allmaras formulations predicted similar results for the lift and drag coefficients as a function of Reynolds number and AOA. The different - and - turbulence models predicted quite different results for both the lift and the drag coefficients, as previously noted from the grid-independence tests. In particular, the Standard -, the Standard -, and the Transition -- models could not capture the lift and the drag forces dependence on the Reynolds number since almost identical values were obtained for each computational experiment. However, the remaining models were able to capture a Reynolds number effect but in different magnitude.
Interestingly, the drag coefficient predictions of the Standard - model are almost the double in magnitude as compared with the SST - model predictions, while considering large Reynolds numbers (e.g., see Figures 15–17). This behavior is similar to the one observed when contrasting the outputs of the Standard - and the Realizable - turbulence models. It should be noted that the drag predictions of the Transition-based models are substantially lower when comparing with the outputs of full-turbulence turbulence models and are in better agreement with the experimental measurements (e.g., see Figures 11–17). This behavior is consistent with the findings observed in previous works [26, 36, 38]. However, despite the improvements in the drag coefficient estimation, the predictions of the lift coefficient by the Transition-based models tend to be larger as the AOA increases and differ when comparing with the outputs of any other full-turbulence model or when comparing with the experimental measurements. In the following sections, separated descriptions and interpretations of this problematic and the accuracy of the remaining flow models on the prediction of aerodynamic forces are provided.
3.3.1. Description of the Accuracy of the Tested Flow Models for Predicting the Lift Coefficient
Based on the results shown in Figures 11 to 17, it can be observed that at the highest simulated Reynolds number (i.e., Re = ) the agreement between the experimental data and all the flow models predictions is relatively good for nonstall conditions. However, as the Reynolds number becomes smaller, the discrepancy becomes larger, except for the RS LP-S model whose predictions remain close to the experimental data for all Reynolds numbers. The RS QP-S model predicted the lowest lift coefficient values on prestall conditions among all the tested models and at all Reynolds numbers, while XFOIL and the Transition -- models predicted the largest lift coefficient values.
When stall conditions were considered, the Transition -- model failed to predict stall at all. Both XFOIL and the Transition SST model performed somewhat better but significantly overpredicted the lift coefficient for AOA beyond the experimental stall angle, similar to the observations provided by Yuhong and Congming [36]. Out of the Standard models, both Spalart-Allmaras models, the NRG -, the Realizable -, the SST -, and the RS Low-Re S- models overpredicted the poststall lift at all Reynolds numbers. For the Standard - and the - models, the poststall lift is somewhat overpredicted at low Reynolds numbers, while for high Reynolds numbers the poststall lift is underpredicted, with good agreement at intermediate Reynolds numbers (e.g., < Re < ). The RS LP-S model comes out the best in lift prediction for all the considered inflow conditions. The remarkable aspect of this model is that it predicts the lift coefficient with reasonable accuracy even in stall states. The RS QP-S and the RS Low-Re S- underpredict and overpredict lift, respectively, at stall conditions and at all Reynolds numbers.
3.3.2. Discussion of the Lift Coefficient Predictions by Different Fluid Flow Models
The prestall lift coefficient overprediction observed at low Reynolds numbers (e.g., see Figures 11–14) can be largely attributed to the inadequate prediction of the laminar-to-turbulent flow transition (i.e., not only close the airfoil surface) and the inadequate turbulence modeling at the outer layers of the BL. Both factors affect the pressure distribution around the airfoil, thus directly affecting the prestall lift force computation, as discussed next. The differences observed while contrasting the prestall lift predictions of the different Transition-based models, which showed the largest discrepancies in the prediction of the experimentally observed prestall lift, are due to their differences while predicting the laminar-to-turbulent flow transition, as exemplified in Figure 8. It can be observed that the laminar-to-turbulent flow transition predictions of XFOIL and the Transition SST model differ quite substantially at low Reynolds numbers. Moreover, an inadequate modeling of turbulent forces in the outer layers of the BL, together with an inadequate laminar-to-turbulent flow transition prediction, can lead to an inadequate modeling of the local flow speed, which in turn affects the pressure distribution around the airfoil, thus affecting the computation of the main component of the lift force. It was observed that Transition-based models typically predict higher flow speeds close to the top surface of the airfoil (e.g., see Figure 21), thus resulting in larger lift forces as compared to full-turbulence models or the experimental data. An extended description of this issue is given in Section 3.4, where a validation study of the predicted pressure coefficient distributions is presented.
As the Reynolds number increased, the airfoil’s superficial laminar-to-turbulent transition locations moved towards the leading edge in a nonlinear fashion, as shown in Figure 8. An increased agreement between the models’ predictions of the laminar-to-turbulent flow transition locations, as well as an increased concordance while predicting the experimental prestall lift force (e.g., see Figures 15–17), was observed as the Reynolds number increased. In such cases, a large part of the airfoil’s BL was turbulent and the effects of the laminar-to-turbulent flow transition were less relevant in the prediction of the lift force since turbulent forces were the most important mechanism regulating the transversal diffusion of momentum between the different layers of the airfoil’s BL, hence affecting the flow speed magnitude in most of the airfoil’s vicinity and, thus, the pressure distribution along the airfoil’s surface. Therefore, the increased agreement between the prestall lift force predictions by all tested models and the experimental measurements at large Reynolds numbers is understood given that the experimental laminar flow became turbulent sooner at the airfoil’s surface (i.e., closer to the leading edge of the airfoil) as the Reynolds number increased. The differences between the predictions of the prestall lift by the different models at large Reynolds numbers essentially depended on the way turbulent forces were modeled. Nonetheless, since all the models’ predictions of the prestall lift were similar in those conditions and showed good agreement with the experimental data, it can be inferred that almost negligible changes on the prestall lift predictions occur at large Reynolds numbers when considering either turbulent inflow conditions (i.e., when employing full-turbulence models) or laminar inflow conditions (i.e., when employing Transition-based models).
When the AOA increases, the laminar-to-turbulent flow transition location moves towards the leading edge (i.e., upstream) at the airfoil’s top surface whereas at airfoil’s bottom surface it moves towards the trailing edge (i.e., downstream), as shown in Figure 8. Moreover, opposite effects occur when the AOA is decreased to negative values. As the absolute value of the AOA increases, a stall condition is reached independently of the considered Reynolds number. Nonetheless, the AOA at which the stall condition occurs depend on the Reynolds number and the geometric complexity of the airfoil. Under stall conditions, almost half of the experimental airfoil’s BL was turbulent, and the flow detachment condition occurred on the turbulent side (e.g., see Figures 21–24). Therefore, the flow detachment modeling is strongly dependent on the modeled turbulent effects. Following this line of reasoning, the most important deviations in the prediction of the lift force by the different turbulence models were found in stall conditions.
The overprediction of the lift force under stall conditions by full-turbulence models can be understood since free-stream turbulence promotes the earlier formation of turbulent BLs, which typically prevent flow detachment conditions (i.e., surface flow attachment/reattachment is enhanced in turbulent BL). Therefore, the modeled delay on the onset of stall was a consequence of the modeled turbulent BL, as opposed to the partial laminar BL of the experimental data. However, the simulated free-stream turbulence intensity was rather low (1%) and the magnitude of the computed delay on the onset of stall seems to be unrealistic and it is not expected in real free-stream turbulent experiments, which unfortunately were not available in the literature. In this regard, the RS LP-S model properly regulates the delay on the onset of stall, as shown in Figures 11–17, due to the proper modeling of anisotropic turbulent forces. Moreover, since the RS LP-S model considered a fully turbulent inflow condition, the effects of the laminar-to-turbulent flow transition and the modeling of a turbulent BL at the bottom surface of the airfoil played a minor role in the prediction of the poststall lift force.
It should be noted that the performance of the Standard - and the Standard - models for predicting lift coefficients can be good at certain Reynolds numbers, but the agreement can be shown to be fortuitous. As illustrated in Figures 7 and 11–17, their predictions are literally independent of the Reynolds number, whereas the experimental lift and drag coefficients increased and decreased with increasing Reynolds numbers, respectively. Therefore, the observed agreement can be regarded as an occasional coincidence. This observation contradicts with the findings of Kumar et al. [29] who, for a fixed Reynolds number of , concluded that the Standard - model provided the best predictions of the NACA 4412 airfoil aerodynamics.
3.3.3. Description of the Accuracy of the Tested Flow Models for Predicting the Drag Coefficient
The drag coefficient was overpredicted by all the RANS-based approaches, being the Transition-based models the best performers at all Reynolds numbers. The best drag predictions occurred at low Reynolds numbers and at low absolute AOA, with the discrepancy monotonically increasing towards higher Reynolds numbers and higher AOA, as shown in Figures 11 to 17.
In all cases, the Standard - and the Standard - models predicted the largest drag coefficients. Similar predictions were found for the SST - and the RS Low-Re S- models; both models perform poorly at lower Reynolds numbers and when stall conditions were considered. All Reynolds stress models overpredicted the experimental drag coefficient values with the RS Low-Re S- predicting the lowest drag values among the three different models. Both the RS LP-S and the RS QP-S predicted similar drag values at all Reynolds numbers and AOA. XFOIL comes out best, quite accurately predicting the drag for all the considered inflow conditions.
3.3.4. Discussion of the Drag Coefficient Predictions by Different Fluid Flow Models
In prestall conditions, the friction-based drag is the dominant component of the total drag force; thus any increase in the magnitude of the wall shear stresses will significantly affect the total drag force. Figure 9 shows the drag components predicted by one full-turbulence model (the RS LP-S) and two Transition-based models (the Transition SST and XFOIL) at a Reynolds number of . The friction-based drag prediction of the RS LP-S model is significantly larger when compared to the other models’ predictions at any given AOA. This is because the RS LP-S model predicts a fully turbulent BL around the airfoil, similarly to other full-turbulence models. Therefore, the overprediction of the prestall drag force by the different full-turbulence models can be attributed to the modeled small-scale turbulence found near the surface of the airfoil (i.e., the airfoil’s full-turbulent BL), which produces larger wall shear stresses and therefore larger values of the skin friction coefficient (), as shown in Figure 25. As described in [1, 27], it is well known that free-stream turbulence increases the total drag force and in some cases doubles the drag force compared to laminar cases. Therefore, in prestall conditions, Transition-based models do a better job in predicting drag because they consider a laminar inflow condition, predict a laminar-to-turbulent flow transition, and solve the flow as turbulent in the airfoil’s downstream vicinity and in the resultant wake. From Figure 9, it can be observed that important discrepancies between XFOIL and the Transition SST model predictions of the laminar-to-turbulent flow transition occur only at prestall conditions, which in turn result in the observed discrepancy of both models for predicting the same values of the friction-based drag. The Transition SST model predicted a lower friction-based drag component (consequence of predicting a delayed laminar-to-turbulent flow transition), resulting in a slight underprediction of the drag coefficient at low AOA when compared with the XFOIL predictions or the experimental measurements, as shown in Figures 11 to 17.
As the Reynolds number increases, the laminar-to-turbulent flow transition points at both the top and the bottom surfaces of the airfoil move towards the leading edge, as shown in Figure 8. However, as opposed to intuition, the relocation of the laminar-to-turbulent flow transition locations does not lead to larger prestall total drag coefficients. Figure 10 shows the evolution of the predicted drag components by various models, while considering different Reynolds numbers. It can be observed that both the friction-based drag coefficient and the pressure-based drag coefficient decreased as the Reynolds number increased and hence the total drag coefficient decreased. This fact can be understood since the rate of increase of the net drag force was lower compared to the rate of increase of the free-stream dynamic pressure at different Reynolds numbers. Moreover, the observed reduction of the pressure-based drag was also a consequence of the modeled turbulence, which prevented the formation of adverse pressure gradients in prestall conditions. The predictions of the pressure-based drag at low AOA were similar among all the different models.
The discrepancies observed in prestall conditions between the Transition-based models’ predictions and the experimental measurements, as the Reynolds number increased, were due to the modeling of the pressure-based drag, which was heavily influenced by the AOA and the modeled turbulence effects, as shown in Figure 10 and as discussed in detail in Section 3.4. Equivalently, the differences between the full-turbulence models’ predictions of the prestall drag were due to the modeling of the pressure-based drag, since almost no changes in the friction-based drag values were found when varying the AOA, as shown in Figures 9 and 10.
When stall conditions were reached, the pressure-based drag coefficient became the dominant component of the total drag coefficient (e.g., see Figure 9) and the friction-based drag coefficient decreased due to the flow detachment conditions. Thus, the growth rate of the total drag force increased dramatically in stall conditions. As it can be observed from Figures 11 to 17, the drag discrepancies between the full-turbulence models’ predictions and the experimental measurements increased more rapidly when the AOA increased, as compared to the discrepancies observed between the Transition-based models and the experimental measurements. This fact implies that the drag force is very sensitive to the modeled turbulence effects. Therefore, under stall conditions, the laminar-to-turbulent flow transition becomes less relevant and the proper modeling of the turbulent effects and the flow detachment/reattachment conditions becomes crucial in the prediction of the total drag force.
The large discrepancies of the drag coefficient predictions by the different full-turbulence models, observed at large AOA, depended on the modeled medium-to-large-scale turbulence, which induced larger values of the pressure-based drag (), which in turn produced a premature flow detachment, as discussed in detail in Section 3.4. Moreover, their predictions were influenced by the lift force overprediction described above, which created an induced drag component (or lift-induced drag ), consequence of the modification on the pressure distribution around the airfoil due to the trailing vortex system that accompanies the lift generation.
Based on the findings presented before, the common modeling pitfalls are related to the inadequate modeling of the laminar-to-turbulent flow transition, the turbulent effects, and the flow detachment conditions. Moreover, there appears to be no current methodology that yields accurate results for both aerodynamic forces for any given inflow condition. On the one hand, the turbulence modeling approach (with the RSM LP-S model as the top performer) is capable of predicting the experimental lift force with reasonable accuracy but tends to overpredict drag force, as the flow is considered turbulent in the entire spatial domain. On the other hand, Transition-based models (as exemplified by XFOIL or the Transition SST model) are able to accurately predict the experimental drag characteristics but tend to overpredict lift forces. The main difference between the RANS-based full-turbulence models and the RANS-based Reynolds stress models is the way in which Reynolds stresses are approximated (e.g., the Boussinesq isotropic eddy viscosity assumption is considered in most --based models [32]). Moreover, the difference between the Reynolds stress models is the way they model the pressure-strain term in the equations for the transport of the Reynolds stresses (), which models transport due to pressure and the mean turbulent strain rate interactions [24]. Since the best lift force predictions were achieved by the RS LP-S model in contrast with the predictions of both the RS QP-S and the RS Low-Re S- models, it can be inferred that the pressure-strain term plays an important role in the prediction of lift coefficient in stall conditions.
In order to enhance the level of understanding of the effects that the modeled turbulence produces on the airfoil aerodynamics, the next section contrasts the main flow variables solved within the airfoil’s vicinity.
3.4. Main Differences between the Flow Models’ Predictions
In order to understand what mechanisms are driving the modeling procedures, a comparison of the main variables solved (i.e., the local flow speed, the pressure coefficient distributions, and the turbulent variables) was conducted. The test cases used for such comparison are (1) an airflow at a Re = and = 0° and (2) an airflow at a Re = and = 20°. Based on the observed agreements and disagreements with the experimental measurements, the contrasted models were XFOIL, the Standard -, the SST -, the Transition SST, and the RS LP-S models.
3.4.1. Test Case 1 (Re = and = 0°)
The first test case was reasonably well predicted by all models in terms of the lift coefficient. In this case, XFOIL and the Transition SST model slightly overpredict the lift coefficient while the RS LP-S slightly underpredicts it. The best lift coefficient prediction is given by the Standard - and the SST - models. The drag coefficient is accurately predicted by XFOIL and the Transition SST models. The worst drag predictions are given by the Standard - and the RS LP-S models (e.g., see Figure 15).
As expected, the predicted pressure coefficient profiles along the airfoil’s surface are very similar among all the considered models and are similar to the experimental measurements, as shown in Figure 18. In this case, both XFOIL and the Transition SST model underpredict the pressure coefficient at all chord positions at the top surface and conversely at the bottom surface. This is consistent with the already noted lift coefficient overprediction. The transition point predicted by XFOIL differs slightly from the one predicted by the Transition SST model, as noted in Figures 8 and 9, where XFOIL predicts the transition point sooner than the Transition SST model. The pressure coefficient shape profiles for both the Standard - and the SST - models are almost identical. At the top surface and in most of the bottom surface, the best performer model is the RS LP-S. The disagreement between the full-turbulence models’ predictions and the experimental measurements is mostly observed within the range 0 ≤ ≤ 0.40 at the top surface and is understood given that the experimental flow was laminar in that section. An increased agreement between the full-turbulence models’ predictions and the experimental measurements is observed after , as the experimental flow became turbulent. Furthermore, an increased agreement between all the flow models’ predictions and the experimental measurements is observed at both surfaces close to the trailing edge of the airfoil.
3.4.2. Test Case 2 (Re = and = 20°)
The second test case is challenging to solve because it considers a stall condition. The best lift force predictors are the RS LP-S and the Standard - models; however, they overpredict the drag force when compared with the remaining flow models’ predictions. From Figure 19, it can be observed that the predictions of the pressure coefficient distribution along the top surface of the airfoil are substantially different between flow models. Interestingly, the pressure coefficient distributions from XFOIL and the Transition SST turbulence model are very similar in both surfaces. However, XFOIL and the Transition SST models overpredicted the absolute value of the pressure coefficient at most of the airfoil’s top surface. This is consistent with the fact that both models predicted similar values for the lift and drag coefficients, as shown in Figures 9 and 15. It is conspicuous that the Standard - model predicted large values of the pressure coefficient (i.e., lower absolute values of the pressure coefficient) along the leading edge of the airfoil, which led to the already observed discrepancy in the predicted aerodynamic forces, as shown in Figures 10 and 15. The overprediction of the pressure coefficient by the Standard - model increased with increasing Reynolds number, as shown in Figure 20, which presents the resultant pressure coefficient distributions for the case of simulating an airflow at a Re = and at an = 20°. This behavior was also found in the predictions of the Standard - and NRG - models. The best performer turbulence model is, again, the RS LP-S showing the best agreement with the experimental pressure coefficient distributions.
3.4.3. Discussion of the Pressure Coefficient Predictions by Different Flow Models
The pressure coefficient and the lift coefficient overprediction of XFOIL has been noted previously [40, 54] and it is attributed to (1) the inadequate prediction of the local flow speed within the laminar section, close to the airfoil surface, and (2) the lacking ability of XFOIL to simulate the turbulent BL development, rotational effects, and the proper flow attachment/detachment conditions. Equivalently, the Transition SST model typically overpredicts the local flow speed at the airfoil’s top surface, close to leading edge, compared to other models that predict more accurate aerodynamics (e.g., see Figures 21 and 23).
The overprediction of the pressure coefficient (or the underprediction of the absolute value of the pressure coefficient) by full-turbulence models (e.g., Standard -) is a direct consequence of the modeled turbulent forces found near the leading edge of the airfoil, as shown in Figures 22 and 24 which illustrate the contours of turbulent kinetic energy, which in turn are limited to the range of values found for the same situation but considering the RS LP-S model (i.e., the uncolored zones represent zones where the magnitude of the turbulent kinetic energy is larger or lower than the maximum or minimum values found for the same situation but considering the RS LP-S model). The large turbulent kinetic energy and the large turbulent viscosity (not shown) predicted by the Standard - model within the airfoil’s vicinity produce an enhanced transversal diffusion of momentum, mostly noted at the upper surface of the airfoil, resulting in both a substantial reduction of the flow speed magnitude and an affectation of the resultant flow direction, as shown in Figures 21 and 23, when compared to the RS LP-S, SST -, and the Transition SST models' results. The induced disturbances on the wind field, product of the large modeled turbulent forces, lead to premature flow detachment conditions, thus producing thicker and larger wakes and promoting the prediction of large values of pressure-based drag coefficient. This is observed phenomenon is consistent with the already studied drag predictions of each full-turbulence model, as shown in Figure 10. Therefore, the prediction large pressure-based drag coefficient values by most full-turbulence models, with the exception of the Transition models, can be partially attributed to the large modeled turbulent forces, which produced premature flow detachment conditions. Nonetheless, proper premature flow detachment conditions are required for the adequate prediction of the lift reduction that accompanies the stall regime, which is well described by the RS LP-S model.
The prediction of large turbulent forces within the airfoil’s vicinity is mostly observed when employing turbulence models that are based on the Boussinesq isotropic eddy viscosity hypothesis [55]. Nonetheless, this modeling issue does not affect the predictions of the Reynolds stress-based models as they solve additional transport equations for the six independent Reynolds stresses and the eddy viscosity isotropic assumption is avoided, thus denoting the main reason why the turbulent forces are better described by the Reynolds Stress-based models compared to the Standard - or the Standard - models.
In order to improve the level of understanding of the already noted drag coefficient overpredictions, an illustration of the skin friction coefficient () is shown in Figure 25. It is worth noting that in the frontal top and bottom surfaces of the airfoil the skin friction coefficient is larger when considering full-turbulence models. As previously discussed, these larger values are expected since the airfoil’s BL is turbulent. Both XFOIL and the Transition SST models predicted lower values of the skin friction coefficient since the flow was laminar in that section. When the laminar-to-turbulent flow transition occurred, the computed increased and became comparable to the predicted values by full-turbulence models. Interestingly, there is reasonable agreement between the predictions of XFOIL, the Transition SST, and the RS LP-S models of the at the airfoil’s top surface. It can be observed that the skin friction coefficient tends to decay rapidly at the top forward half of the airfoil surface, when considering the case = 20°. The observed reduction is attributed to a flow detachment condition. Interestingly, as already discussed, full-turbulence models predict total flow detachment sooner than XFOIL or the Transition SST models, the Standard - being the one that predicts more turbulent kinetic energy in the airfoil’s vicinity and the one that predicts total flow detachment sooner. Equivalent to , the friction-based drag tends decay due to the flow detachment experienced at larger AOA, as previously shown in Figure 9.
Figure 26 shows the magnitude and the shape profile of the solved independent superficial Reynolds stresses, considering the RS LP-S model. In 2D cases, four independent stresses are solved, corresponding to the normal stresses and the off-diagonal Reynolds stresses . It should be noted from this figure that the normal stresses are, in magnitude, more important than the off-diagonal Reynolds stresses at both surfaces of the airfoil and for both studied cases. However, the off-diagonal Reynolds stress cannot be neglected. For both test cases ( = 0° and = 20°) the superficial Reynolds stresses tend to decay rapidly towards the tip of the airfoil, with the exception of the tip zone at the airfoil’s bottom surface for the = 20° case, where the formation of vortexes takes place. Furthermore, the Reynolds stresses shape profiles are significantly affected by flow detachment conditions, as noted from the case = 20°. It was observed that the magnitude of all the Reynolds stresses decayed rapidly in the near wake as a result of the turbulent dissipation rate found in such zone. Finally, as mentioned before, the modeling of such independent stresses is important to get accurate predictions of the lift coefficient in the stall regime, thus denoting the importance of modeling anisotropic turbulent flows during the airfoil aerodynamics prediction.
4. Conclusions and Future Research
This work presented a comprehensive assessment of the accuracy of thirteen state-of-the-art fluid flow models applied to the problem of predicting aerodynamic forces on two-dimensional airfoils. The primary aims of the work were to (1) extend previous literature conclusions by considering additional flow models and a wide range of inflow conditions, (2) identify the best performer models and their key modeling features, (3) study the effects of different free parameters on the simulations’ outcomes, (4) provide a complete set of physical interpretations based on validated computational results, and (5) provide recommendations for future work regarding the modeling of two-dimensional airfoil aerodynamics. The specific improvements achieved during the work development are summarized next.(i)One of the main concerns in previous literature works results is related to the quality of the numerical solutions. During the development of the present computational assessment, it was possible to achieve an adequate numerical convergence, while considering any of the tested turbulence models and even while considering stall conditions. This was possible due to the proper setup of the computational experiments, which considered model-tailored spatial domain discretizations. In this regard, extensive guidelines of the procedures and calculations were provided to allow the reader to replicate them in different scenarios. Furthermore, additional recommendations, based on a boundary topology analysis and on a set of grid-independence tests, were given with the aim of improving the quality and the computational efficiency of the solution process. Through the extensive study of the effects of different free parameters, it was possible to spare computational resources without compromising the modeling accuracy. The grid-independence tests revealed that at least 120,000 elements, with their corresponding refinements towards the airfoil’s surface, are required to simulate stable airfoil aerodynamics at large Reynolds numbers and at large AOA. Moreover, among the most relevant strategies for improved numerical convergence, it was found that the use of finite volume-based solution methods incorporating preconditioned systems, together with the proper variation of the underrelaxation factors and the Courant number values, was crucial to significantly speed up the numerical convergence process.(ii)Contradictory findings were distinguished in the literature with regard to which models, or combination of models, are the most effective in terms of solution quality and computational efficiency while predicting two-dimensional airfoil aerodynamics. Through the development and validation of a variety of computational experiments that encompassed a wide range of inflow conditions, the best modeling approaches were objectively identified. However, it was found that the prediction of aerodynamic forces considering laminar inflow conditions is complex and there is still no methodology that by itself provides accurate predictions of both the lift and the drag forces. On the one hand, the turbulence modeling approach (with the RSM LP-S model as the top performer) was able to predict lift forces that agreed with the experimental results, even when stall conditions were considered, but drag forces were substantially overpredicted. By assuming a free-stream turbulent inflow condition, it was demonstrated that the frictional (or viscous) component of the drag was larger in the frontal zone of the airfoil, since the airfoil’s BL was turbulent, as opposed to the partially laminar BL layer found in the experimental measurements. In addition, premature flow detachment conditions were observed when considering free-stream turbulent inflow conditions in contrast with considering free-stream laminar inflow conditions, which were modeled with the aid of XFOIL and the Transition SST models. The premature flow detachment conditions predicted by full-turbulence models were a consequence of the large turbulent forces predicted at the leading edge of the airfoil, which produced a substantial modification of the local flow speed magnitude and its direction within the airfoil’s vicinity, which in turn produced larger values of the pressure-based drag due to the induced adverse pressure gradients. It was observed that as the magnitude of the turbulent kinetic energy found close to the airfoil’s surface increased the absolute value of the pressure coefficient decreased. This modeled effect increased the drag force predictions discrepancy between the full-turbulence models predictions and the experimental results. On the other hand, XFOIL and the Transition SST turbulence models predicted near identical drag forces as compared to the literature results but overpredicted lift forces even in prestall conditions. This is because both XFOIL and the Transition SST models predicted large flow speed magnitudes within the airfoil’s vicinity, thus affecting the prediction of the pressure coefficient at all chord positions. The overprediction of the flow speed magnitude by XFOIL is believed to be due its lacking ability to properly simulate turbulent BLs, rotational effects, and flow detachment conditions. Furthermore, the overprediction of the flow speed magnitude by the Transition SST model is related to the prediction of the laminar-to-turbulent flow transition in the outer layers of the BL (i.e., not only close to the airfoil’s surface). Important differences were found when contrasting the laminar-to-turbulent flow transition predictions of XFOIL and the Transition SST models in prestall conditions.(iii)The worst performance was observed when full-turbulence models based on the Boussinesq isotropic eddy viscosity hypothesis were employed, which include all the --based models, as they were unable to model a Reynolds number dependence while predicting both aerodynamic coefficients. Moreover, they predicted the largest drag coefficients, which differed with the literature results at all AOA and Reynolds numbers. For these full-turbulence models, large values of the turbulent kinetic energy and the turbulent viscosity were found in the airfoil’s vicinity, which led to a greatly underpredicted absolute value of the pressure coefficient.(iv)During the computational assessment, it was found that significant variations in all the model’s predictions occur by not explicitly modeling the inner layers of the airfoil’s BL. The proper modeling of the complete airfoil’s BL leads to the accurate estimation of wall shear stresses and flow detachment conditions, which are required to match the experimental results. The modeling of anisotropic turbulence employing RANS-based approaches complemented with Reynolds stress turbulence models, considering a linear pressure-strain term in the equations for the transport of the Reynolds stresses, was crucial to obtain accurate predictions of the lift coefficient in the stall regime since the prediction of flow detachment conditions is heavily influenced by the modeled turbulent diffusion of momentum across the airfoil’s vicinity. However, detailed experimental studies on flow detachment/reattachment and its relation with the pressure-based drag coefficient and the free-stream turbulence intensity are required in order to improve the lift predictions in the stall regime.
Based on the validated computational results, the most important modeling pitfalls, which still provide room for improvement and define fertile areas of inquiry, as well as the main research needs, were identified and are presented in summary form.(i)Since the performance of the RS LP-S and the Transition-based models for predicting lift and drag forces was good, respectively, coupling the Transition model with the RS LP-S model could result in a model capable of predicting accurate airfoil aerodynamics when considering laminar/transitional inflow conditions.(ii)Important differences in the laminar-to-turbulent flow transition prediction were found when contrasting the outputs of XFOIL and the outputs of the Transition SST model. A recalibration of the Transition SST model constants could enhance the performance of the model while predicting prestall aerodynamics at different Reynolds numbers.(iii)One important unaddressed question is related to quantification of the effects that different free-stream turbulent conditions have on the prediction of aerodynamic forces (i.e., the determination of the effects of varying the inlet turbulence intensity and the eddy length scale on the aerodynamic predictions), while considering different inflow conditions (e.g., different Reynolds numbers and AOA). In many real applications (e.g., aircrafts flying within the earth’s BL, automobiles or wake-affected wind turbines, etc.), free-stream turbulent conditions are experienced and their aerodynamics should be studied assuming fully turbulent conditions, while considering different turbulence intensities and eddy length scales. In order to validate the outputs of the cases that considered free-stream turbulent inflow conditions, it is required to perform wind tunnel tests over the NACA 4412 airfoil considering free-stream turbulent flows containing different levels of turbulent intensity. Such experimental information is not available in the literature.(iv)Other relevant cases to be studied are (1) the computation of the two-dimensional airfoil aerodynamics up to , (2) assessment of the ability of different turbulence models for predicting transonic and supersonic aerodynamics, (3) assessment of the accuracy of advanced techniques based on LES and different subgrid-scale models for the prediction of two-dimensional airfoil aerodynamics, (4) evaluation of the ability of the Transition-based models for predicting airfoil aerodynamics at very low Reynolds numbers (e.g., the prediction of laminar separation bubbles), and (5) assessment of the net effects the superficial roughness produces on the lift and drag forces prediction.Although RANS-based and FANS-based approaches are often used for comparison purposes during the development of alternative flow models, it is clear that these approaches still exhibit deficiencies while predicting two-dimensional airfoil aerodynamics. The findings of this work, together with the provided physical interpretations, are believed to be useful and can be considered as a building block for the development of alternative flow models.
Nomenclature
| : | Cell face area (m2) |
| : | Airfoil chord (m) |
| : | Drag coefficient (—) |
| : | Friction-based drag coefficient (—) |
| : | Pressure-based drag coefficient (—) |
| : | Skin friction coefficient (—) |
| : | Lift coefficient (—) |
| : | Pressure coefficient (—) |
| : | Critical amplification factor for the transition prediction method (—) |
| : | Pressure (kPa) |
| Re: | Reynolds number [—] |
| : | Temperature (K) |
| TI: | Turbulence intensity (%) |
| : | Reynolds stress (m2/s2) |
| : | Friction velocity (m/s) |
| : | Local kinematic viscosity (m2/s) |
| : | Velocity (m/s) |
| : | Cell volume (m3) |
| : | Axial distance from the airfoil leading edge (m) |
| : | Normalized transition point at the top surface of the airfoil (—) |
| : | Normalized transition point at the bottom surface of the airfoil (—) |
| : | Distance to the nearest wall (m) |
| : | Dimensionless wall distance (—) |
| : | Angle of attack () |
| : | Dynamic viscosity (kg/m-s) |
| : | Fluid density (kg/m3) |
| : | Airfoil density (kg/m3) |
| : | Shear stress (Pa) |
| : | Difference of absolute drag coefficient (—) |
| : | Difference of absolute lift coefficient (—) |
| : | Time step (s). |
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
Support from the Research Chair for Wind Energy at Instituto Tecnológico y de Estudios Superiores de Monterrey de Monterrey (ITESM), under Contract CAT158, and support from the Mexican Council for Science and Technology (CONACYT) are gratefully acknowledged. Special thanks are due to José Santiago Martínez Torres for his valuable aid while conducting computational experiments.