See discussions, stats, and author profiles for this publication at: [Link]

net/publication/283341589

Numerical Simulation of Airfoil Flow at High Angle of Attack

Article · January 2015

DOI: 10.4172/2090-8369.1000116

CITATION READS

1 992

1 author:

Argyris G. Panaras
Independent Researcher
49 PUBLICATIONS 1,088 CITATIONS

SEE PROFILE

All content following this page was uploaded by Argyris G. Panaras on 20 January 2016.

The user has requested enhancement of the downloaded file.

 Sc
ien
ce and Te
ch Panaras, J Vortex Sci Technol 2015, 2:2

no
Vortex Science and Technology [Link]
x
Vorte

ISSN: 2090-8369 logy

Research Article Open Access

Numerical Simulation of Airfoil Flow at High Angle of Attack

Argyris G Panaras*
Aerospace Engineering Consultant, Athens, Greece

Abstract
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.

Keywords: Computational Fluid Dynamics; Vortex flows; Turbulence will be used to test the accuracy of various turbulence models. This
models configuration has been experimentally tested by Coles and Wadcock
[2] and is used for the validation of turbulent models for many decades.
Introduction What makes this test case unique is the appearance of a stable separation
Computational Fluid Dynamics (CFD) has undergone a remarkable vortex on the upper surface of the airfoil, near the trailing edge. This
development over the last three decades and is used routinely, together test case is included in NASA’s Turbulence Modeling Resource web
with experimental techniques, in the design of flight vehicles. Second- site.
order Reynolds Average Navier Stokes (RANS) codes are mainly used, Description of the Experiments
while higher order codes (LES and DNS) gradually enter the field, as
the power of computers increases. Presently, the application of higher Coles and Wadcock [2], in an attempt to describe the trailing-
order codes is restricted to rather low Reynolds numbers. In general, edge separation process on an airfoil operating near maximum lift,
the discretization of the flow equations adds numerical dissipation, performed hot-wire measurements in the boundary layer, the separated
which diminishes as the order of a scheme increases. Practically, region, and the near wake for a flow past an NACA 4412 airfoil at
second order codes achieve accepted accuracy, in cases of optimized M=0.09 and α=13.87 degrees. The Reynolds number based on chord
aerodynamic shapes, where the flow is attached or mildly separated. was 1,520,000. Special care was taken to achieve a two-dimensional
However, second order codes are not appropriate for simulating mean flow. The main instrumentation was a flying hot wire; that is,
unsteady vortices, due to their numerical dissipation. Recourse to a hot-wire probe mounted on the end of a rotating arm. The tests
higher order schemes or use of a locally very fine mesh improves the were performed at the GALCIT 10-ft wind tunnel. According to the
diffusion problem. For example, for the simulation of aircraft wakes the authors, the main conclusion from the Reynolds-stress data is that the
application of RANS codes is restricted around the aircraft, while the separation process is relatively regular up to the trailing edge of the
roll-up of the calculated vorticity at a downstream cross section and the airfoil. The real challenge to understanding lies in the merging process
formation of the trailing vortices are computed by DNS or LES codes for the two shear layers just downstream of the trailing edge and in
or vortex methods. the subsequent rapid relaxation toward the final state of a conventional
wake far downstream. Figure 1 is a display of mean-velocity vectors
The performance of the turbulence models which are used for
and contours of (u1' u2' ) at the separated region near the trailing edge
closing the RANS equations is another factor, which affects the accuracy
of the model.
of second-order CFD codes adversely. It is common experience that
different turbulence models result in different predictions, when In the experiments, the upper and lower boundary layers were
applied to a particular complex flow. Jameson [1] has stated a generally tripped (2.5%c upper surface and 10.3%c lower surface). However, in
accepted fact: “it is doubtful whether a universally valid turbulence CFD simulations fully turbulent computations are performed. Also
model, capable of describing all complex flows, could be devised.” Most the calculations are performed here on grids with a farfield outer
of the turbulence models are based on the Boussinesq hypothesis, boundary extending to 20c, but the experiment was in a relatively small
according to which the apparent turbulent shear stresses are related wind tunnel, which may have had some effect. Coles and Wadcock [2]
linearly to the rate of mean strain through an apparent scalar turbulent provide surface Cp and normalized velocity field data. It is important
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 non- *Corresponding author: Argyris G Panaras, Aerospace Engineering Consultant,
linear corrections have been proposed. In general, non-linear models Athens, Greece, E-mail: [Link]@[Link]
perform better than linear ones. Received May 25, 2015; Accepted June 03, 2015; Published June 18, 2015
There is continuous progress towards the improvement of the Citation: Panaras AG (2015) Numerical Simulation of Airfoil Flow at High Angle of
accuracy of CFD codes, particularly in turbulence modeling. To Attack. J Vortex Sci Technol 2: 116. doi:10.4172/2090-8369.1000116
validate emerging new concepts experimental data around simple Copyright: © 2015 Panaras AG. This is an open-access article distributed under
configurations are used. In the present article, the case of a NACA 4412 the terms of the Creative Commons Attribution License, which permits unrestricted
airfoil tested in wind tunnel at flow conditions close to maximum lift use, distribution, and reproduction in any medium, provided the original author and
source are credited.

J Vortex Sci Technol

ISSN: 2090-8369 VST, an open access journal Volume 2 • Issue 2 • 1000116
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

Page 2 of 5

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: ui , T , ρ . The Boussinesq equation is,
1 2
− ρτ ij = µt (Sij − Skkδ ij ) − ρ kδ ij (1)
3 3
where k is the turbulent kinetic energy and Sij is the strain rate tensor
given by,
1 ∂ui ∂u j
=Sij ( + ) (2)
2 ∂x j ∂xi
The k-ω turbulence model of Wilcox [4] is a linear eddy viscosity
model, which includes equation (1) in its formulation. The eddy
viscosity, μt, is related to mean and turbulent quantities by,
ρ k (3)
µt = Cµ
' '
Figure 1: Experimental mean velocity vectors and contours of (u1u2 ) in the ω
vicinity of trailing edge. where ω is the specific dissipation rate.
The turbulent kinetic energy and specific dissipation rate are
to note that the experimental velocity data were non-dimensionalized calculated by solving transport equations similar to those that express
with respect to a non-traditional velocity at a location only about 1 the flow parameters [4].
chord below and behind the airfoil. This is different from the freestream
velocity value. For posting velocity profile results in the NASA’s The explicit algebraic stress model (EASM) of Rumsey and Gatski
Turbulence Modeling Resource, some authors do an ad hoc correction. [5] replaces the linear Boussinesq approximation with a non-linear
In the present article, no velocity-profile comparisons are included. relationship of the strain rate and rotation rate tensors,
− ρτ ij =
µt* f (Sij ,Wij , k ) (4)
Description of Code and Turbulence Model
The CFD code ISAAC developed by Morrison [3] is used in this where Wij is the rotation rate tensor,
study. ISAAC is a second-order, upwind, finite-volume method 1 ∂ui ∂u j
where advection terms in the mean and turbulence equations are =
Wij ( − ) (5)
2 ∂x j ∂xi
solved by using Roe’s approximate Riemann solver coupled with the
MUSCL scheme. Viscous terms are calculated with a central difference The eddy viscosity is given by,
approximation. Mean and turbulence equations are solved coupled by
using an implicit spatially split, diagonalized approximate factorization
µ * = − ρ kα1
t
(6)
solver. Because of the high order spatial discretization of the MUSCL and α1 is obtained by a solution of a cubic equation. The non-linear
scheme, ISAAC uses flux limiters, to avoid spurious oscillations which relation of Rumsey and Gatski [5] is coupled with the regular k-ω
would otherwise occur in shock waves, discontinuities or sharp changes equations of Wilcox [4].
in the solution domain. The user is able to switch on or off a selected
The Baldwin-Lomax [6] model is a two-layer algebraic zero-
limiter. ISAAC has been developed to test a large range of turbulence
equation model. The eddy viscosity coefficient is given by algebraic
models. Algebraic models, various k-ε and k-ω formulations, and
equations. In the inner layer it follows the Prandlt-Van Driest
Reynolds stress transport models are included in ISAAC.
formulation:
In the present article, the NACA 4412 flow is calculated with the
(µt )inner ρ (κ Dη )2 Ω
= (7)
k-ω turbulence model of Wilcox [4], the Explicit Algebraic Reynolds
Stress Model (EASM) of Rumsey and Gatski [5] and a modification where κ is the von Karman constant (equal to 0.41), D is the van
of the algebraic turbulence model of Baldwin and Lomax [6] done by Driest damping factor, Ω is the absolute value of the vorticity and η is
the present author, in order to improve its accuracy in separated flows. the distance normal to the wall. The damping factor is, D = 1−exp(−
ηuτ/26νw), where uτ = (|τw| /ρw)1/2, τw is the wall shear stress and νw the
It is known that the functional form of the Navier-Stokes equations
wall kinematic viscosity.
is the same for laminar and turbulent flows. In the latter case, however,
time average has been applied to equations, since in turbulent flows, In the outer layer the following equations are used:
each flow or thermodynamic parameter has a mean value and a
(µt )outer = Ccp (0.0168 ρ Fwakeγ ) (8)
random turbulent fluctuation (for example: u= i ui + ui' ). The time
averaged flow equations are known as Reynolds Averaged Navier-
Stokes equations (RANS). In them, new apparent shear stresses, known F1 = ηmax Fmax (9)
as Reynolds stresses: τ ij = ( ρ ui' u 'j ) , have appeared, which need to be
cwkηmax udif
2

calculated. The calculation of the Reynolds stresses is not easy; transport F2 = (10)
equation for each of them must be solved, increasing the total CPU Fmax
cost. Alternatively, the Boussinesq hypothesis is applied, according to Fwake = min ( F1 , F2 ) (11)
which the apparent turbulent shear stresses are related to the rate of
mean strain, through an apparent scalar turbulent or “eddy” viscosity The quantity Fmax is the maximum value of the moment of vorticity:
coefficient, μt. For the calculation of the turbulent viscosity coefficient,

J Vortex Sci Technol

ISSN: 2090-8369 VST, an open access journal Volume 2 • Issue 2 • 1000116
 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

Page 3 of 5

F(η) = ηΩD (12) No grid refinement study has been performed, but comparisons of the
present results with similar ones included in TMR and run with much
The parameter ηmax is the value of η at which F(η) , equation (12),
finer grids (897×257), indicated that for the k-ω and EASM models
is maximum.
the results were similar. Calculations were performed by using the k-ω
The Klebanoff intermittency factor is: and EASM models, as well as the previously described modification
to the Baldwin-Lomax model by Panaras. The experiments indicate
γ = [1 + 5.5(η / δ )6 ]−1 (13)
that the formed separation vortex stays constantly at the same position
The quantity udif is the difference between maximum and minimum (Figure 2d). Thus, a steady-state calculation procedure is sufficient for
velocity in the velocity profile. The thickness of the boundary layer the simulation of the flow. However, for completeness, the steady-state
is defined by: δ = ηmax/Ckleb. The constants appearing in the previous calculations were continued by time-accurate runs, by employing the
relations are: Ccp=1.6, Cwk=0.25, CKleb=0.3. τ-ts scheme of Rumsey et al. [9]. Furthermore, some runs employed
the τ-ts scheme solely, assuming free-stream conditions for the
The wake function (Fwake) is equal to F1 in attached flows and equal initialization of the flow.
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 Figure 2 the standing separation vortices predicted by the tested
in equation (10) is replaced by a reference value, ηref. The value of ηref turbulence models are compared to the experimental evidence, given
differs from flow to flow and it is based on existing semi-empirical by NASA at the Turbulence Modeling Resource web site. It is observed
relations that define the boundary layer parameters of an equivalent in Figure 2a that the k-ω model predicts a very small separation vortex,
flat plate flow of the same Reynolds number. quite different from the experimental one posted by the TMR page
curator. On the contrary, the EASM model predicts a separation vortex
For the estimation of ηref along a flat-plate flow, the semi-empirical of size comparable to that of the experimental one, although its shape
analysis of Falkner [7] is used, which is valid for a Reynolds number is different and it is located further upstream (Figure 2b). Actually,
between 105 and 1010. According to Falkner [7] the boundary layer the separation vortex has an inclined delta shape. Also, the upstream
growth along a flat plate is, boundary layer, which separates and folds around the vortex, bends
and forms an abrupt U-turn between the vortex and the shear layer
0.1285x
δ= (14) of the lower surface. At this point we mention that the shape of the
Re1x/7 prediction of the k-ε-SST turbulence model shown in TMR has very
If this equation is combined with the relation: δ = ηmax/CKleb, then similar shape and behavior. The separation region predicted by the
algebraic model has a shape similar to that given by the EASM model,
n 0.1285 × CKleb 0.03855 (15)
but within the U-turn of the shear-layer a smaller elongated vortex
=max
=
x Rex1/7 Rex1/7 has appeared (Figure 2c). The small vortex is generated by the shear
layer of the lower airfoil surface and it is counter-rotating. Probably
Defining ηref at a characteristic separation length of each examined the algebraic model predicts less turbulent flow than the EASM model.
configuration, Lsep, the final calculation scheme is, However, we note that the two regions of Reynolds stress concentrations
Lsep 0.03855 shown in Figure 1b give the impression of existence of two vortices.
nref = (16) Also, the distribution of the velocity vectors in Figure 1a supports the
L ReL1/7 sep existence of an elongated vortex in area B, or of a turning shear layer,
where ηref has been non-dimensionalized by the length of the body, L. as that predicted by EASM (Figure 2b). More accurate experiments are
Then equation (10) is replaced by, 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
cwkηref udif
2
model is closer to the experimental evidence than that of the baseline
F2 = (17)
Fmax Baldwin-Lomax and of the k-ω model.
For aerospace configurations we propose the equality: Lsep=L. But if Since flux limiters introduce numerical dissipation and the
there is extensive crossflow separation, as in the case of slender bodies presently examined flow is incompressible, the calculations shown in
at incidence, more accurate results are obtained if the characteristic Figure 2 were performed with the flux limiter off. To test their effect,
length is equal to a crossflow length (Lsep=d, for axisymmetric bodies). the calculations were repeated with limiter on. The results are shown
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 [8].

Results
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 Figure 2: Visualization of the separation vortex; (a) k-ω model, (b) EASM, (c)
airfoil and 81 points in the normal to the surface direction (257×81). Algebraic modified, (d) Experimental evidence (taken from TMR).

J Vortex Sci Technol

ISSN: 2090-8369 VST, an open access journal Volume 2 • Issue 2 • 1000116
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

Page 4 of 5

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 [5], the formed non-linear model (EASM) gives improved
results comparable to those posted in TMR by the top-rated SST and
Figure 3: Pressure-coefficient distribution around the surface of the NACA SA models. In addition an improved version of the classical algebraic
4412 airfoil, for the tested models.
zero-equation turbulence model of Baldwin and Lomax [6] 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 [8]. 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
Figure 4: Results of calculations employing a flux limiter in ISAAC code.
the simulated one.
Before closing we wish to add a comment on the convergence of the
in Figure 4. It is seen in this figure that the flux limiter introduces a numerical calculations. Since the separation vortex does not convect
significant dissipation effect. The k-ω model does not predict separation downstream but stays at the same position, the flow is steady. Thus,
at all. The vortex of the EASM has shrunk, and the second vortex of a steady-state calculation procedure is sufficient for the numerical
the algebraic model has disappeared. A comparison of Figures 2 and 3 simulation of the flow. However, because of the effect of the involved
leads to the conclusion that indeed flux limiters introduce considerable turbulence model, it is not known a priori whether a particular
numerical dissipation. The predicted flows are more turbulent than simulation leads to a steady flow. Hence, the involvement of a time-
those given by the same calculation code, but with limiters off. accurate calculation procedure is more appropriate. Furthermore, in
the course of this study we discovered that a time-accurate scheme
In TMR results for the k- ω and EASM models are posted. The
leads smoothly and faster to the converged solution. An example for
separation vortex predicted by the k-ω model has shape and size similar
the calculations that involved the EASM is given in Figure 5. It is
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 time-
accurately, 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”.

Discussion and Conclusions

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 [2], 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 Figure 5: Convergence rate for the EASM calculations, applying steady-state
fine meshes, based on the top-rated turbulence models and posted or time-accurate run procedure.

J Vortex Sci Technol

ISSN: 2090-8369 VST, an open access journal Volume 2 • Issue 2 • 1000116
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

Page 5 of 5

seen in this figure that the time-accurate calculation converges fast 3. Morrison JH (1992) A Compressible Navier-Stokes Solver with Two-Equation
and Reynolds Stress Turbulence Closure Models.
and smoothly. Of course, each iteration of the time-accurate scheme
includes 15 sub-iterations. Still the convergence is faster compared to 4. Wilcox DC (1988) Reassessment of the Scale Determining Equation for
the steady-state runs. The final solutions of the two calculation schemes Advanced Turbulence Models. AIAA Journal 26: 1299-1310.
were found to be identical. 5. Rumsey CL, Gatski TB (2000) Recent Turbulence Model Advances Applied to
Multielement Airfoil Computations. Journal of Aircrafter 38: 904-910.
Acknowledgement
6. Baldwin BS, Lomax H (1978) thin layer approximation and algebraic model for
The author wishes to express his sincere thanks to Professor Spyros Voutsinas separated turbulent flows. AIAA Paper 78-257.
of NTUA for his permission to use the computational facilities of Fluid Section for
computing the flows shown in this article. 7. Falkner VM (1943) A New Law for Calculating Drag: The Resistance of a
Smooth Flat Plate with Turbulent Boundary Layer. Aircraft Engineering and
References Aerospace Technology 15: 65-69.
1. Jameson A (2001) Perspective on Computational Algorithms for Aerodynamic 8. Panaras AG (2015) Turbulence modeling of flows with extensive crossflow
Analysis and Design. Progress in Aerospace Sciences 37: 197-243. separation. Aerospace, to appear.
2. Coles D, Wadcock AJ (1979) Flying-Hot-Wire Study of Flow Past an NACA 9. Rumsey CL, Sanetric MD, Biedron RT, Melson ND, Parlette EB (1996)
4412 Airfoil at Maximum Lift. AIAA Journal 17: 321-329. Efficiency and accuracy of time-accurate turbulent Navier-Stokes computations.
Computers and Fluids 25: 217-236.

Submit your next manuscript and get advantages of OMICS

Group submissions
Unique features:

• User friendly/feasible website-translation of your paper to 50 world’s leading languages

• Audio Version of published paper
• Digital articles to share and explore
Special features:

• 400 Open Access Journals

• 30,000 editorial team
• 21 days rapid review process
• Quality and quick editorial, review and publication processing
• Indexing at PubMed (partial), Scopus, EBSCO, Index Copernicus and Google Scholar etc
• Sharing Option: Social Networking Enabled
• Authors, Reviewers and Editors rewarded with online Scientific Credits
Citation: Panaras AG (2015) Numerical Simulation of Airfoil Flow at High • Better discount for your subsequent articles
Angle of Attack. J Vortex Sci Technol 2: 116. doi:10.4172/2090-8369.1000116 Submit your manuscript at: [Link]/editorialtracking/vortex-science/[Link]

J Vortex Sci Technol

ISSN: 2090-8369 VST, an open access journal Volume 2 • Issue 2 • 1000116

View publication stats