Argyris G Panaras*
Aerospace Engineering Consultant, Athens, Greece
Received Date: May 25, 2015 Accepted Date: June 03, 2015 Published Date: June 18, 2015
Citation: Panaras AG (2015) Numerical Simulation of Airfoil Flow at High Angle of Attack. J Vortex Sci Technol 2:116. doi:10.4172/2090-8369.1000116
Copyright: © 2015 Panaras AG. This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
Visit for more related articles at Fluid Mechanics: Open Access
A NACA 4412 airfoil tested in wind tunnel at flow conditions close to maximum lift is used for testing the accuracy of various turbulence models. What makes this test case unique is the appearance of a stable separation vortex on the upper surface of the airfoil, near the trailing edge. The linear k-ω turbulence model, a non-linear Explicit Algebraic Stress Model and a modified version of the algebraic Baldwin-Lomax model are tested using the same grid and input file. It is found that the tested turbulence models capture the physics of unsteady separated flow. However, the size and shape of the predicted separation vortex are different for each tested turbulence model. Also, good agreement between computational and experimental surface pressures is observed. The results support the view that such a simple configuration is appropriate to be used as benchmark for validating turbulence and LES models.
Computational Fluid Dynamics; Vortex flows; Turbulence models
Computational Fluid Dynamics (CFD) has undergone a remarkable development over the last three decades and is used routinely, together with experimental techniques, in the design of flight vehicles. Secondorder Reynolds Average Navier Stokes (RANS) codes are mainly used, while higher order codes (LES and DNS) gradually enter the field, as the power of computers increases. Presently, the application of higher order codes is restricted to rather low Reynolds numbers. In general, the discretization of the flow equations adds numerical dissipation, which diminishes as the order of a scheme increases. Practically, second order codes achieve accepted accuracy, in cases of optimized aerodynamic shapes, where the flow is attached or mildly separated. However, second order codes are not appropriate for simulating unsteady vortices, due to their numerical dissipation. Recourse to higher order schemes or use of a locally very fine mesh improves the diffusion problem. For example, for the simulation of aircraft wakes the application of RANS codes is restricted around the aircraft, while the roll-up of the calculated vorticity at a downstream cross section and the formation of the trailing vortices are computed by DNS or LES codes or vortex methods.
The performance of the turbulence models which are used for closing the RANS equations is another factor, which affects the accuracy of second-order CFD codes adversely. It is common experience that different turbulence models result in different predictions, when applied to a particular complex flow. Jameson  has stated a generally accepted fact: “it is doubtful whether a universally valid turbulence model, capable of describing all complex flows, could be devised.” Most of the turbulence models are based on the Boussinesq hypothesis, according to which the apparent turbulent shear stresses are related linearly to the rate of mean strain through an apparent scalar turbulent or “eddy” viscosity coefficient, μt. However, in strongly separated flows, the actual dependence of the modeled turbulent shear stresses to the mean strain is non-linear. For alleviating this problem, various nonlinear corrections have been proposed. In general, non-linear models perform better than linear ones.
There is continuous progress towards the improvement of the accuracy of CFD codes, particularly in turbulence modeling. To validate emerging new concepts experimental data around simple configurations are used. In the present article, the case of a NACA 4412 airfoil tested in wind tunnel at flow conditions close to maximum lift will be used to test the accuracy of various turbulence models. This configuration has been experimentally tested by Coles and Wadcock  and is used for the validation of turbulent models for many decades. What makes this test case unique is the appearance of a stable separation vortex on the upper surface of the airfoil, near the trailing edge. This test case is included in NASA’s Turbulence Modeling Resource web site.
Coles and Wadcock , in an attempt to describe the trailingedge separation process on an airfoil operating near maximum lift, performed hot-wire measurements in the boundary layer, the separated region, and the near wake for a flow past an NACA 4412 airfoil at M=0.09 and α=13.87 degrees. The Reynolds number based on chord was 1,520,000. Special care was taken to achieve a two-dimensional mean flow. The main instrumentation was a flying hot wire; that is, a hot-wire probe mounted on the end of a rotating arm. The tests were performed at the GALCIT 10-ft wind tunnel. According to the authors, the main conclusion from the Reynolds-stress data is that the separation process is relatively regular up to the trailing edge of the airfoil. The real challenge to understanding lies in the merging process for the two shear layers just downstream of the trailing edge and in the subsequent rapid relaxation toward the final state of a conventional wake far downstream. Figure 1 is a display of mean-velocity vectors and contours of at the separated region near the trailing edge of the model.
In the experiments, the upper and lower boundary layers were tripped (2.5%c upper surface and 10.3%c lower surface). However, in CFD simulations fully turbulent computations are performed. Also the calculations are performed here on grids with a farfield outer boundary extending to 20c, but the experiment was in a relatively small wind tunnel, which may have had some effect. Coles and Wadcock  provide surface Cp and normalized velocity field data. It is important to note that the experimental velocity data were non-dimensionalized with respect to a non-traditional velocity at a location only about 1 chord below and behind the airfoil. This is different from the freestream velocity value. For posting velocity profile results in the NASA’s Turbulence Modeling Resource, some authors do an ad hoc correction. In the present article, no velocity-profile comparisons are included.
The CFD code ISAAC developed by Morrison  is used in this study. ISAAC is a second-order, upwind, finite-volume method where advection terms in the mean and turbulence equations are solved by using Roe’s approximate Riemann solver coupled with the MUSCL scheme. Viscous terms are calculated with a central difference approximation. Mean and turbulence equations are solved coupled by using an implicit spatially split, diagonalized approximate factorization solver. Because of the high order spatial discretization of the MUSCL scheme, ISAAC uses flux limiters, to avoid spurious oscillations which would otherwise occur in shock waves, discontinuities or sharp changes in the solution domain. The user is able to switch on or off a selected limiter. ISAAC has been developed to test a large range of turbulence models. Algebraic models, various k-ε and k-ω formulations, and Reynolds stress transport models are included in ISAAC.
In the present article, the NACA 4412 flow is calculated with the k-ω turbulence model of Wilcox , the Explicit Algebraic Reynolds Stress Model (EASM) of Rumsey and Gatski  and a modification of the algebraic turbulence model of Baldwin and Lomax  done by the present author, in order to improve its accuracy in separated flows.
It is known that the functional form of the Navier-Stokes equations is the same for laminar and turbulent flows. In the latter case, however, time average has been applied to equations, since in turbulent flows, each flow or thermodynamic parameter has a mean value and a random turbulent fluctuation (for example: ). The time averaged flow equations are known as Reynolds Averaged Navier- Stokes equations (RANS). In them, new apparent shear stresses, known as Reynolds stresses: have appeared, which need to be calculated. The calculation of the Reynolds stresses is not easy; transport equation for each of them must be solved, increasing the total CPU cost. Alternatively, the Boussinesq hypothesis is applied, according to which the apparent turbulent shear stresses are related to the rate of mean strain, through an apparent scalar turbulent or “eddy” viscosity coefficient, μt. For the calculation of the turbulent viscosity coefficient,various turbulence models have been developed, in which the Reynolds stresses and other terms of turbulent fluctuation parameters are related to mean values of the flow: The Boussinesq equation is,
where k is the turbulent kinetic energy and Sij is the strain rate tensor given by,
where ω is the specific dissipation rate.
The turbulent kinetic energy and specific dissipation rate are calculated by solving transport equations similar to those that express the flow parameters .
The explicit algebraic stress model (EASM) of Rumsey and Gatski  replaces the linear Boussinesq approximation with a non-linear relationship of the strain rate and rotation rate tensors,
where Wij is the rotation rate tensor,
The eddy viscosity is given by,
The Baldwin-Lomax  model is a two-layer algebraic zeroequation model. The eddy viscosity coefficient is given by algebraic equations. In the inner layer it follows the Prandlt-Van Driest formulation:
where κ is the von Karman constant (equal to 0.41), D is the van Driest damping factor, Ω is the absolute value of the vorticity and η is the distance normal to the wall. The damping factor is, where is the wall shear stress and νw the wall kinematic viscosity.
In the outer layer the following equations are used:
The quantity Fmax is the maximum value of the moment of vorticity:
The parameter ηmax is the value of η at which F(η) , equation (12), is maximum.
The Klebanoff intermittency factor is:
The quantity udif is the difference between maximum and minimum velocity in the velocity profile. The thickness of the boundary layer is defined by: The constants appearing in the previous relations are: Ccp=1.6, Cwk=0.25, CKleb=0.3.
The wake function (Fwake) is equal to F1 in attached flows and equal to F2 in separated ones. The present author has found that the accuracy of Baldwin-Lomax model is increased in separated flows if the ηmax in equation (10) is replaced by a reference value, ηref. The value of ηref differs from flow to flow and it is based on existing semi-empirical relations that define the boundary layer parameters of an equivalent flat plate flow of the same Reynolds number.
For the estimation of ηref along a flat-plate flow, the semi-empirical analysis of Falkner  is used, which is valid for a Reynolds number between 105 and 1010. According to Falkner  the boundary layer growth along a flat plate is,
If this equation is combined with the relation: then
Defining ηref at a characteristic separation length of each examined configuration, Lsep, the final calculation scheme is,
where ηref has been non-dimensionalized by the length of the body, L.
Then equation (10) is replaced by,
For aerospace configurations we propose the equality: Lsep=L. But if there is extensive crossflow separation, as in the case of slender bodies at incidence, more accurate results are obtained if the characteristic length is equal to a crossflow length (Lsep=d, for axisymmetric bodies). This assumption is reasonable, since in high-alpha flows, the flowfield is dominated by the separated crossflow. The described above modification of the Baldwin-Lomax turbulence model is applied in flows with extensive crossflow separation in .
The examined test case is also included in NASA’s Turbulence Modeling Resource (TMR) as well as in the User’s Manual of ISAAC. For running ISAAC, the input file and the grid provided by Morrison were used. This is a low speed flow: M=0.09, Rec=1.52×106, α=13.87°. At these flow conditions a small separation vortex is formed on the upper surface of the airfoil, near the trailing edge. Some turbulence models are not able to predict the vortex. The grid consists of 257 points around the airfoil and 81 points in the normal to the surface direction (257×81). No grid refinement study has been performed, but comparisons of the present results with similar ones included in TMR and run with much finer grids (897×257), indicated that for the k-ω and EASM models the results were similar. Calculations were performed by using the k-ω and EASM models, as well as the previously described modification to the Baldwin-Lomax model by Panaras. The experiments indicate that the formed separation vortex stays constantly at the same position (Figure 2d). Thus, a steady-state calculation procedure is sufficient for the simulation of the flow. However, for completeness, the steady-state calculations were continued by time-accurate runs, by employing the τ-ts scheme of Rumsey et al. . Furthermore, some runs employed the τ-ts scheme solely, assuming free-stream conditions for the initialization of the flow.
In Figure 2 the standing separation vortices predicted by the tested turbulence models are compared to the experimental evidence, given by NASA at the Turbulence Modeling Resource web site. It is observed in Figure 2a that the k-ω model predicts a very small separation vortex, quite different from the experimental one posted by the TMR page curator. On the contrary, the EASM model predicts a separation vortex of size comparable to that of the experimental one, although its shape is different and it is located further upstream (Figure 2b). Actually, the separation vortex has an inclined delta shape. Also, the upstream boundary layer, which separates and folds around the vortex, bends and forms an abrupt U-turn between the vortex and the shear layer of the lower surface. At this point we mention that the shape of the prediction of the k-ε-SST turbulence model shown in TMR has very similar shape and behavior. The separation region predicted by the algebraic model has a shape similar to that given by the EASM model, but within the U-turn of the shear-layer a smaller elongated vortex has appeared (Figure 2c). The small vortex is generated by the shear layer of the lower airfoil surface and it is counter-rotating. Probably the algebraic model predicts less turbulent flow than the EASM model. However, we note that the two regions of Reynolds stress concentrations shown in Figure 1b give the impression of existence of two vortices. Also, the distribution of the velocity vectors in Figure 1a supports the existence of an elongated vortex in area B, or of a turning shear layer, as that predicted by EASM (Figure 2b). More accurate experiments are needed. As regards the surface pressure, it is observed in Figure 3 that the Cp distribution of the EASM and of the modified Baldwin-Lomax model is closer to the experimental evidence than that of the baseline Baldwin-Lomax and of the k-ω model.
Since flux limiters introduce numerical dissipation and the presently examined flow is incompressible, the calculations shown in Figure 2 were performed with the flux limiter off. To test their effect, the calculations were repeated with limiter on. The results are shown in Figure 4. It is seen in this figure that the flux limiter introduces a significant dissipation effect. The k-ω model does not predict separation at all. The vortex of the EASM has shrunk, and the second vortex of the algebraic model has disappeared. A comparison of Figures 2 and 3 leads to the conclusion that indeed flux limiters introduce considerable numerical dissipation. The predicted flows are more turbulent than those given by the same calculation code, but with limiters off.
In TMR results for the k- ω and EASM models are posted. The separation vortex predicted by the k-ω model has shape and size similar to that shown in Figure 2a. As for the EASM prediction, a separation vortex similar to that shown in Figure 4c is posted, but the enveloping shear layer is wavy, indicating a lack of convergence. The curator of TRM mentions that: “for this particular case the EASMko2003-S model does not converge readily to a steady-state result when using this code (CFL3D) on this refined grid (897x257). However, when run timeaccurately, the solution settles down and becomes reasonably steady (quasi-steady) with only very small oscillations in drag coefficient. Note that these are compressible code results at "essentially incompressible" conditions of M=0.09. There may be a very small influence of compressibility”.
The presented results indicate that the flow around an airfoil at high-angle-of-attack is quite appropriate for validating turbulence models. At conditions close to the maximum lift, a large separation vortex is formed on the upper surface of the airfoil close to its trailing edge. According to Coles and Wadcock , who many years ago studied experimentally the examined configuration: “the real challenge to understanding lies in the merging process for the two shear layers just downstream of the trailing edge”. Actually, neither their pioneered measurements by using a flying hot wire, nor accurate simulations on fine meshes, based on the top-rated turbulence models and posted in NASA’s Turbulence Modeling Resource web site, have answered this challenge. The measured velocity field is not sufficiently accurate downstream of the separation vortex. Also, the evaluated turbulence models predict small or large separation vortices, of oval or of inclined delta shape. But each prediction is different. However, one may confidently argue that non-linear turbulence models have higher accuracy than linear ones.
According to the present results, if the linear Boussinesq relation in the k-ω model is replaced by the non-linear equation of Rumsey and Gatski , the formed non-linear model (EASM) gives improved results comparable to those posted in TMR by the top-rated SST and SA models. In addition an improved version of the classical algebraic zero-equation turbulence model of Baldwin and Lomax  model is tested in the present study. The modified B-L model is very robust and appropriate for the simulation of extensively separated flows, like those generated around flight vehicles flying at high-incidence and in components of supersonic/hypersonic air vehicles, when swept-shock waves generated on their surfaces interact with the surface boundary layers . The application of the modified algebraic model to the flow around the NACA 4412 airfoil, introduced an additional source of uncertainty. A small elongated and counter-rotating vortex, generated by the shear layer of the lower airfoil surface, appears into the U-turn region of the separation shear-layer. Since some turbulence models predict the U-turn of the shear layer, but not the small vortex, it is quite probable that the algebraic model predicts a less turbulent flow than the simulated one.
Before closing we wish to add a comment on the convergence of the numerical calculations. Since the separation vortex does not convect downstream but stays at the same position, the flow is steady. Thus, a steady-state calculation procedure is sufficient for the numerical simulation of the flow. However, because of the effect of the involved turbulence model, it is not known a priori whether a particular simulation leads to a steady flow. Hence, the involvement of a timeaccurate calculation procedure is more appropriate. Furthermore, in the course of this study we discovered that a time-accurate scheme leads smoothly and faster to the converged solution. An example for the calculations that involved the EASM is given in Figure 5. It is seen in this figure that the time-accurate calculation converges fast and smoothly. Of course, each iteration of the time-accurate scheme includes 15 sub-iterations. Still the convergence is faster compared to the steady-state runs. The final solutions of the two calculation schemes were found to be identical.
The author wishes to express his sincere thanks to Professor Spyros Voutsinas of NTUA for his permission to use the computational facilities of Fluid Section for computing the flows shown in this article.