Effect of trailing edge shape on the separated flow characteristics around an airfoil at low Reynolds

18 Pages·2017·14.24 MB·English

Effect of trailing edge shape on the separated flow characteristics around an airfoil
at low Reynolds number: A numerical study
Nikitas Thomareis and George Papadakis
Citation: Phys. Fluids 29, 014101 (2017); doi: 10.1063/1.4973811
View online: http://dx.doi.org/10.1063/1.4973811
View Table of Contents: http://aip.scitation.org/toc/phf/29/1
Published by the American Institute of Physics
Articles you may be interested in
Compressibility effects on the flow past a rotating cylinder
Phys. Fluids 29, 016101016101 (2017); 10.1063/1.4973564
Vortex identification from local properties of the vorticity field
Phys. Fluids 29, 015101015101 (2017); 10.1063/1.4973243
Effects of grid geometry on non-equilibrium dissipation in grid turbulence
Phys. Fluids 29, 015102015102 (2017); 10.1063/1.4973416
Scalar mixing and strain dynamics methodologies for PIV/LIF measurements of vortex ring flows
Phys. Fluids 29, 013602013602 (2017); 10.1063/1.4973822
PHYSICSOFFLUIDS29,014101(2017)
Effect of trailing edge shape on the separated flow characteristics around
an airfoil at low Reynolds number: A numerical study
NikitasThomareisa) andGeorgePapadakis
DepartmentofAeronautics,ImperialCollegeLondon,LondonSW72AZ,UnitedKingdom
(Received8July2016;accepted26December2016;publishedonline17January2017)
Direct numerical simulations of the flow field around a NACA 0012 airfoil at Reynolds number
50000andangleofattack5◦with3differenttrailingedgeshapes(straight,blunt,andserrated)have
beenperformed.Bothtime-averagedflowcharacteristicsandthemostdominantflowstructuresand
theirfrequenciesareinvestigatedusingthedynamicmodedecompositionmethod.Itisshownthatfor
thestraighttrailingedgeairfoil,thismethodcancapturethefundamentalaswellasthesubharmonicof
theKelvin-Helmholtzinstabilitythatdevelopsnaturallyintheseparatingshearlayer.Thefundamental
frequencymatcheswellwithrelevantdataintheliterature.Theblunttrailingedgeresultsinperiodic
vortexshedding,withfrequencyclosetothesubharmonicofthenaturalshearlayerfrequency.The
shedding,resultingfromaglobalinstability,hasanupstreameffectandforcestheseparatingshear
layer.Duetoforcing,theshearlayerfrequencylocksontothesheddingfrequencywhilethenatural
frequency(anditssubharmonic)issuppressed.Thepresenceofserrationsinthetrailingedgecreatesa
spanwisepressuregradient,whichisresponsibleforthedevelopmentofasecondaryflowpatterninthe
spanwisedirection.Thispatternaffectsthemeanflowinthenearwake.Itcanexplainanunexpected
observation,namely,thatthevelocitydeficitdownstreamofatroughissmallerthanthedeficitafter
a protrusion. Furthermore, the insertion of serrations attenuates the energy of vortex shedding by
de-correlatingthespanwisecoherenceofthevortices.Thisresultsinweakerforcingoftheseparating
shearlayer,andboththesubharmonicsofthenaturalfrequencyandthesheddingfrequencyappear
inthespectra.PublishedbyAIPPublishing.[http://dx.doi.org/10.1063/1.4973811]
I. INTRODUCTION chordratiolargerthan25%.Thickairfoilswithsharptrailing
edgeshoweverhavepooraerodynamiccharacteristics,while
Performance improvement of lifting devices has always
blunttrailingedgesofferstructuralbenefitswithoutsacrificing
beenachallengefortheaerospaceindustry.Geometricmodifi-
theaerodynamicperformance.2,5
cationofthewingsectionisoneofthemostcommonmethods
The exposed blunt part of the airfoil however leads to
usedtoattainthisgoal.Theabsenceofacomplicatedactua-
increased pressure drag due to shedding. For wind turbines,
tionmechanismhasmadepassivecontrolmethodsattractive,
increaseddragisnotsuchanimportantissuefortheinboard
the most well known and widely used being the extendable
regionoftheblade.2Nevertheless,thebenefitsofbluntairfoils
flapsandslatsemployedsincethefirstdaysofaviation.Since
willincreaseifthepenaltyofincreaseddragcanbemitigated.
the1950s,othermethodshavebeenstudied;forexample,the
Biomimeticsinspiredseveralinvestigatorstoimitatetheser-
truncationoftherearpartofthewingtocreateablunttrailing
ratedgeometriesobservedinanimalssuchaswhalesorowls.
edge,1 leading to what is known as a flatback profile. These
Tanner6 used M-shaped serrations at the trailing edge of a
earlyinvestigationsshowedthepotentialimprovementofthe
bluntprofileandfoundadecreaseinthebasepressureofupto
maximumliftcoefficientforthesameflowconditions.More
64%comparedtothestraightblunttrailingedge.Gai7 found
recentinvestigationshavedemonstratedadditionalbenefitsof
similarresultswithTannerandstudiedtheeffectoftriangular
the flatback airfoil, both experimentally and numerically,2–4
serrationswithanglesof60◦and120◦,whichwerealsoshown
suchasincreasedliftcurveslope,optimizedstructuralcharac-
tohavepositiveeffects.Rodriguez8 investigatedtheeffectof
teristics,anddecreasedsensitivitytoleading-edgetransition.
squaredserrationsona2Dbodyandfoundareductionofupto
TherebenefitsareimportantforlowtomediumReynoldsnum-
40%ofthetotaldrag,withanadditionalstudyofthelongitu-
berapplications,suchasunmannedairvehiclesorsmallwind
dinalvorticesinthenear-wakeofthebodyandtheireffecton
turbines.
theaerodynamicsmechanisms.AmorerecentworkbyKrentel
Forsmallwindturbinesinparticular,blunttrailingedge
andNitsche9onatruncatedNACA0012airfoilusingvarious
airfoilsoffersignificantadvantages.2 Thecyclicgravitational
geometricalpatternsshowedadecreaseinthedragby29%for
andaerodynamicloadresultsinlargebendingmomentsinthe
aspecificReynoldsnumber(44000).Thelatest(experimen-
inboardregionoftheblades.Thisnecessitatesthickandstruc-
tal)studyinthisareabyNedic´andVassilicos10hasshownthat
turallyrobustairfoilsectionsclosetothehub,withthicknessto
theeffectofserratedtrailingedgescanbefurtherimprovedby
addingaself-similarrepeating(i.e.,fractal)patterninsmaller
a)Author to whom correspondence should be addressed. Electronic mail: scales.Itwasshownthatthesinglescaletriangularserrations
[e-mail] reducedvortexsheddingandimprovedthelift-to-dragratiofor
1070-6631/2017/29(1)/014101/17/$30.00 29,014101-1 PublishedbyAIPPublishing.
014101-2 N.ThomareisandG.Papadakis Phys.Fluids29,014101(2017)
higheranglesofattack,andthattheadditionofafractalmulti- suitableforunstructuredgridsforpredictingtheincompress-
scalepatternreducedevenmorethesheddingenergywithout ibleflowintransitionalseparatingbubbles.ANACA0012air-
affectingtheimprovedaerodynamiccharacteristics. foilwithastraighttrailingedgewasexamined.DNSstudiesat
A thorough study on the flow physics behind a bluff Re = 50000 were performed at two angles of attack, 5◦ and
bodywithasinusoidaltrailingedgeofconstantthicknesswas 8◦, and used as a reference to assess the performance of the
madebyTombazis.11 Hefoundthatthethree-dimensionality sub-gridscalemodels.Recently,DNSwithReynoldsnumber
effects in the wake can reduce the correlation length of the 400000wasreported.27
basepressure,causingadecreaseinthepressuredrag.Further Simulations around airfoils withmodified trailing edges
investigation of the wake of bodies with spanwise undula- are even more rare. Jones and Sandberg28 performed DNS
tion12,13 showedthatthebasepressureaftertheblunttrailing simulations for a NACA 0012 with serrated trailing edge
edgecanbecorrelatedwiththefrequenciesofthevortexdis- at the same Reynolds and Mach number as their previous
locations arising due to the spanwise inhomogeneity of the investigation for a straight trailing edge.25 A serrated plate
body. The aforementioned cases of bluff bodies with wavy with very small and uniform bluntness (equal to 1.2×10−3
trailing edges were also studied numerically with similar times the airfoil chord) was appended to the straight trailing
conclusions.14 edge.Althoughthemainfocusofthepaperwasthestudyof
Mostoftheaforementionedinvestigationsareexperimen- the acoustic behavior of the airfoil, the analysis of the flow
tal.Asalreadymentioned,airfoilswithmodifiedtrailingedges field close to the trailing edge revealed the breakup of large
aresuitableforlowtomediumReynoldsnumberapplications structuresfromtheboundarylayeraswellasthecreationof
forwhichDirectNumericalSimulation(DNS)isfeasible.An horseshoe vortices from the serrations. The same authors in
important characteristic in the aerodynamic behavior of air- anotherpaper29foundthatthepresenceofserrationsdoesnot
foilsatmoderateReynoldsnumbersisthelaminarseparation, changethehydrodynamicfieldontheairfoilupstreamofthe
whichoccursshortlyaftertheleadingedgeoftheairfoil,giving serrations, including the behavior of the laminar separation
risetoahighlycomplexanddynamicfield.15Thebehaviorof bubble.
theseparatedflowdependsontheangleofattack,theReynolds Inthepresentpaper,weinvestigatetheeffectoftrailing
number,andthetypeoftheairfoil.Kotapatietal.16distinguish edgemodificationstothetimeaverageaswellasthedynamic
3 different scenarios: attached flow at small angles of attack featuresoftheflow.Thecentralaimistostudytheinteraction
(caseA),laminarseparation,transition,andreattachment(case of the separating shear layer with trailing edges of different
B), and massively separated post-stall flow at high angles of shapes.Weareinterestedonlyinthehydrodynamiceffectof
attack(caseC). serrationsandnottheireffectoftheacousticfield.Themodi-
Most of the DNS studies of separating and transitional fiedtrailingedgesareobtainedbyremovingmaterialfromthe
flowshoweverhavebeenperformedonplanargeometries.16–23 mainbodyoftheairfoil,therebyexposingafinitebluntness,
The adverse pressure gradient required for the formation of which is around 30 times larger compared to that of Jones
a laminar separation bubble (encountered in real airfoils) is and Sandberg.28 As will be seen later, this difference affects
reproduced by a suitable boundary condition at the upper significantly the dynamics of the separating shear layer. We
boundaryofthecomputationaldomain.Severalofthesestud- use the Dynamic Mode Decomposition (DMD) to identify
ies are concerned with the control of the laminar separation. the structures generated by the shear layer and the exposed
There are only a few DNS studies on actual airfoils. The bluntness.
latter simulations are more challenging and require curvi- The paper is organized as follows: Section II pro-
linear (or unstructured grids), but are more realistic as the vides details on the airfoils examined and computational
pressuregradientthatresultsfromtheinteractionofthesep- methodused;SectionIIIbrieflydescribestheDMDmethod;
aration bubble with the potential flow around the airfoil and Section IV presents the results for the 3 airfoils exam-
the Kutta condition is part of the solution and not imposed ined; and Section V summarises the main findings of this
externally. work.
Shanetal.24performedDNSaroundaNACA0012at4◦
and Re based on the chord length equal to 105. The flow at
II. COMPUTATIONALSETUPANDVALIDATION
these conditions corresponds to case B scenario. Indeed the
authors captured the laminar separation, transition, and reat- Theflowaroundthreeairfoils,allbasedonNACA0012,
tachmentoftheshearlayer.Thesheddingfromtheseparated wasinvestigatedwithDNS.Theairfoilshaddifferenttrailing
shearlayerwasattributedtotheKelvin-Helmholtzinstability edgeshapes:straight,truncated(blunt),andserrated.Theair-
butnocharacteristicfrequencywasmentioned.Jonesetal.25 foilwiththestraighttrailingedgehadastandardNACA0012
investigatedtheeffectofforcingonthebehaviourofalami- crosssection.Theblunttrailingedgeairfoil(alsoreferredto
nar separation bubble for a NACA 0012 at Re = 50000 and asflatbackairfoilbelow)wasformedbytruncatingtheNACA
Mach 0.4. The forcing used in Ref. 25 was a single fre- 0012sectionatadistanceh =0.133C fromthetrailingedge
quency volume forcing term applied in a predefined region, of the unmodified airfoil, thus exposing a uniform bluntness
withaninducedspanwisevariation.Theyfoundthatforcing ofthickness(cid:15) =0.037C.Theairfoilwiththeserratedtrailing
improvestheaerodynamicperformance,requiringlittleenergy edge is shown in Figure 1(a). The serrations had triangular
input. A mechanism by which turbulence can sustain itself shape,exposingabluntnessthattaperedlinearlyinthespan-
aftertheremovalofforcingwasproposed.Lehmkuhletal.26 wisedirection;itwasmaximumatthetroughandzeroatthe
investigated the capabilities of two sub-grid scale models peak.Themaximumbluntnesswasequaltothatoftheflatback
014101-3 N.ThomareisandG.Papadakis Phys.Fluids29,014101(2017)
employedtogetherwiththeBoomerAMGalgebraicmultigrid
preconditioneroftheHyprepackage.32
Thesizeofthecomputationaldomainin3Dmustbelarge
enough to capture the potential flow as accurately as pos-
sible, while at the same time minimising the computational
expense. The domain of Jones et al.25 extends a radius 5.3C
away from the airfoil and computations in a larger domain
with radius 7.3C showed a small variation in the pressure
distribution. This modest domain size was found to be suf-
ficienttocapturethepotentialflowandthiswasattributedto
FIG.1. Serratedtrailingedgegeometryandcharacteristicdimensions.(a)
thecharacteristics-basedboundaryconditionsemployed.The
Geometryoftheairfoilwithtriangularserrations.(b)Planarviewofthetrailing
edge. domainofShanetal.24extended3Cintheupstreamanddown-
streamdirectionsand4Cinthecross-streamdirection;againa
airfoil.AplanarviewofthetrailingedgeisshowninFigure compressibleflowsolverwasusedwithnon-reflectivebound-
1(b). Three geometric variables define the shape of the ser- ary conditions in all boundaries. Zhang and Samtaney33 and
rations:theheight(h),theperiod(L),andthehalf-angle(φ), Lehmkuhletal.26 employedlargerdomains,witharadiusof
withvaluesh=0.133C,L=0.266C andφ=45◦. 10Cand20C,respectively.Inthepresentwork,forthestraight
For the discretization of the computational domain, a trailingedgeairfoil,3Dsimulationsintwodomainswereper-
C-gridtopologywasemployed.TheX,Y,andZ coordinates formed. In the small domain, the inlet boundary is located
are defined as follows: X is the free stream direction, Y is 6.0C upstreamoftheleadingedge,theoutletboundary4.0C
thecross-stream,andZ isthedirectionalongthespanofthe downstreamofthetrailingedge,andthespanwiseextentwas
domain. The origin of the coordinate system is the leading Lz = 0.1C. In the larger domain, the inlet, top, and bottom
edgeoftheairfoil. boundaries are located 18C away from the airfoil, the outlet
The chord Reynolds number was Re = 50000 and the boundary20Cdownstreamofthetrailingedge,andthespan-
C
angle of attack α=5◦. This combination results in laminar wiselengthis0.2C.Comparisonbetweentheresultsobtained
separationandturbulentreattachmentandenablesthestudyof fromthetwodomainswillbepresentedlater.
theeffectofthemodifiedtrailingedgeontheseparatedflow For the blunt trailing edge simulations, the spanwise
features and the wake. Unless otherwise stated, frequencies extent was kept at Lz=0.2C in order to capture the 3D
arenon-dimensionalisedusingfreestreamvelocityU andthe instabilities of the vortices emanating from the sharp trail-
∞
chordlengthC.Fortheairfoilswiththemodifiedtrailingedges, ing edge. These vortices scale with the edge thickness, (cid:15),
thenominalchordlengthofthestandardNACA0012profile and are expected to develop 3Dinstabilities at the Re exam-
fromwhichtheyoriginatewasusedasareferencelength. ined. Note that the Re based on (cid:15) is 1850, and 3D instabil-
An in-house code, called PantaRhei, was employed for ities develop (at least for the flow around a cylinder) when
thenumericalsimulations.TheincompressibleNavier-Stokes Re is larger than about 190.34 There are two types of insta-
equations bilities (termed modes A and B) that have different span-
wise length scales (4(cid:15) and (cid:15), respectively).34 The spanwise
∂ui + ∂ujui =−∂p + 1 ∂2ui, (1) extent of Lz=0.2C is therefore wide enough to accommo-
∂t ∂xj ∂xi ReC ∂x2 date both instabilities. For the serrated trailing edge air-
j
foil, the spanwise distance was two serration periods, i.e.,
∂u
i =0 (2) Lz=0.532C.
∂x
i Thetopandbottomsurfaceswereassignedfree-slipcon-
aresolvedusingtheFiniteVolumeMethod(FVM)onacol- ditions, while the side boundaries were periodic. The flow
located,unstructuredgrid.Theconvectivetermsareapproxi- velocitywasprescribedattheinlet,whileattheoutletaone-
matedwitha2ndordercentralscheme.Fortimeadvancement, dimensionalconvectiveboundaryconditionthatemploysthe
thesecondorderbackwarddifferencingschemeisused.Vis- localcellvelocitywasused.Thisboundaryconditionallows
cous terms are treated fully implicitly. The convective terms the vortices to exit smoothly the computational domain with
arelinearizedbyextrapolatingu fromthetwoprevioustime minimumreflection.
j
steps and treating u implicitly. Continuity is enforced via Grid resolution was very fine close to the airfoil, with
i
thePressureImplicitwithSplittingOperator(orPISO)algo- 740 nodes around the surface (with appropriate clustering at
rithm.30ThecodeisparallelisedusingthePETSc31library.For the suction side to adequately resolve separation, transition,
thesolution of thepressure equation, the GMRES methodis andreattachment),220 nodes at thewall-normal distanceup
FIG.2. Contour plots of the ratio of
thecharacteristicgridsizetothelocal
Kolmogorovscale.
014101-4 N.ThomareisandG.Papadakis Phys.Fluids29,014101(2017)
FIG.3. Spanwiseautocorrelationfunc-
tionsforthestraightandblunttrailing
edgeairfoils.Thefirstpointisinthesep-
arationbubble,the2ndand3rdpoints
areinthereattachedturbulentboundary
layer,andthelastpointsareinthevery
nearwake.Thecrossstreamcoordinates
ofthepointsareforclaritymentioned
here: for (a) Y/C = 0.038, 0.00037,
(cid:12) (cid:12)
0.031, 0.075andfor(b)Y/C=0.025,
(cid:12) (cid:12) (cid:12)
0.019, 0.041, 0.062goingfromthe
firsttolastpoint.
to0.15C and500nodesatthewakeupto2C downstreamof and the blunt trailing edge airfoil at 4 points (3 inside the
thetrailingedge.Thenumberoflayersalongthespanvaried boundary layer and 1 in the near wake). The correlation
from 100 (for the straight trailing edge) to 140 (for the ser- curvesdependontheselectedpoint.Thebestdecorrelationfor
ratedtrailingedge),leadingtoatotalgridsizebetween50and both airfoils is obtained for the most upstream point, which
70 × 106 cells. The resolution in the wall-normal direction is located before transition. For the other points, the values
is at worse ∆y+<1, while for the streamwise and spanwise reducetosmallvaluesatthemiddleofthedomain.Zhangand
directions,itis∆x+<10,∆z+<10,respectively.Afterasta- Samtaney37computedthespanwisecorrelationsforspanwise
tistically steady behavior is reached, data are collected for domainsizesthatrangefrom0.1Cto0.8Candfoundthatthe
a period of up to 40 time units (a time unit is defined as results depend not only on the location selected but also on
tu=tU /C). the velocity component examined. In Figure 4 we compare
∞
InordertoassessthesuitabilityofgridresolutionforDNS thecomputedstatisticswithavailableresultsintheliterature
studies, we computed the ratio of the characteristic cell size for the small and large domain for the straight trailing edge
(definedasthecubicrootofthecellvolume)tothelocalKol- airfoil.
√
mogorovlengthscale,i.e., 3V.Contourplotsforthisratiofor TheDNSresultsofLehmkuhletal.26 areidealforsuch
η
theblunttrailingedgeandtheserratedtrailingedgeairfoil(in comparison.Figures4(a)and4(b)showthedistributionofthe
a plane through the peak) are shown in Figure 2. It must be pressureCp andskin-frictionCf coefficientsalongtheairfoil
mentionedthatthecellsizesintheXZplanearethesameforall surfaceforthesmallandlargedomains.Itcanbeseenthatthe
3configurations(thereareonlyminordifferencestoaccount effectofthedomainsize(betweenthetwodomainstested)is
for the different trailing edge geometry). As it can be seen, indeednegligible.
the ratio varies significantly and is maximal in the region of Asfarasthecomparisonbetweenthepresentresultsand
turbulentreattachmentandthewake.Inthelargestpartofthe thoseofLehmkuhletal.26isconcerned,therearesomesmall
domain, the ratio is around 1 (or smaller) and there are very discrepancies,especiallyintheCf coefficient,butoverallthe
smallregionsforwhichtheratioisaround2orslightlylarger matchingisreasonablygood.TheshapeoftheCp profilesis
(the max value is 2.5). Resolution requirements for a proper verysimilartothatmeasuredexperimentallyforlaminarsep-
DNShavebeenreportedintheworkofMoinandMahesh35 arationbubbleswithtransitionfollowedbyreattachment.38–40
andDonzisetal.36 Thelatterauthorsmentionthatastandard O’MearaandMueller38 identifythelocationofseparationas
resolution for the 2nd order quantities (such as rms veloci- thestartoftheconstant-pressureregion(alsoknownaspres-
tiesortheenergyspectrumwhichtakesverysmallvaluesand sure“plateau”),thefreeshearlayertransitionpointastheend
fallsoffrapidlyatwavenumbersk > 1/η)is∆x/η ≈ 2.Itis of the “plateau,” and the reattachment as the point at which
expectedthereforethatthepresentmeshissufficientlyfinefor pressure recovery exhibits a sharp decrease. These observa-
theresultstobetrueDNS. tions are broadly consistent with the Cp and Cf plots. The
InFigure3,thespanwiseautocorrelationsofthestream- pressure“plateau”appearstostartalittledownstreamofthe
wisevelocity,R ,are shown forthe straight(smalldomain) separationpoint, while thereattachement pointmatcheswith
uu
FIG.4. Comparison between the cur-
rentresultsforCp andCf coefficients
andtheDNSresultsfromtheworkof
Lehmkuhletal.26
014101-5 N.ThomareisandG.Papadakis Phys.Fluids29,014101(2017)
TABLEI. Comparisonoftime-averagedliftanddragcoefficientsCL,CD
andseparationandreattachmentpointsXsep,Xre.
Xsep Xre CL CD
Present 0.145C 0.564C 0.589 0.0271
ZhangandSamtaney37 0.141C 0.58C 0.562 0.0282
Lehmkuhletal.26 0.169C 0.566C 0.569 0.0291
the position at which the rate of pressure recovery changes.
InspectionoftheReynoldsstressdistributions(showninthe
left column of Figure 10) indicates that transition starts at
around X/C=0.45, slightly upstream of the end of the pres-
sure“plateau.”Theskinfrictionhasverysmallvaluesinside
therecirculationregionupstreamoftransition,whichiscon-
sistent with the “dead air region” of Horton.41 Following FIG.6. Qcontoursintherange[20,100]×(U∞/C)2forthestraighttrailing
transition, Cf reaches a minimum negative value, indicative edgecase.Thecolorscalevaluesarebasedonthestreamwisevelocity.
ofthe“reverseflowvortex.”41 Increasedmomentumtransfer
duetoturbulentmixingeventuallyeliminatesthereverseflow
III. EXTRACTIONOFDOMINANTMODES
andtheflowreattachesatX/C =0.564.
The separation and reattachment points as well as the Model reduction methods aim at reducing the complex-
forcecoefficientscomputedinthepresentworkarecompared ity of flows by extracting and analyzing their most domi-
with those from two other studies in Table I. There is good nantmodes(orstructures).ProperOrthogonalDecomposition
agreementwithpreviouslyreportedresults.Smalldiscrepan- (POD) is the most widely used reduction method.43 POD
ciesdoexist,butitmustbeborneinmindthatthetransition modes are mutually orthogonal and their amplitudes have
andreattachmentlocationsareverysensitivetothenumerical multi-frequency content. This indicates that the spatial vari-
discretisation and mesh resolution. The lift and drag coeffi- ationisdecoupledfromthetemporalvariationofthemodes.
cientsarealsoinreasonablygoodagreementwiththeresults Inthecurrentwork,weusetheDynamicModeDecompo-
ofLehmkuhletal.26andZhangandSamtaney.37 sition(DMD)methodproposedbySchmid.44,45 Themethod
Figures 5(a) and 5(b) show the variation of streamwise providesasetofnon-orthogonalmodes,eachhavingachar-
velocityandturbulentkineticenergyacrosstheboundarylayer acteristic frequency. The temporal and spatial variations are
at three locations after reattachment (results obtained using therefore fully coupled. This method is more suitable than
thesmallerdomain).Theprofilesarecomparedwiththoseof PODfortheflowfieldsexaminedinthepresentpaper,which
Lehmkuhletal.26andagainverygoodagreementisobserved. arecharacterizedbyasetoffewfundamentalfrequencies.Itis
In Figure 6 instantaneous iso-surfaces of the Q crite- thereforepossibletoextractdirectlythestructuresthatoscil-
rion are shown for the visualization of the flow separation late at these frequencies. DMD has already been applied to
and transition. The Q criterion was introduced by Jeong and differentflowssuchasswirlingjet46 andflowsaroundcylin-
Hussain42 for the visualisation of a vortex and it is defined ders of different diameters47 or around airfoils.48,49 A brief
asQ=1(u2 −u u )=1(||Ω||2−||S||2),whereΩandSare description of the algorithm will be presented below, but for
2 i,i i,j j,i 2
rotationandstrainratetensors,respectively.WhenQ>0,the amorerigoroustreatment,wereferthereadertotheoriginal
rotationratedominatesthestrainrateandservesasawayto worksofSchmid44,45andtheanalysisofJovanovic´etal.50The
identifyavortexcore.Theperiodicformationandbreakupof reviewpaperofBagheri51providesamoregeneraldiscussion
Kelvin-Helmholtzvorticesareclearlyseeninthefigure.The onmodelreductionmethods.
breakup,atthisparticulartimeinstant,startstoappearclearly We assume a series of N time snapshots of the veloc-
ataroundx/C=0.50,attheendofthepressure“plateau,”and ityfieldseparatedby∆t,XN = {v ,v ,...,v }.Eachvector
1 1 2 N
thestartoftherapidpressurerecovery. v is a column with the field data, for example, the u and (cid:51)
i
FIG.5. Comparisonofprofilesacrosstheboundarylayer.
014101-6 N.ThomareisandG.Papadakis Phys.Fluids29,014101(2017)
velocitycomponentsina2Dflow.Weassumethateachvector associated frequencies. For these structures, the growth rate
v hassizeM,usually M (cid:29)N.The methodassumesa linear (i.e., the magnitude of λ) should be very close to 1, which
i
dependencebetweenconsecutivesnapshotsoftheform, indicatesthattheyareneutral.Alargevalueoftheamplitude,
α, is not a reliable indicator for the selection of such struc-
vi+1 =Avi, (3) tuires,asitmaycharacteriseamodethatwilleventuallydecay.
Forthisreason,theamplitudesα aremultipliedby(λ )N−1,so
whereAistheunderlying(unknown)constantsystemmatrix i i
thatthemodeswithhighα but|λ | <1(i.e.,higherdamping
thatdescribesthedynamicbehaviorofthesystem.Ifthedif- i i
ferential evolution equation for variable v is dv=Bv, then ratios)arenotprominent.
dt
A=eB∆t and the aim of DMD is to extract the eigenvalues
and eigenvectors of matrix B using the field snapshots. If IV. RESULTSANDDISCUSSION
we define matrix XN−1 = {v ,v ,...,v }, then because
1 1 2 N−1 In this section, results are presented starting with the
of(3),
straighttrailingedge,followedbytheflatbackandendingwith
XN−1 = {v ,Av ,A2v ,...,AN−2v }. (4)
1 1 1 1 1 the serrated edge. This order of exposition helps understand
Combining(3)and(4)weget better the process of interaction of the trailing edge with the
separatingshearlayer.
XN =AXN−1, (5)
2 1
A. Airfoilwithstraighttrailingedge
whereXN = {v ,v ,...,v }.
2 2 3 N
The economy-size SVD (singular value decomposition) The DMD method was applied on a Z plane at the
ofXN−1 isXN−1 = UΣV∗ (∗ denotestheconjugatetranspose mid span, in an X-Y window with dimensions [−0.2,2]
1 1
ofamatrix),whereΣisadiagonal(r ×r)matrixthatcontains × [−0.27,018]C. In total N = 300 time snapshots separated
the r non-zero singular values (i.e., r is the rank of XN−1), by∆T = 0.01C/ wereused,whichcorrespondstoapprox-
1 U∞
and U and V are the matrices with ortho-normal columns imately 12 shedding cycles. A polar plot of the computed
(U∗U=IandV∗V=I)anddimensions(M×r)and(r×N), eigenvalues is shown in Figure 7(a). The modes appear in
respectively. complexconjugatepairsandtheoneswiththelargestampli-
If AisexpressedasA = UFU∗,thenitcanbeshown50 tudeα(λ )N−1aremarkedwithyellowcolour.Twodominant
i i
thatmatrixF = U∗XNVΣ−1 minimisestheFrobeniousnorm modes can be identified from the DMD spectrum shown in
2
(cid:107) XN −AXN−1 (cid:107)2. The eigenvector φ of F is related to the Figure 7(b): a low frequency mode located at f1 = 3.8 and a
2 1 F i
eigenvectoryiof Abyyi =Uφi.Thereconstructedflowfield highfrequencyatf2=7.4.
v˜ (j =1...N −1)isgivenby Inordertoverifythatthesefrequenciesareindeedpresent
j
and they are dominant, velocity time-signals from two dif-
(cid:88)r ferent locations at the mid span were analysed: one on the
v˜ = y (λ)jα, (6)
j i i i suctionsideoftheairfoilatX/C =0.53,Y/C =0.032(inside
i=1 thetransitioningshearlayerbutupstreamofreattachment)and
(cid:12)
where α is the amplitude of the ith mode and λ the corre- one in the wake at X/C = 2.0, Y/C = 0.17. In Figures 8(a)
i i
spondingeigenvalueofF.Inordertocomputetheamplitudes, and8(b)thepowerspectraldensityofthecross-streamveloc-
asecondoptimizationproblemissolved.50 ityisplotted.Thevelocitysignalatthefirstpoint(Fig.8(a))
When applied to a linear system, Equation (3) is exact. showsapeakatthehighfrequencyf observedinFigure7(b).
2
In this case, the computed modes/eigenvalues represent the Note that in Figure 8(a) the horizontal axis is linear and re-
physically correct structures and their growth rates, frequen- plottingusinglogarithmicaxis(figurenotshown)revealsthat
cies.Fornon-linearsystems,(3)representsabestlinearmap theenergycontentofthesignalisdistributedamongarange
that links all snapshots. If the flow is characterized by peri- offrequencies,centeredaroundabroadpeakatfrequencyf .
2
odicallyrepeatableandpersistent(i.e.,neutrallystable)struc- Suchabroadpeakhasalsobeenobservedexperimentally39,52
tures,theDMDmethodshouldbeabletodetecttheseandthe and is due to the fact that the shear layer amplifies a broad
FIG.7. Amplitude distribution of the
eigenvalues calculated from the DMD
algorithm.
014101-7 N.ThomareisandG.Papadakis Phys.Fluids29,014101(2017)
FIG.8. Power spectral density spec-
trumofthecrossstreamvelocitycom-
ponent V at two different locations,
X/C=0.53C,Y/C=0.032(a)andX/C
(cid:12)
=2.0C,Y/C= 0.17(b).
rangeoffrequencies.53Similarly,thevelocitysignalatX/C= etal.52haveidentifiedthestreamwisedistanceoftheshedvor-
2.0(Figure8(b))showsadominantpeaklocatedatthelower tices,ψ ,andthevelocityattheedgeoftheboundarylayerat
0
frequencyf ;againthisisabroadpeakduetothefactthatthe thepointofseparation,U ,astherelevantlengthandvelocity
1 es
pointisimmersedinaturbulentwake. scales,respectively.Thedistanceψ canbedirectlyestimated
0
Thespatialstructureofthetwomostdominantmodesas fromthespatialstructureoftheDMDmode(Figure9(b))and
wellasthetimeaveragedstreamwisevelocityfieldisdepicted itisfoundtobeψ ≈ 0.08C,whileU = 1.4U .Thecorre-
0 es ∞
inFigure9.Inthetimeaveragedfield,showninFigure9(a),a spondingnon-dimensionalfrequencyisf∗= f ψ /U ≈0.44,
2 2 0 es
thinrecirculationzonecanbeclearlyseenonthesuctionside veryclosetotheobservedrangeofvalues0.45–0.5.52
oftheairfoil.Thehighfrequency,f2=7.4,isassociatedwith Itisknownthatastheamplitudeoftheperturbationgrows
theKelvin-Helmholtzinstabilitytriggeredbytheinflectional spatially, non-linear interactions result in a lower frequency,
velocity profile across the separating shear layer.54 As Fig- which is a subharmonic of the fundamental instability fre-
ure9(b)shows,thismodeisdominantintheseparatingshear quency (for which the mechanism of the generation of the
layer, approximately between X/C = 0.4–0.8. In the region subharmonic refers to the review paper of Ho and Huerre53
whereitismostdominant,theshearlayerrollsupandsheds andreferencestherein).TheDMDmethodcapturesalowfre-
vortices,asshowninFig.6. quencyatf =3.8.Thefootprintofthisfrequencyisshownin
1
It is difficult to compare directly this frequency against Figure9(c).Themodeisactivatedatx ≈0.45C(slightlydown-
results from the literature because the flow conditions are streamcomparedwiththehighfrequencymode)andisseen
different. For example, Jones55 performed one-dimensional tobepresentinthereattachingshearlayerandthenearwake,
spatialstabilityanalysisonanumberoftime-averagedveloc- before decaying slowly further downstream. This frequency
ityprofilesalongthesuctionsideforthesameairfoil,angleof
attack,andReynoldsnumber,butwithMachnumberequalto
0.4,whichmakescompressibilityeffectsnon-negligible.The
frequency with the strongest spatial growth slightly changes
from one profile to the other, but overall the most amplified
frequencyintheshearlayerwasfoundtobef =8.49,whichis
closetothepresentf .Thedifferencecanbeattributedtotwo
2
reasons:first,Jones’sanalysisisone-dimensionalandisbased
ontheOrr-Sommerfeldequationand,second,compressibility
effectsareexpectedtobecomemoreimportantintheareaof
highaccelerationaroundtheleadingedgeandthereforeaffect
thevelocitydistributionoftheseparatingshearlayerandthe
mostamplifiedfrequencies.BoutilierandYarusevych39exam-
inedtheflowaroundaNACA0018atRe=100000andthey
also performed one-dimensional stability analysis using the
measured time-averaged profiles. For the 5◦ angle of attack,
thedominantfrequencywasfoundtobe13.0(TableIIintheir
paper).Thedifferencewiththepresentworkcanbeattributed
tothedifferentReynoldsnumber(thereisapowerlawdepen-
dency of the shear layer frequency and Reynolds number)40
andthedifferentthicknessoftheairfoilthataffectsthevelocity
distributioninthesuctionside.
Thepredictedfrequency, f ,howeverisinagreementwith
2
previousresultswhennon-dimensionalizedwiththeappropri-
FIG.9. Time-averagedstreamwisevelocityfield(a)andtwomostdominant
atelengthandvelocityscales.HuangandHo56andYarusevych DMDmodesforastraighttrailingedge(b)and(c).
014101-8 N.ThomareisandG.Papadakis Phys.Fluids29,014101(2017)
FIG.10. Streamwise,cross-stream,andspanwiseReynoldsstressdistributionsforthestraightandblunttrailingedgeairfoils.
is indeed very close to the subharmonic of the fundamental fluctuations appear denoting the rapid 3D breakdown of the
frequencyf /2=3.7(thedifferenceis2.7%). vortices.TwodimensionalsimulationsforthisReynoldsnum-
2
IntheleftcolumnofFigure10,contourplotsofthethree ber and angle of attack therefore do not describe accurately
normalcomponentsoftheReynoldsstressesareshown.The theflow.
samecolorscaleisusedtofacilitatecomparisonoftheplots.
RapidgrowthoffluctuationsstartsataboutX/C =0.45indi- B. Airfoilwithblunttrailingedge
catingtransitiontoturbulence,asalreadymentioned.Shortly
The DMD method was applied again on the mid-span
afterwardsarethelocationofthestartofrapidpressurerecov-
plane in the same X-Y window as in the straight trailing
eryshownintheC plot(Figure4(a))andthelocationofthe
p
edgecase,i.e.,[−0.2,2]×[−0.27,0.18]C.Intotal,240time
minimumvalueoftheskinfriction(Figure4(b)).Closetothe
snapshots, separated by ∆T=0.01C/ , were used, which
reattachmentpoint(X/C=0.56),allthreenormalcomponents U∞
correspond to 11 shedding cycles. The DMD spectrum in
attainverylargevalues.Afterthepeak,theReynoldsstresses
this window is shown in Figure 11. One very dominant fre-
decay downstream. Note also that strong spanwise velocity
quency appears at f =4.5. Its first harmonic at f =2f
bl,1 bl,3 bl,1
also appears. The shear layer natural frequency identified
previouslyisstillpresentatf =7.5,albeithighlyattenuated.
bl,2
FIG.12. Contourplotsofthetime-averagedvelocityfieldU/U∞(a)andthe
FIG.11. DMDspectrumoftheflowfieldaroundtheblunttrailingedgeairfoil. mostdominantflowmode(b)atfbl,1.
014101-9 N.ThomareisandG.Papadakis Phys.Fluids29,014101(2017)
(cid:12) (cid:12)
minimumC is 1.72forthestraighttrailingedge,and 1.59
P
fortheblunttrailingedge.Toidentifytheoriginofthisdiffer-
ence,potentialflowsimulationswithXfoil57 wereperformed
andexactlythesamebehaviorwasobserved;theoriginthere-
fore is inviscid. This difference has been also observed in
the past by Chen and Agarwal58 and it is a direct effect
of the truncation of the trailing edge, which increases the
radiusofcurvatureoftheleadingedgeinrelationtothechord
length.
Coherent structures originating from the truncated part
of the airfoil are clearly illustrated in Figure 12(b). This
mode, corresponding to the frequency peak f = 4.5, is a
bl,1
characteristicofbluffbodyvonKarmanvortexshedding.Well-
FIG.13. Comparisonofthepressurecoefficientsforthestraightandblunt
organized structures are shed from the truncated part of the
trailingedgecases.
airfoil and convect downstream without significant attenua-
tion,atleastinthecurrentfieldofview.TheStrouhalnumber
It is clear that the truncation of the trailing edge affects sig- ofthismodebasedonthetrailingedgebluntnessandthefree
nificantlytheseparatingshearflowanditsfrequencycontent; streamvelocityisSt =0.17.TheStrouhalnumbermeasured
bl,1
thisinteractionwillbeexaminedinmoredetailbelow. byNe´dicandVassilicos10was0.189foranidenticalgeometry
In Figure 12(a), contour plots of time-averaged stream- butwiththechordbasedReynoldsnumberthreetimeshigher
wise velocity are shown. The flow in the suction side is (Re=150000)andwithoutlaminarseparation(bothboundary
qualitatively similar compared to the straight trailing edge layersweretripped).KrentelandNitsche9foundtheshedding
airfoil:laminarseparation,transitiontoturbulence,andreat- Strouhalnumbertobe0.2,buttheReynoldsnumberbasedon
tachment. The locations of the separation and reattachment thetrailingedgethicknesswas44000,morethananorderof
pointsarehoweverdifferent(at0.23c and0.58c respectively magnitude higher compared to the one in the present paper
compared to 0.165c and 0.56c for the straight trailing edge) (1850).Althoughthismodehasthestrongestpresenceinthe
leadingtoanetreductionoftherecirculationzonesize.The wake,itsfootprintextendsfarupstreaminthesuctionsideof
most striking difference in the mean flows between the two theairfoil,ascanbeseeninFigure12(b).Thisisnotsurpris-
airfoilsisinthenearwake:thepresenceofbluntnesscreatesa ing.Thewakesheddingistheresultofaglobalflowinstability
smallrecirculationzone,whichisabsentinthestraightairfoil andthereforepresentinthewholeflowfield.Thiswakemode
case. isthereforepartiallycollocatedwiththeseparatingshearlayer
Inspection of Figure 13 demonstrates that the pressure andprovidesanexcitation(forcing)toit.Thishasimportant
distributionintheseparatingshearlayerisalsodifferent.The implicationsaswillbediscussedbelow.
FIG.14. Spectrafortheblunttrailing
edgeairfoilfordifferentcrossstreams
(cid:2) (cid:3)
(Y/C −0.3,0.1 )andstreamwise(X/C
=1.05,1.4,2.0)locations.