RANS Simulations of Aerodynamic Performance of NACA 0015 Flapped Airfoil

27 Pages·2017·10.5 MB·English

fluids
Article
RANS Simulations of Aerodynamic Performance of
NACA 0015 Flapped Airfoil
SohaibObeid1,*,RatneshwarJha2,*andGoodarzAhmadi1,*
1 Mechanical&AeronauticalEngineeringDepartment,ClarksonUniversity,Potsdam,NY13699,USA
2 AerospaceEngineering,MississippiStateUniversity,Starkville,MS39762,USA
* Correspondence:[email protected](S.O.);[email protected](R.J.);[email protected](G.A.);
Tel.:[phone](S.O.);[phone](R.J.);[phone](G.A.)
AcademicEditor:AsteriosPantokratoras
Received:8August2016;Accepted:23December2016;Published:5January2017
Abstract: Ananalysisof2DsubsonicflowoveranNACA0015airfoilwitha30%trailingedgeflap
ataconstantReynoldsnumberof106 forvariousincidenceanglesandarangeofflapdeflections
ispresented. Thesteady-stategoverningequationsofcontinuityandmomentumconservationare
solvedcombinedwiththerealizablek-εturbulencemodelusingtheANSYS-Fluentcode(Version13.7,
ANSYS,Inc.,Canonsburg,PA,USA).Theprimaryobjectiveofthestudyistoprovideacomprehensive
understandingofflowcharacteristicsaroundtheNACA0015airfoilasafunctionoftheangleofattack
andflapdeflectionatRe=106usingtherealizablek-εturbulencemodel. Theresultsarevalidated
throughcomparisonofthepredictionswiththefreefieldexperimentalmeasurements. Consistent
with the experimental observations, the numerical results show that increased flap deflections
increasethemaximumliftcoefficient,movethezero-liftangleofattack(AoA)toamorenegative
value, decreasethestallAoA,whiletheslopeoftheliftcurveremainsunchangedandthecurve
justshiftsupwards. Inaddition,thenumericalsimulationsprovidelimitsforliftincrement∆C and
l
C valuestobe1.1and2.2,respectively,obtainedataflapdeflectionof50◦. Thisinvestigation
l,max
demonstratesthattherealizablek-εturbulencemodeliscapableofpredictingflowfeaturesoveran
airfoilwithandwithoutflapdeflectionswithreasonableaccuracy.
Keywords: aerodynamics;lift;drag;NACA0015;flap;k-ε
1. Introduction
The aerodynamic shape of lifting objects is one of their most important design parameters.
Thisparameteraffectstheamountoftheliftandthedragforcesthattheairfoilcangenerate. Aircraft
wingsaremainlyintendedtoprovidethemaximumvalueofthelift-to-dragratioL/D.TheL/Dratio
isdeemedasameasureoftheefficiencyofanaircraft[1]andcanbestatedastheamountofpower
(thrust) that is required to impel an aircraft of a certain weight. The aircraft wing performance is
seriouslyimpairedifflowseparationoccurs[2]. Commonly,flowseparationresultsinalossoflift,an
increaseofdrag,diminishedpressurerecovery,etc. Therefore,aconsiderableamountofresearcheffort
hasbeendevotedtothecontrolofflowseparation,andmanymethodstoattainseparationcontrol
havebeensuggested. Amongthemostpromisingactiveflowcontrolapproachesaresuction,blowing,
syntheticjetactuationandwallmovement.
Properflowcontrolmethodshavetheabilitytoreduceskinfrictionandformdrag,increaselift
andimproveflightcontrollability. Thepotentialbenefitsofflowcontrolincludeenhancedrangeand
payload,shortrunwaylandingcapability,improvedecologicalcomplianceandsavingsinoverallfuel
consumption. Forinstance,maintaininglaminarflowovertheentirewingsurfacecanreducetotal
aircraftdragbyasmuchas15%[3].
Fluids2017,2,2;doi:10.3390/fluids2010002 www.mdpi.com/journal/fluids
Fluids2017,2,2 2of27
Recently, there has been considerable interest in flow control, especially in the field of
aerodynamicswiththeintentofincreasingliftanddecreasingthedragofairfoils. Toallowlanding
andtake-offfromshortrunwaysatreducedgroundspeeds,somemodernairplanesareequippedwith
multi-elementhigh-liftdevicesthatgeneratetherequiredhighlift. Slatandsingleormultipleflapsare
typicalexamplesofsuchdevices[4]. Multi-elementwingdesigns,however,arefoundunfavorable
fromaweightandcomplexitypointofview. Thatisthereasonforreplacingthemulti-elementflap
withasingle-hingedflapinthecurrentdesignstoreducethecomplexitywhileincreasingtheefficiency
ofthewing. Whilereducingthecomplexityofthewing,thesingle-hingedflapincreasesthechancefor
flowseparationontheflapatlargedeflectionangles. Topreventoratleasttominimizetheeffects
offlowseparation,theairflowingoverthewingnearthesurfacemustbeenergizedsothatitcould
overcometheeffectsoftheadversepressuregradientencounteredalongtheflap.
Inrecentyears,thepossibleapplicationoftrailingedgeflapsasapotentialvibrationandnoise
controltoolhasgainedconsiderableattentionforapplicationtothehelicopterbladesandwindturbines.
Applicationsoftrailingedgeflapsystemsinhelicopterrotorshavebeenthesubjectofmuchresearch,
bothinthecontextofindividualbladecontrol[5],aswellassmartrotordevelopment[6–8],forthe
purposesofvibrationreduction[9].
The employment of plain flaps in wind turbines has been also a topic of interest to many
researchers. The efficacy of plain flaps in wind turbines is considered favorable both in relation
toloadalleviation,aswellaspowergeneration[10]. Furthermore,smallflapdeflectionscoulddelay
thelaminar-turbulenttransition[11].
Awiderangeofengineeringapplicationsutilizeairfoilsoperatingatrelativelylowandmedium
chord Reynolds numbers ranging from 104 to 106. As noted before, helicopter rotor blades,
small-to-medium-scalewindturbinesandunmannedaerialvehiclesaretypicalexamples[12]. Inthis
Reynoldsnumberrange,theboundarylayerontheuppersurfaceofanairfoilispronetoexperience
flowseparation,evenatalowincidenceangle. Whilethereisanumberofpublishedworksonairfoil
performance at low Reynolds numbers, still there are not many studies, either experimentally or
numerically,thatcovertherangeofRe=106. Inaddition,despitenumerouspublicationsonthelift
anddragofNACAairfoils,betterunderstandingofairfoilswithonehingedflapisstillofinterest[13].
SymmetricNACAairfoilswithathicknessrangefrom9%to18%alsohavemanyapplications
inindustryandfordemonstrationpurposes[14]. Gault[15]classifiedthestallingcharacteristicsof
thelowspeedNACAairfoilsectionsinto: thinairfoilstall,leadingedgestall,trailingedgestallanda
combinationofleadingedgeandtrailingedgestall. Basedonthisclassification,theNACA0015airfoil
wasselectedforthisstudyasitisclassifiedasamediumthicknessairfoil,whichissusceptibleonlyto
trailingedgestallatRe=106.
NumericalstudiesoftheaerodynamicperformanceoftheNACA0015airfoilwithanintegral-type
trailing edge flap were reported by Hassan [16] using the 2D ARC2D Navier–Stokes flow solver
developedattheNASAAmesResearchCenter. Zhangetal.[17]studiedtheaerodynamicperformance
oftheNACA0015airfoilatdifferentspeedsusingtheFLUENTcode. Theyreportedtheflowpatterns,
pressuredistributions,velocityvectorfieldsandturbulenceintensitiesaroundtheairfoil. Numerical
studies for plain (un-flapped) airfoils were presented by Srinivasan et al. [18]. They investigated
variousturbulencemodelsforunsteadyflowsaroundaNACA0015oscillatingairfoil.
The main goal of the numerical simulations presented in this paper is to provide a thorough
understandingofflowfeaturesaroundtheNACA0015airfoilasafunctionofbothangleofattack
andflapdeflectionatRe=106byusingtherealizablek-εturbulencemodel. Theotherobjectiveisto
validatethecomputationalmodelbycomparisonoftheresultswiththeexperimentaldataandearlier
numerical simulation results. The flow problem here is of a boundary layer nature; therefore, the
fluidmotionsneartheairfoilsurfacesareofinterest. TheReynoldsAveragedNavier–Stokes(RANS)
modelshavebeenusedextensivelyforwall-boundedflows[19]. Ontheotherhand,manyresearchers
considertherealizablek-εmodelasuitablechoiceforanalyzingtheboundarylayerflowsunderstrong
adversepressuregradientsorwithseparation[18].
Fluids2017,2,2 3of27
Therealizablek-εmodelsharesthesameturbulentkineticenergyequationasthestandardk-ε
model,buthasanimprovedequationfortheenergydissipationterm,ε,andusesavariablecoefficient
Cµinsteadoftheconstantvaluethatisusedinthestandardk-εmodel.
The realizable k-ε model has been extensively validated for a wide range of flows [20,21].
Therealizablek-εmodelisfoundtoprovidemoreaccuratepredictionscomparedtothestandardk-ε
modelforthecasesofrotatinghomogeneousshearflows,freeshearflowsinjetsandmixinglayers,
channelandboundarylayerflowsandseparatedflows. Inaddition,therealizablek-εmodelproperly
predictsthespreadingrateforaxisymmetricjets,aswellasthatforplanarjets[18].
Inthepresentstudy, the2DsubsonicflowsoveraNACA0015flappedairfoilwereanalyzed
usingtherealizablek-εturbulencemodeloftheANSYS-Fluent(Version13.7)computationalcode.
Thestaticpressuredistribution,theliftandtheskinfrictioncoefficient,aswellasthedragforceover
theflappedairfoilswerestudied. Attentionwasalsogiventothedetailsoftheairflowturbulence
intensitydistributionaroundtheairfoilwithandwithoutflapdeflection.
AnanalysisoftheflowovertheNACA0015airfoilwithzeroflapdeflectionisfirstpresented.
Thedetailsofthemathematicalmodel,meshingschemesutilizedandthecomputationalanalysesare
described. Thisisfollowedbysimulationsofflappedairfoils. Validationsofthecomputationalresults
forcaseswithandwithoutflapdeflectionarealsopresented. Itisconcludedthattherealizablek-ε
turbulencemodeliscapableofcapturingtheflowconditionsovertheairfoilswithandwithoutflaps
withreasonableaccuracy.
2. ComputationalModeling
2.1. Formulation
Inthisstudy,theANSYS-Fluent(Version13.7)codeisusedasthesolveralongwiththerealizable
k-ε turbulence model for solving the Reynolds Average Navier–Stokes (RANS) equations for flow
aroundtheNACA0015airfoilwithandwithoutahingedflap. Thecodesolvestheequationsofthe
conservationofmassandbalanceofmomentum. Additionally,thetransportequationsarealsosolved
forturbulenceproperties. Therealizablek-εmodelhasthefollowingtransportequations[20]
(cid:34) (cid:35)
∂ (ρk)+ ∂ (cid:0)ρku (cid:1) = ∂ (µ+ µt) ∂k +P +P −ρε−Y +S (1)
∂t ∂x j ∂x σ ∂x k b M k
j j k j
(cid:34) (cid:35)
∂∂t(ρε)+ ∂∂xj(cid:0)ρεuj(cid:1) = ∂∂xj (µ+ σµεt)∂∂xεj +ρC1Sε−ρC2k+ε√2 γε +C1εεkC3εPb+Sε (2)
(cid:104) (cid:105) (cid:113)
whereC =max 0.43+ η , η= Sk, S = 2S S .
1 η+5 ε ij ij
Intheseequations, P representsthegenerationofturbulencekineticenergyduetothemean
k
velocitygradient. P isthegenerationofturblencekineticenergyduetobuoyancy. Theturbulent
b
(cid:113)
eddy viscosity is given as µt = ρCµεk2, where Cµ = Ao+ 1AskUε∗, U∗ ≡ SijSij+ΩˇijΩˇij, Ωˇij =
Ω −2ε ω , Ω = Ω −ε ω . Here, Ω is the mean rate of the rotation tensor viewed in a
ij ijk k ij ij ijk k ij
rotatingreferenceframewiththeangularvelocityω .Themodelconstants A and A aregivenas:
Ao =4.04,As = √6cosφ,whereφ= 13cos−1(cid:16)√6W(cid:17),kW = SijSS(cid:101)j3kSki, (cid:101)S = (cid:113)SijSoij, Sij =s 12(cid:104)∂∂uxij + ∂∂uxji(cid:105);
andthemodelconstantsare: C1ε =1.44,C2 =1.9,σk =1.0andσε =1.2.
2.2. EvaluationofC ,C andC
p l f
Thepressurecoefficient,C ,andtheliftcoefficient,C,areevaluatedas,
p l
p−p∞
C = (3)
p
q∞
Fluids 2017, 2, 2 4 of 27
1 𝑥
Fluids2017,2,2 4of27
𝐶 ≈ ∫ [ [𝐶 ] −[𝐶 ] ] 𝑑( ) (4)
𝑙 𝑝 Pressure Side 𝑝 Suction Side 𝑐
0
In these equations, p is the s(cid:90)u1rf(cid:2)ace static pressure; 𝑝 is the free(cid:3) str(cid:16)exam(cid:17) static pressure; 𝑞 is the
C ≈ [C ] −[C∞] d ∞ (4)
free stream dynamic preslsure,0 whicphP reisss urdeeSfiidneed asp S𝑞ucti=on1Si𝜌d𝑣e2 ; acnd x/c is the normalized
∞ 2 ∞
chordI-nwtihsees epoeqsiutiaotnio. nTsh,pe issktihne fsruicrtfiaocne sctoaetifcficpireensts uCrfe ;isp ∞thies trhaetiofr eoef stthreea smursftaacteic sphreeassr usrter;eqs∞s, iτs, thtoe
ffrreeeesstrtreeaamm ddyynnaammicic ppreressssuurere. ,Twhahti cish, i sdefinedasq∞ = 12ρv2∞;andx/cisthenormalizedchord-wise
position. TheskinfrictioncoefficientC istheratioofthesurfaceshearstress,τ,tofreestreamdynamic
f τ
pressure. Thatis, 𝐶𝑓 =𝑞 (5)
τ∞
C = (5)
f q∞
2.3. Computational Domain and Boundary Conditions
2.3. ComputationalDomainandBoundaryConditions
To facilitate the grid generation process, as well as the analysis, similar to the earlier study of
Tofacilitatethegridgenerationprocess,aswellastheanalysis,similartotheearlierstudyof
Hassan [16], the airfoil with the deflected trailing edge flap is treated as a single-element airfoil with
nHoa sgsaapn [b1e6t]w,teheena itrhfoe ilflwapit’hs tlheeaddienflge ceteddgetr aainlidn gthede gbeaflsea poifs tthreea tfeodrwasaards inpgolret-ieolne moef ntthaei rafoiriflowili.t hFnoor
gapbetweentheflap’sleadingedgeandthebaseoftheforwardportionoftheairfoil.Fordeflectingthe
deflecting the flap, solid body rotations were assumed, and a four-point spline smoothing was made
flap,solidbodyrotationswereassumed,andafour-pointsplinesmoothingwasmadefortheresulting
for the resulting airfoil at the chord-wise position corresponding to the location of the flap hinge
airfoilatthechord-wisepositioncorrespondingtothelocationoftheflaphingepoint. Figure1shows
point. Figure 1 shows the schematic geometry of the NACA 0015 with various trailing edge
theschematicgeometryoftheNACA0015withvarioustrailingedgeflapdeflections.
flap deflections.

Figure 1. Geometry of the NACA 0015 airfoil with a 30% trailing edge deflected flap.
Figure1.GeometryoftheNACA0015airfoilwitha30%trailingedgedeflectedflap.
To develop a high efficiency simulation, a C-structured grid topology of a semi-elliptical shape
Todevelopahighefficiencysimulation,aC-structuredgridtopologyofasemi-ellipticalshape
with semi-major diameter 43.5 c and semi-minor diameter 10 c, where c is the chord length, was
with semi-major diameter 43.5 c and semi-minor diameter 10 c, where c is the chord length, was
generated using the preprocessor Gambit (Version 2.4.6, ANSYS, Inc., Canonsburg, PA, USA). The
generated using the preprocessor Gambit (Version 2.4.6, ANSYS, Inc., Canonsburg, PA, USA).
elliptic shape of the grid provides major advantages over other conventional shapes, such as
The elliptic shape of the grid provides major advantages over other conventional shapes, such as
rectangular and semicircle upstream and rectangular downstream. The main advantage is providing
rectangularandsemicircleupstreamandrectangulardownstream. Themainadvantageisproviding
flexibility in using the same mesh for different angles of attack by only shifting the tail part of the
flexibilityinusingthesamemeshfordifferentanglesofattackbyonlyshiftingthetailpartofthe
mesh at the airfoil trailing edge in accordance with the specific angle of attack. For this purpose, it is
meshattheairfoiltrailingedgeinaccordancewiththespecificangleofattack. Forthispurpose,itis
made sure that the edge of the elliptic domain has sufficient vertical length to contain the incoming
madesurethattheedgeoftheellipticdomainhassufficientverticallengthtocontaintheincoming
flow for angles of attack up to 20°. The approach also leads to considerable reduction in the needed
flowforanglesofattackupto20◦. Theapproachalsoleadstoconsiderablereductionintheneeded
number of cells in the mesh in the far-field; thus allowing for the majority of the cells in the mesh to
numberofcellsinthemeshinthefar-field;thusallowingforthemajorityofthecellsinthemeshtobe
be concentrated around the airfoil [22].
concentratedaroundtheairfoil[22].
The distance of the rear pressure far field from the airfoil is 30 c, and the majority of the cells are
Thedistanceoftherearpressurefarfieldfromtheairfoilis30c, andthemajorityofthecells
clustered toward the airfoil surface. This sufficiently large domain was chosen to cover flow
are clustered toward the airfoil surface. This sufficiently large domain was chosen to cover flow
disturbances created by the airfoil and avoid unphysical reflections from the outer edges of the grid.
disturbancescreatedbytheairfoilandavoidunphysicalreflectionsfromtheouteredgesofthegrid.
The schematic of the computational domain and mesh around the airfoil is shown in Figure 2.
TheschematicofthecomputationaldomainandmesharoundtheairfoilisshowninFigure2.
Fluids2017,2,2 5of27
Fluids 2017, 2, 2 5 of 27
Fluids 2017, 2, 2 5 of 27

Figure 2. Domain of calculations and boundary conditions.
Figure2.Domainofcalculationsandboundaryconditions.
The grids construFicgteudre f 2o.r D tohmis asintu odfy c ahlcauvlea taiobnosu at n1d0 4b,o0u0n0d caerlyls c ownidthit ioan fso. ur-node quadrilateral
Theelegmriednst. cToon ssitmruuclatteed thfoe rwtahkies asrteuad ycorhraecvtley,a bDooullet 1[2024] ,0re0c0omcemllesndwedit husainfgo ua rf-innoe dgeridq uwaidthr ilateral
elementq.uTaThdeor ilgasrtiiemdrasul clcaoetlnless tirtnhu tcehteewsde a fakorree atsah rrisae tashteucrdo thyrra ehnca otvltehy ,earDb toyoupletl eo1f0[ 4c2e,20ll]0s0.r Techceeollrsme fwomrieteh, nr edaf eifndoeudru -qsnuionaddgeri laqautefiarndalre icleagltlresir dal with
elemweernet .p lTaoce dsi omnu tloapte o ft hthee wboaukned aarrye ala yceorr rgercidtl yo,n Dthoel luep p[2e2r ]s irdeec aonmdm loewndere sdi due soifn tgh ea a ifrifnoei l gouritdli nwe.i th
quadrilateral cells in these areas rather than other type of cells. Therefore, refined quadrilateral
quaFdigruilraet e3r aslh coewllss itnh eth deestea ailrse aosf rtahteh estrr uthcatunr eodth mere tsyhp, ew ohfi cchel liss. hTihgehrleyf ocrlue,s treerfeinde tdo wquaardd rtihlaet eariarflo ciel lls
cellswereplacedontopoftheboundarylayergridontheuppersideandlowersideoftheairfoil
surface, and it is fine in the vicinity of wake regions and deflected flaps.
were placed on top of the boundary layer grid on the upper side and lower side of the airfoil outline.
outline. Figure3showsthedetailsofthestructuredmesh,whichishighlyclusteredtowardtheairfoil
Figure 3 shows the details of the structured mesh, which is highly clustered toward the airfoil
surface,anditisfineinthevicinityofwakeregionsanddeflectedflaps.
surface, and it is fine in the vicinity of wake regions and deflected flaps.

(a)

(a)

(b)

(b)

(c)
Figure 3. Mesh clustering toward airfoil surface and deflected flap regions. (a) Airfoil with 0° flap
deflection; (b) airfoil with 10° flap deflection; (c) airfoil with 40° flap deflection.

(c)
The boundary conditions for this simulation were the pressure far-field at the computational
FigduorFemig3au.inrMe p3e.es rMhipechslhue rscytle urasitnnedrgi ntsgot awttoiawornadarrday ia rifwrofaoillills sucurorfnfaadccieet iaoannnd d addt eeftfllheeecct etadeid rffloaflipal prseurgreifoganicose.n. (saT.)h (Aea i)rpfAorieilrs fwsouiitrlhew 0fi°at hrfl-af0ipe◦ ldfl ap
boduenfldeactriyo nis; (ab )n aoinrf-orielf wleictthin 1g0 °b folaupn ddeafrlyec ctioonnd; i(tci)o ani rbfoaisl ewdi tohn 4 0“°R fileamp adnenfl eicntvioanr.i ants” used to model a
deflection;(b)airfoilwith10◦flapdeflection;(c)airfoilwith40◦flapdeflection.
free-stream condition at infinity, with the free-stream Mach number and static conditions being
spTechiefi ebdo.u Antd athrey wcoanlld, itthioe nsst afnodr atrhdi sw saimll ufulantcitoionn w beoruen tdhaer yp rceosnsduirteio fna rw-faise luds aedt .t hTeh ec ocmalpcuultaattioionn al
Tdhomepraboiconeu dpnuedrreaispr hyate crtyoh nea dpnirdtei sossnutasrteif oofnarar-rtfyhie ilswd sabilmol ucunoldnaadtriiioetisno, nwa sa ewtr eethltlh ea esa pitrhrfoee islss huseuraerrf-asftacrere.-s fisT echladel cuapltraettsihoseunsrc eoa tmf awrp-auflilet aldti onal
domabionupnedrairpyh iesr ya nanond-rsetfalteicotninagr ybowuanldlacroyn cdointidointioant tbhaeseadir ofoni l“sRuiermfaacnen. Tinhveapriraenstssu” ruesfeadr -tfioe mldobdoelu an dary
isanforene--rsetflreeacmti ncgonbdoituionnd aart yincfoinnidtyi,t iwonithb athseed froene-s“tRreiaemm aMnanchin nvuamriabnert sa”nuds setdattico cmonoddietiloansf rbeeei-nsgtr eam
specified. At the wall, the standard wall function boundary condition was used. The calculation
conditionatinfinity,withthefree-streamMachnumberandstaticconditionsbeingspecified. Atthe
procedures at the pressure far-field boundaries, as well as the shear-stress calculations at wall
wall, the standard wall function boundary condition was used. The calculation procedures at the
pressurefar-fieldboundaries,aswellastheshear-stresscalculationsatwallboundariesweredescribed
intheFluentlecturenotes[23].Insomecases,themeshisadaptedbasedonthestaticpressuregradient,
Fluids2017,2,2 6of27
usingthedefaultmeshadaptationcontrolsettings,sothesolverperiodicallyrefinesthemeshinthe
regionsofhighpressuregradients.
2.4. SettingupoftheNumericalSimulationParameters
Time independent pressure-based solver is used within ANSYS-Fluent (Version 13.7) for the
analysis. Therealizablek-εturbulencemodelisselectedforanalyzingtheboundarylayerflowover
the airfoil. The airflow is assumed to be incompressible. A simple scheme with the Green-Gauss
cell-based gradient implicit formulation of pressure velocity coupling is utilized [23]. For spatial
discretization, the second order upwind differencing scheme which offers several advantages
over a central-differencing formulation for computing viscous flows is used in this work [24].
The ANSYS-Fluent code solves the coupled governing equations of fluid motion simultaneously
andprovidesupdatingcorrectionforthepressurevalueineachiteration[25]. Aconvergencecriterion
of1×10−8wasusedforthecontinuity,x-velocity,y-velocity,kandε. Allsolutionsconvergedwith
thestandardinterpolationschemeforcalculatingcell-facepressureandsecondorderup-winddensity,
momentum,turbulentkineticenergy,turbulentdissipationrateandenergyinterpolationschemesfor
turbulentflow.
3. ResultsandDiscussion
3.1. MeshIndependenceTests
Toensurethatthesimulationresultsareindependentofgridsize,differentcomputationalmeshes
wereinspected. Thisisdonebyrunningcaseswithincreasingnumberofgridcellsuntilthesimulation
results did not change with the use of progressively finer grids. Table 1 lists the properties of six
differentgridswithvaryingdensitythathavebeeninspectedfortheflowpatternaroundtheNACA
0015 airfoil with zero flap deflection. This table lists some details of the grids, including the total
numberofcells,thenumberofnodesN andN placed,respectively,inthexandydirectionsandthe
x y
maximum,minimumandaveragevaluesofnon-dimensionalnormaldistancefromthewall,y+,for
eachgrid. Itisobservedthatallofthegridsinspectedhaveconsiderablylowy+ values,particularly
GridsIII–VI,tosufficientlyresolvetheviscoussub-layer.
Table1.Detailsofgridsusedinmeshsensitivitytesting.
Grid No.ofCells Nx Ny Maxy+ Miny+ Avery+
I 24,910 690 35 32.5 4.8 13.85
II 53,040 884 50 16.5 3.4 06.55
III 76,128 976 65 12 1 05.50
IV 103,192 1184 75 9.2 0.8 04.20
V 141,168 1384 85 8.8 0.7 04.05
VI 367,235 1850 90 1.01 0.01 0.500
Table2providesmoredetailsonthedistributionofcellsfordifferentinspectedmeshesalongthe
xandydirections. Tohaveabettercontrolofthedistributionofthemesh,theairfoilbodywasdivided
intotwoparts;thefrontportionoftheairfoil,whichextendsfromtheleadingedgeupto30%ofthe
chordlength,andtherearportionoftheairfoil,whichincludestherestoftheairfoilbody. Thewake
regionextendsfromtheairfoiltrailingedgeuptotherightfar-fieldend. Ineachcase,thefirst30%of
theairfoilwasmeshedwithnon-uniformgridspacingalongthex-axis,whilethelast70%oftheairfoil
meshedwithuniformgridspacing. Thenumberofgridcellslocatedintheairfoilrearportionismore
thanthree-timesthatlocatedatthefrontpart.
Fluids2017,2,2 7of27
Table2.Detailsofthegridcellsdistributionalongthex-axisandy-axis.
Grid LeadingEdgeupto30%ofChord TrailingEdgeupto70%ofChord WakeRegion Height
I 90 325 275 35
II 110 374 400 50
III 112 384 480 65
IV 116 468 600 75
V 120 544 720 85
VI 128 597 1125 90
Themeshsizeneartheairfoilsurfaceisacriticalparameterforpropersimulationofboundary
layerflowproperties. Thesizeofthefirstcellheightnearthewall∆ywasestimatedbasedonthe
physicalpropertiesofthefluidusedandtheselectedvaluesofthenon-dimensionalnormaldistance
fromthewally+. Table3givesvaluesforestimatedsizeofthisfirstcellheightasafractionofairfoil
chordlength. Themeshgrowthfactorthatcontrolsthetransitionofcellsizeandspecifiestheincrease
inelementsizeforeachsucceedinglayerisfoundtobe1.2.
Table3.Estimatedsizeofthefirstcellheightnearthewallasafractionofairfoilchordlength.
Grid Maxy+ Miny+ Avery+ Max∆y Min∆y Aver∆y
I 32.5 4.8 13.85 7.65×10−4c 1.00×10−4c 3.17×10−4c
II 16.5 3.4 6.55 3.78×10−4c 7.80×10−5c 1.50×10−4c
III 12 1 5.50 2.75×10−4c 2.30×10−5c 1.26×10−4c
IV 9.2 0.8 4.20 2.11×10−4c 1.83×10−5c 9.93×10−5c
V 8.8 0.7 4.05 2.02×10−4c 1.60×10−5c 9.27×10−5c
VI 1.01 0.01 0.500 2.38×10−5c 2.35×10−7c 1.18×10−5c
In this section, the grid inspection studies for flows at Reynolds numbers of 105 and 106 are
performed. Theresultsofthetwo-equationrealizablek-εmodelofShihetal.[20]werecomparedwith
thoseoftheone-equationclosuremodeloftheSpalart–Allmaras(SA).TheSAmodelwasdeveloped
mainlyforsimulatingaerodynamicflows. Thismodelsolvesthemodeledtransportequationforthe
kinematiceddy(turbulent)viscosity[26]. TheresultsobtainedfromSAandR-k-εnumericalmodels
werecomparedwiththeexperimentaldataofRethmel[27].Tables4and5showthepredictedvaluesof
theliftcoefficientsforsixdifferentmeshesbyusingtheSpalart–Allmaras(SA)andtheR-k-εturbulence
modelsforRe=105and106,respectively. Itisapparentthattherealizablek-εsolutionsobtainedby
GridsIV,VandVIareingoodagreementwiththeexperimentaldata. Despitethefactthatthenumber
ofcellsinGridVIisapproximately3.7-timesthatofGridIV,thesegridsprovidednearlythesame
resultsfortheliftcoefficient. ItisalsoseenthattheSAmodelyieldssignificantlyhigherC values
l
comparedtotheexperiments. Regardlessoftheturbulencemodelused,thenumericalresultsforthese
lastthreegridsremainalmostthesamewithincreasingthenumberofcellsfromGridIV–GridVI.
Table4.LiftcoefficientvaluesfortheNACA0015airfoilatα=12◦,δ=0◦andRe=105usingRANS
models.SA,Spalart–Allmaras.
Grid C Exp. C (R-k-εModel) C (SAModel)
l l l
I 0.90±0.005 0.7654 0.9845
II 0.90±0.005 0.8292 1.0276
III 0.90±0.005 0.8412 1.0988
IV 0.90±0.005 0.9091 1.1169
V 0.90±0.005 0.9016 1.1199
VI 0.90±0.005 0.9073 1.1178
Fluids2017,2,2 8of27
Fluids 2017, 2, 2 8 of 27
TTaabbllee 55.. LLiifftt ccooeefffificciieenntt vvaalluueess ffoorr tthhee NNAACCAA 00001155 aaiirrffooiill aatt αα == 1122◦°,, δδ == 00◦° aanndd RRee == 110066 uussiinngg
RRAANNSS mmooddeellss..
Grid C Exp. C (R-k-ε Model) C (SA Model)
Grid lC Exp. l C (R-k-εModel) l C (SAModel)
l l l
I I 1.11.212 ±± 0.00.00055 8.5486. 5×4 610×−11 0−1 1.00341. 0034
II II 1.11.212 ±± 0.00.00055 9.6595.×65 150×−11 0−1 1.16421. 1642
III III 1.11.212 ±± 0.00.00055 9.8991.×89 110×−11 0−1 1.25211. 2521
IV IV 1.11.212 ±± 0.00.00055 1.14613.1 463 1.3161 .316
V 1.12±0.005 1.1491 1.3230
V 1.12 ± 0.005 1.1491 1.3230
VI 1.12±0.005 1.1484 1.3201
VI 1.12 ± 0.005 1.1484 1.3201
Thelaststepinthegridconvergenceinspectionwasfocusedontheanalysisofthedistributionof
The last step in the grid convergence inspection was focused on the analysis of the distribution
thepressurecoefficientalongtheairfoilchord,aswellasthevelocityprofilesontheuppersurface
of the pressure coefficient along the airfoil chord, as well as the velocity profiles on the upper surface
oftheairfoilinsomeselectedsections. Figure4presentsthechordwisedistributionsofthepressure
of the airfoil in some selected sections. Figure 4 presents the chordwise distributions of the pressure
coefficient (C ) profile for the airfoil predicted by Grids I, IV and VI using the R-k-ε turbulence
coefficient (Cpp) profile for the airfoil predicted by Grids I, IV and VI using the R-k-ε turbulence model
model for α = 12◦, δ = 0◦ and Re = 105 (a) and Re = 106 (b). Here, the predicted distributions of
for α = 12°, δ = 0° and Re = 105 (a) and Re = 106 (b). Here, the predicted distributions of the pressure
thepressurecoefficientarecomparedwiththeexperimentaldataofRethmel[27]. Itisseenthatthe
coefficient are compared with the experimental data of Rethmel [27]. It is seen that the model
modelpredictionsforthepressurecoefficientaregenerallycomparablewiththeexperimentaldata.
predictions for the pressure coefficient are generally comparable with the experimental data. In
In particular, the predictions obtained by Grids IV and VI are very close to each other and to the
particular, the predictions obtained by Grids IV and VI are very close to each other and to the
distributionoftheexperimentaldataforbothReynoldsnumbers.
distribution of the experimental data for both Reynolds numbers.

FFiigguurree 44.. CCoommppaarriissoonno offp prersessusurerec oceofefifcfiiecinetnctu cruvrevseasl oanlogntgh ethseid seisdoesf aoifr faoiirlfforiol mfrodmiff edrieffnetrgenritd gsraitdαs a=t
1α2 ◦= a1n2d° aδn=d 0δ◦ =w 0i°t hwtihthe tehxep eerxipmeerinmtaelndtaalt adoaftaR oeft hRmetehlm[2e7l] [a2t7(]a a)tR (ae)= R1e0 =5 1a0n5d a(nbd) (Rbe) =Re1 0=6 1.06.
Figure 5 shows the velocity profiles on the upper surface of airfoil at the 20% cord location
Figure5showsthevelocityprofilesontheuppersurfaceofairfoilatthe20%cordlocationfrom
tfhroemle athdein lgeaeddigneg aesdpgree adsi cptereddbicytetdh ebrye athliez arebalelizka-εblme ok-dεe mluosdinelg uGsirnidgs GI,riIdVs aI,n IdVV aInfdo rVαI f=or1 α2◦ =a 1n2d°
δan=d0 δ◦ =fo 0r°R feor= R1e0 5= a1n0d5 aRned= R1e0 =6 .1T0h6. eTvheer vtiecratlicaaxli sax(yis) i(ny)F iing uFrigeu5rde e5n doetensotthees tvheer tvicearltidcaislt adniscteanfrcoem frtohme
the surface of the airfoil as a fraction of the chord length. It is seen that the model predictions for
surfaceoftheairfoilasafractionofthechordlength. ItisseenthatthemodelpredictionsforGridsIV
Grids IV and VI are almost identical. That is, further refinement of Grid IV does not have a
andVIarealmostidentical. Thatis,furtherrefinementofGridIVdoesnothaveanoticeableeffecton
noticeable effect on the velocity profile.
thevelocityprofile.
Fluids2017,2,2 9of27
Fluids 2017, 2, 2 9 of 27

FigurFeig5u.rVee 5lo. Vcietlyocpitryo fiprleosfiloens othn ethuep uppepresr usurfrafacceeo off aaiirrffooiill aatt aa 2200%% cocrodr dlocloatciaotnio frnomfro tmhe tlheaedlienagd eidngge edge
atα=at1 α2◦ = a1n2d° aδnd= δ0 ◦= a0n° danfdo rfo(ra ()aR) Ree= =1 10055 aanndd ((bb)) RRee = =10160. 6.
The presented grid sensitivity analysis shows that the model predictions with the use of
TGhreidps rIeVs eanntded VgIr iadres evnesryit icvliotsye atnoa leyaschis oshthoewr santhda ttot htehem eoxdpeelripmreendtiacl tidoantas. wDiuthe tthoe tuhsee hoigfhG rids
IVancdomVpIuatraetivoenrayl ccolsots aestsooceiaatcehd owthither thaen dusteo otfh Gereidx pVeIr, iGmreidn tIVal wdiattha a. Dtoutael tnoutmhbeehr iogfh 10c3o,m19p2 ucetalltsi onal
costawssaos cuisaetded inw thiteh suthbeseuqsueenotf aGnarilydseVsI. ,GridIVwithatotalnumberof103,192cellswasusedinthe
subsequentanalyses.
3.2. NACA 0015 Airfoil with Zero Flap Deflection
3.2. NACATh0e0 1a5irfAloiwrfo pilrowpitehrtZieesr aorFoluanpdD thefle eNctAioCnA 0015 with zero flap deflection are first studied, and
the corresponding distributions of static pressure and velocity magnitude at different incidence
TheairflowpropertiesaroundtheNACA0015withzeroflapdeflectionarefirststudied,andthe
angles are evaluated and compared with the published numerical results and/or experimental data.
correspondingdistributionsofstaticpressureandvelocitymagnitudeatdifferentincidenceanglesare
Figure 6 shows the static pressure and velocity magnitude contours around the NACA 0015
evaluatedandcomparedwiththepublishednumericalresultsand/orexperimentaldata.
airfoil for a few selected incidence angles. As NACA 0015 is a symmetric airfoil, at zero incidence
Figure6showsthestaticpressureandvelocitymagnitudecontoursaroundtheNACA0015airfoil
angle, the static pressure and velocity distribution over the airfoil are symmetric, which results in
forafzeewro sleiflet cftoerdcei nacnidd ean csetaagnngalteiosn. ApsoNintA, CexAac0t0ly1 5ati stahes ynmosme eotfr icthaei rafoiriflo,ial.t Tzehreorei nacried ernegcieonasn golfe ,the
staticapcrceelsesruarteeda nfldowvse loovceirt yanddis turnibduert iothne oavirefroitlh tehaati rrfeoaiclha rtehes yhmigmheesttr sicp,ewedh iacth three sauirltfsoiiln mzaexriomluifmtf orce
andathsitcakgnneassti pooninpto. Tinhte, veexlaocctitlyy ias thtihghe (nmoasrekoedf tbhye raedir sfopiolt.sT) ihne trhee alorewr pergeisosnusreo rfegaicocne laenrda tveidcefl voerwsas. over
anduTnhdee mrtahxeimauirmfo pilrtehssautrree oaccchutrhs eath tihgeh setsatgsnpaetieodn aptotinhte wahirefno itlhme vaexliomciutym ist hzeicrok.n esspoint. Thevelocity
At an incidence angle of 5°, the contours of static pressure over the airfoil become asymmetric;
ishigh(markedbyredspots)inthelowpressureregionandviceversa. Themaximumpressureoccurs
the pressure on the upper surface becomes lower than the pressure on the lower surface; regions of
atthestagnationpointwhenthevelocityiszero.
high pressure on the airfoil lower surface become dominant; and a lift coefficient of 0.531 is
Atanincidenceangleof5◦,thecontoursofstaticpressureovertheairfoilbecomeasymmetric;
generated due to the pressure imbalance.
thepressureontheuppersurfacebecomeslowerthanthepressureonthelowersurface;regionsof
Figure 6 shows that as the angle of attack increases, the stagnation point is shifted towards the
highpressureontheairfoillowersurfacebecomedominant;andaliftcoefficientof0.531isgenerated
trailing edge on the bottom surface; hence, it creates a low velocity region at the lower surface of the
duetoairtfhoeilp arneds sau hreigihm vbealolacnitcye r.egion on the upper side of the airfoil. Thus, the pressure on the upper
Fsiidgeu roef 6thseh aoiwrfositl hisa tloawsethr ethaanng tlheeo afmatbtiaecnkt ipnrcersesuarsee,s w,thheeresatsa gthnea tpiorensspuorein otni sthshe ilfotewdert oswidae rids sthe
trailinhgigehderg ethoann tthhee abmotbtioemnt spurerfsasucere;.h Tehnecree,foitrec,r einactreesasainlogw thve eilnocciditeynrceeg ainognlea tist haessloocwiateerds uwriftha cteheo fthe
airfoiilnacrnedasae hofi gthhe vlieflto ccoietfyficrieegnito, naso wnetllh aes uthpep ienrcrseidasee ooff tthhee dariargfo ciol.efTfihciuesn,t.t hTheisp rinecsrseuarsee oinn ththe elifut pper
coefficient continues up to a maximum, after which the lift coefficient decreases.
sideoftheairfoilislowerthantheambientpressure,whereasthepressureonthelowersideishigher
It is also seen that the flow field around the airfoil varies markedly with the incidence angle. In
thantheambientpressure. Therefore,increasingtheincidenceangleisassociatedwiththeincrease
addition to the changes in velocity and pressure distributions, the properties of the boundary layer
oftheliftcoefficient,aswellastheincreaseofthedragcoefficient. Thisincreaseintheliftcoefficient
flow along the airfoil surface also change. At low incidence angles up to about 12°, the boundary
continuesuptoamaximum,afterwhichtheliftcoefficientdecreases.
layer is fully attached to the surface of the airfoil, and the lift coefficient increases with angle of
It is also seen that the flow field around the airfoil varies markedly with the incidence angle.
attack, while the drag is relatively low. With the increase of incidence angle, the boundary layer is
Inaddthiitciokennteod.t hWehcehna nthgee sinicnidveenlcoec iatnygalen dofp trhees sauirrfeoidl iiss tirnibcruetaisoends ,toth aebpourot p1e3r°t ioers loarfgtehre, tbhoeu anddvaerrysel ayer
flowalongtheairfoilsurfacealsochange. Atlowincidenceanglesuptoabout12◦, theboundary
layerisfullyattachedtothesurfaceoftheairfoil,andtheliftcoefficientincreaseswithangleofattack,
whilethedragisrelativelylow. Withtheincreaseofincidenceangle,theboundarylayeristhickened.
Fluids2017,2,2 10of27
Fluids 2017, 2, 2 10 of 27
Whentheincidenceangleoftheairfoilisincreasedtoabout13◦orlarger,theadversepressuregradient
pressure gradient imposed on the boundary layers become so large that separation of the boundary
imposedontheboundarylayersbecomesolargethatseparationoftheboundarylayeroccurs.Aregion
layer occurs. A region of recirculating flow over the entire upper surface of the airfoil forms, and the
ofrecirculatingflowovertheentireuppersurfaceoftheairfoilforms,andtheregionofhigherpressure
region of higher pressure on the lower surface of the airfoil becomes smaller. Consequently, the lift
onthelowersurfaceoftheairfoilbecomessmaller.Consequently,theliftdecreasesmarkedly,andthe
decreases markedly, and the drag increases sharply. This is a typical condition in which the airfoil is
dragincreasessharply. Thisisatypicalconditioninwhichtheairfoilisstalled.
stalled.

FigurFeig6u.rSet a6t. iSctaptrice spsruersesuarne danvde vloecloitcyitym magagnnitiutuddee ccoonnttoouurrss aarorouunnd dthteh NeANCAAC 0A01050 a1i5rfoaiilr faot idliafftedreinffte rent
incidence angles.
incidenceangles.
For a further increase of the airfoil incidence angle to 20°, the stagnation point shifts
Forafurtherincreaseoftheairfoilincidenceangleto20◦,thestagnationpointshiftssignificantly
significantly further towards the trailing edge on the bottom surface. The recirculating flow region
furthebrectoomweas rddosmthineantrt aailnidn gcoevdegrse tohne tehnetirbeo uttpopmer ssuurrffaaccee. oTf hthee reaicrifrociul, laatnidn gthfleo awirfrloewgi oins fbuellcyo mes
dominseapnatraatnedd cfroovme rtsheth uepepnetri rseurufapcpee orf stuhref aaicrefooilf. tThheisa lierafodisl ,toa nfudrtthheer arierdfluocwtioins fouf ltlhyes leifpt aforarctee danfrdo am the
uppersesvuerrfea icnecroefatshee oaf itrhfeo dilr.aTgh fiosrcleea. dstofurtherreductionoftheliftforceandasevereincreaseofthe
dragforceT.he air flowing along the top of the airfoil surface experiences a change in pressure, moving
from the ambient pressure in front of the airfoil, to a lower pressure over the surface of the airfoil,
Theairflowingalongthetopoftheairfoilsurfaceexperiencesachangeinpressure,movingfrom
then back to the ambient pressure behind the airfoil. The region where fluid must flow from low to
theambientpressureinfrontoftheairfoil, toalowerpressureoverthesurfaceoftheairfoil, then
high pressure (adverse pressure gradient) could cause flow separation. If the adverse pressure
backtotheambientpressurebehindtheairfoil. Theregionwherefluidmustflowfromlowtohigh
pressure(adversepressuregradient)couldcauseflowseparation. Iftheadversepressuregradientis
toohigh,thepressureforcesovercomethefluidinertialforces,andtheflowseparatesfromtheairfoil