High fidelity numerical simulation of airfoil thickness and kinematics effects on flapping airfoil

21 Pages·2013·16.6 MB·English

by  Meilin Yu

JournalofFluidsandStructures42(2013)166–186
ContentslistsavailableatScienceDirect
Journal of Fluids and Structures
journal homepage: www.elsevier.com/locate/jfs
High fidelity numerical simulation of airfoil thickness and
kinematics effects on flapping airfoil propulsion
Meilin Yua,n, Z.J. Wanga, Hui Hub
aDepartmentofAerospaceEngineering,TheUniversityofKansas,Lawrence,KS66045,USA
bDepartmentofAerospaceEngineering,IowaStateUniversity,Ames,IA50011,USA
a r t i c l e i n f o a b s t r a c t
Articlehistory: High-fidelitynumericalsimulationswiththespectraldifference(SD)methodarecarried
Received8October2012 outtoinvestigatetheunsteadyflowoveraseriesofoscillatingNACA4-digitairfoils.Airfoil
Accepted2June2013 thicknessandkinematicseffectsonthe flapping airfoil propulsion arehighlighted.Itis
Availableonline17July2013
confirmedthattheaerodynamicperformanceofairfoilswithdifferentthicknesscanbe
Keywords: very different under the same kinematics. Distinct evolutionary patterns of vortical
Flappingairfoil structures are analyzed to unveil the underlying flow physics behind the diverse flow
Airfoilthickness phenomena associated with different airfoil thickness and kinematics and reveal the
Kinematics synthetic effects of airfoil thickness and kinematics on the propulsive performance.
HighorderNavier–Stokessimulations
ThicknesseffectsatvariousreducedfrequenciesandStrouhalnumbersforthesamechord
Spectraldifference lengthbasedReynoldsnumber(¼1200)arethendiscussedindetail.Itisfoundthatat
relatively small Strouhal number (¼0.3), for all types of airfoils with the combined
pitchingandplungingmotion(pitchangle201,thepitchaxislocatedatonethirdofchord
lengthfromtheleadingedge,pitchleadingplungeby751),lowreducedfrequency(¼1)is
conduciveforboththethrustproductionandpropulsiveefficiency.Moreover,relatively
thinairfoils(e.g.NACA0006)cangeneratelargerthrustandmaintainhigherpropulsive
efficiencythanthickairfoils(e.g.NACA0030).However,withthesamekinematicsbutat
relativelylargeStrouhalnumber(¼0.45),itisfoundthatairfoilswithdifferentthickness
exhibitdiverse trendonthrust production and propulsiveefficiency,especiallyat large
reducedfrequency(¼3.5).ResultsoneffectsofairfoilthicknessbasedReynoldsnumbers
indicate that relative thin airfoils show superior propulsion performance in the tested
Reynolds number range. The evolution of leading edge vortices and the interaction
betweentheleadingandtrailingedgevorticesplaykeyrolesinflappingairfoilpropulsive
performance.
&2013ElsevierLtd.Allrightsreserved.
1. Introduction
Unsteady flapping-wing aerodynamics has witnessedgreatprosperity in the last three decades. Thisprosperityis not
onlydrivenbythepowerofconsistentcuriosityabouttheeffectivepropulsionexhibitedbynaturalflyers,butalsodueto
theincreasinginterestinthedesignofMicroAirVehicles(MAVs),whichhavebeenconsideredtohavegreatpotentialto
open up new opportunities for surveillance-like missions in the near future. A large amount of experimental and
computational research has been carried out to reveal flapping-wing aerodynamics and many comprehensive reviews
nCorrespondingauthor.Tel.:[phone].
E-mailaddress:[email protected](M.Yu).
0889-9746/$-seefrontmatter&2013ElsevierLtd.Allrightsreserved.
http://dx.doi.org/10.1016/j.jfluidstructs.2013.06.001
M.Yuetal./JournalofFluidsandStructures42(2013)166–186 167
Nomenclature
C chordlength,m
CP pressurecoefficient,ðp−p∞Þ=ð0:5ρU2∞Þ R
C〈CTpow〉 ttihmruestavceoreafgfiecdiepnot,wTehrrucsote=fðf0ic:5ieρnUt2∞,−CðÞ1=TRÞ tt00þT½LiftðtÞy_ðtÞþMomentðtÞθ_ðtÞ(cid:2)dt=ð0:5ρU3∞CÞ
〈C 〉 timeaveragedthrustcoefficient,ð1=TÞ t0þTC ðtÞdt
〈CT 〉 timeaveragedthrustcoefficientcalculatt0edfroTmthepressureforce
T_P
〈C 〉 timeaveragedthrustcoefficientcalculatedfromtheviscousforce
T_v
f oscillationfrequency,Hz
H ,y plungeamplitudeandplungepositionofairfoil,m
0
k reducedfrequency,2πfC/U∞
Ma Machnumber
ReC Reynoldsnumberbasedonthechordlength,ρU∞C/μ
ReT Reynoldsnumberbasedontheairfoilthickness,ρU∞αC/μ
St Strouhalnumber,2fH0=U∞ R
ΔP averagedpressuredifferencesinthestreamwisedirection, ∂ΩpdS/FrontalArea
α thicknessratio(i.e.theratiobetweentheairfoilthicknessandthechordlength)
η propulsiveefficiency,〈C 〉=〈C 〉
pow T pow
θ ,θ pitchamplitudeandpitchangleofairfoil,deg
0
ϕ phaselagbetweenpitchingandplungingmotions,deg
andbookchapters(Hoetal.,2003;Platzeretal.,2008;RozhdestvenskyandRyzhov,2003;Shyyetal.,2010,1999,2008;
Triantaflyllouetal.,2004;Wang,2005)havesystematicallysummarizedtheprogressfromdifferentaspects.However,due
tothe intrinsic complexity of the unsteadyaerodynamics, there are still manyopen questions to be answered. One such
issueisthesyntheticeffectsofthewingthicknessandkinematicsonthethrustproductionandpropulsiveefficiencyduring
the flapping flight. As is widely recognized, natural flyers, e.g. birds and insects, maneuver their thin wings (i.e., wing
thicknessisonlyafewpercentofthechordlength)viathe“pitchingandrolling”motion.However,previousstudiesseldom
investigated the synergy of thin airfoils and sophisticated flapping kinematics but were limited to thick airfoils (airfoil
thicknesslargerthan10%ofchordlength)andrelativelysimplemotions(i.e.pitchingorplunging).Thesepreviousefforts
arebrieflyreviewedasfollows.
The first explanation for the thrust generation with oscillating airfoils was given by Knoller (1909) and Betz (1912)
independentlybasedontheinviscidassumptionandtheeffectiveangleofattack(AOA)concept.Thiswasexperimentally
confirmed by Katzmayr (1922) through mounting a stationary wing in an oscillating flow. After that, von Karman and
Burgers(1935)showedathought-provokingwaytoexplainthethrustordragproductionbycheckingmomentumsurfeitor
deficitinthewakebasedonthewakevorticesorientationandlocation.Thishasgraduallybecomeaqualitativeprincipleto
judgewhetheranoscillatingmechanismgeneratesthrustordrag.Astheaerodynamicperformanceoftheoscillatingfoils/
wings is closely tied to the vortical structures, the vortical structure analysis becomes an effective toolto unveil the key
mechanisminflapping-wingaerodynamics.Withdifferentflowvisualizationtechniques,Freymuth(1988),Koochesfahani
(1989),Jonesetal.(1998),LaiandPlatzer(1999),Godoy-Dianaetal.(2008),Schnipperetal.(2009),BohlandKoochesfahani
(2009)andYuetal.(2012),justtonameafewhere,haveexperimentallyinvestigatedthevortexstructuresaroundplunging
orpitchingfoils.Inparallel,byusingvariousnumericalsimulationortheoreticalanalysismethods,theevolutionofvortex
fields around plungingor pitching foils and the corresponding propulsion features (i.e. thrust production and propulsive
efficiency)havebeenstudiedbymanyresearchers(Triantaflyllouetal.,1993;Jonesetal.,1998;Wang,2000;LewinandHaj-
Hariri,2003;YoungandLai,2004;Yuetal.,2012).Notethatalltheaforementionedstudiesonwakestructuresarerelatedto
simpleairfoilkinematics,namelypureplungingorpitchingmotions,andrelativelylargeairfoilthickness.
Asreportedbymanyresearchers,acombinedpitchingandplungingmotionisamoreelaboratemodeloftheflappingwing
kinematics. Also, there is a growing awareness in this community that the airfoil thickness can affect thrust generation and
propulsive efficiency. Based on extensive numerical study, Jones and Platzer (1997) argued that the combined pitching and
plungingmotioncouldadjusttheeffectiveAOAoftheairfoil,enhancingthethrustgenerationorpowerextractionperformance.
Andersonetal.(1998)experimentallyexaminedthepropulsiveperformanceofaNACA0012airfoilwithacombinedpitchingand
plunging motion, and gave the parameters for optimal thrust production. Isogai et al. (1999), Tuncer and Platzer (2000),
Ramamurti and Sandberg (2001), Read et al. (2003), Schouveiler et al. (2005) and Kang et al. (2009) confirmed that better
propulsiveperformancewasachievedwhenthecombinedpitchingandplungingmotionwithsuitablephaselagwasadopted.
Thestudiesinairfoilthicknesshaveshownlargediscrepancyinitseffectsonaerodynamicperformance.LentinkandGerritsma
(2003)comparedtheperformancesofseveralplungingairfoilswithdifferentthicknessatReynoldsnumber150,andconcluded
that airfoils with larger thickness could generate greater thrust than the airfoils with small thickness. Based on the inviscid
simulation from an unsteady panel code, Cebeci et al. (2004) found that thickness had a negligible effect on the propulsive
efficiency.Anetal.(2009)foundthatthicknessratiowasacrucialparameterforthrustproductionandtheirresultsindicatedthat
168 M.Yuetal./JournalofFluidsandStructures42(2013)166–186
there existed an optimal thickness ratio for the thrust generation at Reynolds number 185. Ashraf et al. (2011) systematically
studiedthethicknesseffectsonthepropulsiveperformancesofflappingairfoilsatdifferentReynoldsnumbersrangingfrom200
to2(cid:3)106.TheyfoundthatatlowReynoldsnumbers,thinairfoilsoutperformedthickairfoils,whileathighReynoldsnumbers
thereexistedanoptimalthicknessforthethrustproductionandpropulsiveefficiency.Recently,Yuetal.(2012)foundthatatlow
Reynoldsnumbers,viscousforcecouldevengeneratethrustonplungingthinairfoils(e.g.NACA0004andNACA0006)duetothe
unsteadycircumnavigationofleadingedgevortices(LEVs).
Basedontheliteraturereview,itisclearthatstudiesonthesynergyoftheairfoilkinematicsandthicknessonflapping-
wing propulsion are limited, although this can be crucial for the MAV design. As is known, flapping-wing propulsion is
closelyrelatedtotheunsteadyvortexdynamics.Highordercomputationalfluiddynamics(CFD)methodsappeartobemore
accuratethantheirlowordercounterpartsforsimulationsofvortex-dominatedflow(Wang,2007).Somerecentworkon
thenumericalsimulationsofbio-inspiredflowwithhigh-ordermethods(e.g.Visbal,2009;Visbaletal.,2009;Perssonetal.,
2010;Liangetal.,2010;Yuetal.,2011;Ouetal.,2011)hasprovidedimpressiveresultsforthecommunityanddemonstrated
the potential of the high-order numerical simulation in unsteady flapping-wing aerodynamics. The simulations in the
presentpaperarecarriedoutusingthehigh-orderaccuratedynamicgridNavier–Stokesflowsolver.
Althoughallnatureflyersusecomplex3Dwingkinematicsintheflappingflight,someintrinsicprinciplesofunsteady
aerodynamicsdofallintothe2Dframeworksandunderstandingthesebasicphysicscanhelpmakeprogressonthestudyof
the3Dflappingwingaerodynamics.Therefore,thisstudyisonlyrestrictedto2Dflappingairfoilsimulations.Theremainder
ofthispaperisorganizedasfollows.InSection2,thedynamicunstructuredgridbasedspectraldifference(SD)methodis
first introduced. The simulationparametersare presented at the end of this section.Numericalresultsaredisplayed and
discussedinSection3.Thesyntheticeffectsofairfoilthicknessandkinematicsonflappingairfoilpropulsionareanalyzedin
detail, and the aerodynamic performance dependency on the reduced frequency and Strouhal number is also discussed
there.Section4brieflyconcludesthiswork.
2. Numericalmethodandsimulationparameters
2.1. Governingequations
NumericalsimulationsarecarriedoutwithanunsteadycompressibleNavier–Stokessolverusingdynamicunstructured
gridbasedhigh-orderSDmethoddevelopedinYuetal.(2011)andSunetal.(2006).The2DunsteadycompressibleNavier–
Stokesequationinconservationformreads,
∂Q ∂F ∂G
þ þ ¼0: ð1Þ
∂t ∂x ∂y
Herein,Q¼(ρ,ρu,ρv,E)Taretheconservativevariables,andF; G arethetotalfluxesincludingboththeinviscidandviscous
fluxvectors,i.e., F ¼Fi−Fvand G ¼Gi−Gv,whichtakesthefollowingform:
8 9 8 9
>>>><pþρuρu2>>>>= >>>>< ρρuvv >>>>=
Fi¼ ; Gi¼ ;
>>>>:uðEρuþvpÞ>>>>; >>>>:vpðEþþρvp2Þ>>>>;
8 9 8 9
>>>>< τ0 >>>>= >>>>< τ0 >>>>=
xx xy
Fv¼>>>>:uτ þvττyxþμCpT >>>>;; Gv¼>>>>:uτ þvττyyþμCpT >>>>; ð2Þ
xx yx Pr x xy yy Pr y
In Eq. (2), ρ is the fluid density, u and v are the Cartesian velocity components, p is the pressure, and E is the total
internalenergy,μisdynamicviscosity,C isthespecificheatatconstantpressure,PristhePrandtlnumber,andTisthe
p
temperature.ThestresstensorsinEq.(2)takethefollowingform
(cid:2) (cid:3) (cid:2) (cid:3)
u þv u þv
τ ¼2μ u − x y ; τ ¼2μ v − x y ; τ ¼τ ¼μðv þu Þ ð3Þ
xx x 3 yy y 3 xy yx x y
Onassumingthattheperfectgaslawisobeyed,thepressureisrelatedtothetotalinternalenergyby E¼ðp=ðγ−1ÞÞþ
ð1=2Þρðu2þv2Þwithconstantγ,whichclosesthesolutionsystem.
Toachieveanefficientimplementation,atime-dependentcoordinatetransformationfromthephysicaldomain(t,x,y)to
thecomputationaldomain(τ,ξ,η),asshowninFig.1(a),isappliedtoEq.(1).Thecoordinatetransformationisamapping
betweencomputationaldomainandthephysicaldomain,whichtakesthefollowingform,
t¼τ;
x¼xðτ;ξ;ηÞ;
y¼yðτ;ξ;ηÞ: ð4Þ
M.Yuetal./JournalofFluidsandStructures42(2013)166–186 169
Fig.1. (a)Transformationfromamovingphysicaldomaintoafixedcomputationaldomain;(b)Distributionofsolutionpoints(asdenotedbycircles)and
fluxpoints(asdenotedbysquares)inastandardquadrilateralelementforathird-orderSDscheme.
herein,(ξ,η)∈[−1,1]2,arethelocalcoordinatesinthecomputationaldomain.ReadersarereferredtothebookbyTannehill
etal.(1997)formoredetailsofthecoordinatetransformation.
ThenthetransformedequationofEq.(1)takesthefollowingform
∂Q~ ∂F~ ∂G~
þ þ ¼0; ð5Þ
∂τ ∂x ∂y
where
8
>><Q~ ¼jJjQ
>>:FG~~¼¼jjJJjjððQQξηtþþFFξηxþþGGξηyÞÞ : ð6Þ
t x y
Inthetransformationshownabove,theJacobianmatrixJtakesthefollowingform:
0 1
xξ xη xτ
∂ðx;y;tÞ B C
J¼ ∂ðξ;η;τÞ ¼@yξ yη yτA: ð7Þ
0 0 1
!
Notethatthegridvelocityvg¼ðxτ;yτÞisrelatedwith(ξt,ηt)by
8
<ξ ¼−!v ⋅∇ξ
:ηt¼−!vg⋅∇η: ð8Þ
t g
2.2. Spacediscretization
TheSDmethodisusedforthespacediscretization.IntheSDmethod,twosetsofpointsaregiven,namelythesolution
andfluxpoints,asshowninFig.1(b)fora2Dquadrilateralelement.Unknownsolutionsordegreesoffreedom(DOFs)are
definedatthesolutionpoints(SPs),andfluxesarecalculatedonfluxpoints(FPs).Inthepresentstudy,thesolutionpoints
arechosenastheChebyshev–Gaussquadraturepoints.ForaPN−1reconstruction,Nsolutionpointsareneededin1Dandare
specifiedas
(cid:2) (cid:3)
2s−1
ξ ¼−cos π ; s¼1;2;⋯;N: ð9Þ
s 2N
IthasbeenprovedinJameson(2010)thattheadoptionoftheLegendre–Gaussquadraturepointsasthefluxpointscan
ensure the stability of the SD method. Therefore, the flux points are selected to be the Legendre–Gauss points with end
points−1and1.Thesepointsaredenotedasξ ; f¼0;1;…;N.
f
TwosetsofLagrangepolynomialsbasedontheSPsandFPsrespectivelycanbespecifiedasfollows.
SPsbasedLagrangepolynomial:
N ξ−ξ
Ls;iðξÞ¼ ∏ ξ−ξs; i¼1;2;…;N: ð10Þ
s¼1;s≠i i s
FPsbasedLagrangepolynomial:
N ξ−ξ
Lf;iðξÞ¼ ∏ ξ−ξf ; i¼0;1;…;N: ð11Þ
f¼0;f≠i i f
170 M.Yuetal./JournalofFluidsandStructures42(2013)166–186
ThereconstructionoftheSDmethodisstatedbrieflyasfollows.Firstofall,theinviscidfluxesarereconstructed.Note
that the fluxes related to the grid movement are incorporated into the inviscid fluxes, e.g., F~i¼jJjðQξτþFiξxþGiξyÞ. The
conservativevariablesQ onthefluxpointsareinterpolatedfromtheconservativevariablesQ onthesolutionpointsviaa
f s
tensorproductionofthe1DLagrangepolynomialsEq.(11),whichtakesthefollowingform
N N
Qfðξ;ηÞ¼ ∑ ∑ Qsðξi;ηjÞLs;iðξÞLs;jðηÞ: ð12Þ
j¼1i¼1
Then the fluxes can be reconstructed at the flux points using Q. Note that this reconstruction is continuous within a
f
standardelement,butdiscontinuousonthecellinterfaces.Therefore,aRiemannfluxorcommonfluxneedstobespecified
ontheinterfacetoensureconservation.Sincetheflowregimeforflappingflightisalmostincompressibleandthepresent
governing equations are compressible N-S equations, the Riemann solver should retain good performance at low Mach
numbers.TheAUSM+-upRiemannsolver(Liou,2006)forallspeedsisimplementedforthepresentsimulationandisproved
tobehavewellatlowMachnumbers.TheprocedurestoreconstructthecommonfluxesfromtheAUSM+-upRiemannsolver
arestatedasfollows.
!
Denotethefacenormalofarbitraryinterfaceby n,thentheinterfacemassflowratem_1=2 reads
(
m_1=2¼a1=2M1=2 ρρL iofthMe1r=w2i4se0; ð13Þ
R
wherethesubscript‘1/2’standsfortheinterface,aandMarespeedofsoundandtheMachnumberrespectively.Notethat
~i ~i
thegridvelocityhasbeenincludedintheinterfaceMachnumberM.ThenumericalnormalfluxesF andG canthenbe
specifiedas
8 ( !
>>>>><F~i¼ m_1=2 ψψL iofthm_e1r=w24ise0 þP1=2 jJjj∇ξjsignð!nU∇ξÞ;
>>>>>:G~i¼ m_1=2(ψψRL iofthm_e1r=w24ise0 þP1=2!jJjj∇ηjsignð!nU∇ηÞ; ð14Þ
R
(cid:4) (cid:5)
whereψ¼ð1;u;v;ðEþpÞ=ρÞT; P¼ 0;pnx;pny;pvg;n T,withnx,nyandvg,nspecifyingthefacenormalcomponentsinxand
ydirectionandthegridvelocityinthefacenormaldirection,respectively.
Afterthis,thederivativesoftheinviscidfluxesarecalculatedonthesolutionpointsusingthefollowingformulas:
~i
∂∂Fξðξ;ηÞ¼ ∑N ∑N F~iðξi;ηjÞLf0;iðξÞLs;jðηÞ;
j¼1i¼0
~i
∂G∂ηðξ;ηÞ¼ ∑N ∑N G~iðξi;ηjÞLs;iðξÞLf0;jðηÞ: ð15Þ
j¼0i¼1
Since the viscous fluxes are functions of both the conservative variables Q and their derivatives ∇Q, slightly more
involvedreconstructionproceduresareneeded.Inthepresentstudy,theapproachproposedinBassiandRebay(1997),also
known!as‘BR1’,isadopted.TheimplementationofthisapproachinSDisbrieflyintroducedasfollows.
Let R ¼∇Q.!On transforming this formula from the physical domain to the computational domain, we obtain the
componentsof R intheconservationformas
(cid:2) (cid:3)
1 ∂jJjQξ ∂jJjQη
Rx¼ xþ x ;
jJj ∂ξ ∂η
(cid:2) (cid:3)
1 ∂jJjQξ ∂jJjQη
Ry¼ yþ y : ð16Þ
jJj ∂ξ ∂η
thenthederivativesinEq.(16)onthesolutionpointscanbecalculatedfromtheconservativevariablesQ onthefluxpoints
f
followingtheprocedureasshowninEq.(15).NotethatthecommonconservativevariablesQcomonelementinterfacesare
usedinthederivativecalculation.InBR1,Qcomistheaverageoftheleftandrightsolutionsontheinterface,
QLþQR
Qcom¼ : ð17Þ
2
Afterthis,thegradientofQistheninterpolatedbacktofluxpointsfollowingtheprocedureasshowninEq.(12)andthe
viscous fluxes can then be calculated on flux points. Again, the gradient of Q from the aforementioned reconstruction is
generally discontinuous on the element interface, and BR1 is used to provide a common gradient ∇Qcom on the element
interfaceasfollows,
∇QLþ∇QR
∇Qcom¼ : ð18Þ
2
~v ~v
ThustheviscousfluxesF andG onfluxpointsareuniquelyspecifiedinalocalcell,andthefluxderivativesonsolution
pointscanthenbecalculatedviatheapproachasshowninEq.(15).
M.Yuetal./JournalofFluidsandStructures42(2013)166–186 171
Fig.2. SketchofthekinematicsoftheNACA0012airfoil.
Fig.3. (a)ConvergencehistoryoftheenergyresidualforthesteadysolutionoftheviscousflowoverastationaryNACA0012airfoilwithimplicit(LU-SGS)
timeintegrationatRe¼5000,Ma∞¼0.05andzeroAOA;(b)pressurecoefficientcontoursfortheconvergedsteadyflow;(c)Machnumbercontoursforthe
convergedsteadyflow.
Onceallfluxderivativesareavailable,theDOFscanbeupdatedwitheitherexplicitorimplicittimeintegrations.
2.3. Simulationparameters
Inthisstudy,thekinematicsoftheairfoilisspecifiedasfollows.
Plungingmotion:y¼H sinð2πftÞ;
0
Pitchingmotion:θ¼θ sinð2πftþϕÞ: ð19Þ
0
Herein, y is the plunge position of the airfoil, H is the plunge amplitude, f is the oscillation frequency, θ is the pitch
0 0
amplitude, and ϕ is the phase lag between pitching and plunging motions. In view of the optimal thrust generation
conditionssuggestedbyAndersonetal.(1998),thepivotpointforthepitchingmotionislocatedone-thirdchordlength
awayfromtheleadingedge,θ issetas201andϕisfixedat751.SketchoftheairfoilkinematicsisdisplayedinFig.2.Grid
0
deformationstrategyassociatedwiththekinematicsoftheairfoilcanbefoundinthepaperbyYuetal.(2011).
Theaerodynamicparametersarespecifiedasfollows.TheReynoldsnumber(Re )basedontheairfoilchordlengthCand
C
thefreestreamvelocityU∞issetas1200formostofsimulations.Meanwhile,theairfoilthicknessbasedReynoldsnumber
(Re )isalsousedtocharacterizedifferentflowfeatures.AlltheseReynoldsnumbersfallwellwithintheinsectflightregime.
T
According to the optimal Strouhal number suggested by Triantaflyllou et al. (1993), the Strouhal number (St¼2fH0/U∞)
basedontheplungeamplitudeH0ischosenas0.3or0.45.Thereducedfrequency(k¼2πfC/U∞)basedontheairfoilchord
lengthCisassignedas1.0or3.5.Asreportedbysomeresearchers(YoungandLai,2004),theinletMachnumberMa∞is
crucialfortheaerodynamicperformanceofaflappingwingifacompressiblesolverisadoptedforthesimulation.Inthe
presentstudy,Ma∞issettobe0.05,whichisinaccordancewiththevaluespecifiedinYoungandLai(2004).
3. Resultsanddiscussion
TheperformanceofthedevelopedsolverforlowMachnumberflowisfirsttestedbycomputingthesteadyviscousflow
overaNACA0012airfoilatRe¼5000,Ma∞¼0.05andzeroAOAwitha3rdorderaccurateSDschemeandanimplicitLU-SGS
timeintegration(Yuetal.,2011)onacoarsemesh.Theresidualconvergencehistory,pressurecoefficient(C )contour,and
P
theMachnumbercontouraredisplayedinFig.3(a)–(c)respectively.ItisobviousthattheNavier–Stokessolverworkswellat
lowMachnumberasnooscillationoftheflowfieldisobserved.
Then h-refinement (grid refinement) and p-enrichment (polynomial degree enrichment) studies are carried out to
determinethesuitablegridandschemeaccuracy,whichensurethattheaerodynamicforceandvorticalstructuresaroundthe
172 M.Yuetal./JournalofFluidsandStructures42(2013)166–186
airfoil are grid and scheme accuracy independent for the simulations of the oscillating airfoils. The case of an oscillating
NACA0006airfoilwiththecombinedpitchingandplungingmotionsatk¼3.5,St¼0.3isusedforthehp-refinementstudy.
TwosetsofgridsinthevicinityoftheairfoilaredisplayedinFig.4(a)and(b).Forthecoarsemesh,100gridpointsareplaced
ontheairfoilsurfaceand35gridpointsareassignedintheradialdirectionfromtheairfoil.Thenumberofgridpointsforthe
finemeshinbothdirectionsis1.5timesthatforthecoarsemesh.AnoverviewofthecomputationalgridisshowninFig.4(c)
with the far field boundary imposed at about 10 chord length away from the airfoil surface. Time histories of the thrust
coefficients(C)fromthehp-refinementstudyareshowninFig.4(d).Itisfoundthattimehistoriesofthrustcoefficientsforthe
T
3rdorderschemewithbothcoarseandfinegridsandthe4thorderschemewiththecoarsegridagreewellwitheachother,
while thatforthe 2ndorder schemewith thecoarsemeshshows marked deviation.Based on theseresults,the 3rdorder
accurate SD scheme with a coarse mesh is chosen to conduct the simulations. Moreover, the non-dimensional time step
(Δt¼ΔtU∞=C)usedinthesimulationsissetas1.2(cid:3)10−5,whichisdeterminedfromatimerefinementstudyonthecasewith
thechosenaccuracyandgridconfiguration.MorevalidationofthepresentflowsolvercanbefoundinYuetal.(2012,2011).
Nextitisfirstlyconfirmedthattheaerodynamicperformanceofairfoilswithdifferentthicknesscanvarydramatically
underthesamekinematics,andthatthecombinedpitchingandplungingmotionwithsuitablephaselagisconducivefor
enhancingtheaerodynamicperformance.Thenweexaminetheairfoilthicknesseffectsonflappingairfoilpropulsionwith
the combined motion at different reduced frequencies, the Strouhal numbers and the Reynolds numbers and unveil the
hiddenflowphysicsviavortexdynamicsanalysis.
3.1. Syntheticeffectsofairfoilthicknessandkinematics
Three kinds of oscillation motions, namely the pure plunging motion, the pure pitching motion and the combined
pitchingandplungingmotion,arestudiedinthissection.Thetimehistoriesof thethrustcoefficientsforNACA0006and
NACA0030withdifferentoscillationmotionsatk¼1;St¼0.3aredisplayedinFig.5.
BycomparingthetimehistoriesofthethrustcoefficientsinFig.5andtimeaveragedthrustcoefficientsinTable1,several
observationscanbesummarizedasfollows.
Fig.4. hprefinementstudyfortheflowovertheNACA0006airfoilwiththecombinedpitchingandplungingmotionatk¼3:5;St¼0.3.(a)Thecoarsegrid
and(b)finegridusedinthegridrefinementstudy.(c)Overviewofthecomputationalgrid.(d)Timehistoriesofthethrustcoefficientsfromthehp
refinementstudy.
M.Yuetal./JournalofFluidsandStructures42(2013)166–186 173
Fig.5. Timehistoriesofthethrustcoefficientsfor(a)theNACA0006and(b)NACA0030airfoilswithdifferentmotionsatk¼1,St¼0.3.
Table1
TimeaveragedthrustcoefficientsfortheNACA0006andNACA0030airfoilswithdifferentmotionsatk¼1; St¼0.3.〈CT〉denotesthetotaltimeaveraged
thrustcoefficient;〈CT_P〉denotesthetimeaveragedthrustcoefficientcalculatedfromthepressureforce;〈CT_v〉denotesthetimeaveragedthrustcoefficient
calculatedfromtheviscousforce.
NACA0006 NACA0030
Plu. Pit. Com. Plu. Pit. Com.
〈CT〉 0.0085 −0.240 0.418 0.020 −0.399 0.105
〈CT_P〉 0.032 −0.198 0.453 0.056 −0.327 0.166
〈CT_v〉 −0.0235 −0.042 −0.035 −0.036 −0.075 −0.061
First,itisfoundthatforbothairfoils,(1)althoughthecombinedmotionisalinearsuperpositionoftheplungingand
pitching motion, the aerodynamic force generated from the combined motion is not the linear superposition of the
aerodynamicforcegeneratedfromtheplungingandpitchingmotion;(2)onlydragisgeneratedwhenpurepitchingmotion
isused;(3)thecombinedmotionoutperformspureplungingmotion.
Second,NACA0006cangeneratelargerthrustthanNACA0030forthecombinedpitchingandplungingmotion.However,
for pure plunging motion, NACA0030 outperforms NACA0006. The pitching motion did not help much in enhancing the
aerodynamic performance of the oscillating NACA0030 airfoil, but it significantly increases the thrust generation of
NACA0006aswhenpureplungingmotionisconducted,theNACA0006airfoilalmostdoesnotgeneratethrust.
Theexplanationforthecommonfeaturesdisplayedbybothairfoilsisgivenasfollows.Thenonlinearsuperpositionof
aerodynamicforcesisattributedtothenonlineareffectsoftheunsteadyvorticalflow.Avisualproofofthisnonlinearityis
thatthevorticalstructuresfromthecombinedmotionarenotlinearsuperpositionofthosefrompureplungingandpitching
motionasclearlydisplayedinFigs.6and8.AccordingtothedefinitionoftheStrouhalnumberforthepitchingmotionas
showninYuetal.(2012),itsvalueinthepresentcaseisStpit¼2f(2Csinθ0/3)/U∞¼0.07.Thisvalueiswellbelowthecritical
Strouhalnumberforthrustproductionanditindicatesthatpurepitchingmotionisnotpowerfulenoughtogeneratethrust.
Instead,theplungingmotionplaysamajorroleinthrustgenerationastheplungingamplitudebasedStrouhalnumberis0.3
here.NotethatthrustgenerationviathepitchingmotionhasbeenpreviouslystudiedbyKoochesfahani(1989),Jonesand
Platzer (1997), Ramamurti and Sandberg (2001) and Bohl and Koochesfahani (2009), just to name a few. There, similar
propulsionperformancefromthepitchingmotionhasbeenreported.Beforeexplainingwhymorethrustisgeneratedfrom
the combined motion than pure plunging one, we recognize that the superior propulsion performance of the combined
motionhavebeenstudiedpreviouslyfromdifferentaspects,suchastheeffectsofthephaselagbetweenplungeandpitch,
theeffectsofAOAandtheeffectsofreducedfrequencyandStrouhalnumberbymanyresearchers(e.g.JonesandPlatzer,
1997;Andersonetal.,1998;TuncerandPlatzer,2000;RamamurtiandSandberg,2001;Readetal.,2003;Schouveileretal.,
2005).Itisreportedfromtheseworkthataerodynamicparametersleadingtogoodperformanceofthecombinedmotion
are generally associated with the favorable dynamic vortex configuration around the airfoils. Here, a comprehensive
explanation for the superior performance of the combined motion is given by relating the aerodynamic force with
kinematics and vorticity fields. First, weintroduce the conceptof ‘frontal area’, which is defined as the projection of the
airfoilinthefreestreamdirectionasshowninFig.2.Obviously,whentheairfoilisparalleltothefreestream,thefrontal
area is actually the airfoil thickness. From Table 1 we find that the drastic change in the thrust coefficients between the
combinedmotionandpureplungingmotionismainlyduetothechangeofthepressureforcecontributiononthethrust.
Recallthatthiscontributionisdeterminedbothbypressuredifferenceinthefree-streamdirectionandfrontalarea.From
thevorticityandpressurefields withdifferentmotionsatdifferentphasesasdisplayed inFigs.6and7(NACA0006) and
174 M.Yuetal./JournalofFluidsandStructures42(2013)166–186
Fig. 6. Vorticityfields for the NACA0006airfoil withdifferentmotionsatk¼1; St¼0.3.(a)–(d)combined pitchingand plunging motion atphases
0; π=2;π;3π=2 respectively; (e)–(h) pitching motion at phases 0; π=2;π;3π=2 respectively; (i)–(l) plunging motion at phases 0; π=2;π;3π=2
respectively.
Figs.8and9(NACA0030),itisclearthatwhenasuitablepitchingmotionisaddedtotheplungingmotion,thefrontalarea
canbeeffectivelyadjustedtoenhancethethrustgenerationfrompressureforce.Forexample,atphases0andπ,thelow
pressureregioninducedbyLEVsisattachedonthewindwardsideofbothairfoils.Atthesamephases,thefrontalareaof
both airfoils almost reaches its maximum value. Thus, the performance in terms of thrust generation is significantly
improvedwiththecombinedflappingmotion.
Now we explain how airfoil thickness affects the aerodynamic performance with different kinematics. The inferior
performanceoftheplungingNACA0006airfoilisattributedtothesmallfrontalarea,whichisalsotheairfoilthicknessin
this case. Specifically, it is observed from Figs. 7(i)–(l) and 9(i)–(l) that at phases 0 and π, which are close to the thrust
generationpeaks,bluntleadingedgeofNACA0030takesfulladvantageofthelowpressureregioninducedbytheleading
edgevortices.FromFig.5itisclearthatthethrustpeakfortheplungingNACA0030airfoilisaboutfourtimesofthatfor
NACA0006.ItisalsofoundfromFig.5thatthedragpeakfortheplungingNACA0030airfoilisonlyabouttwiceofthatfor
NACA0006.ThisexplainswhyNACA0030outperformsNACA0006forpureplungingmotion.Theseresultsalsoagreewith
theobservationsbyYoungandLai(2004)thatleading-edgeeffectsarevitalindeterminingaerodynamicforces.
ThereasonthatthepitchingmotionalmostdoesnotenhancethethrustgenerationoftheoscillatingNACA0030airfoilis
givenasfollows.OnexaminingFig.8(a)and(c),itisobservedthatalthoughthepitchingmotionhelpsenlargethefrontal
areaoftheNACA0030airfoilatphases0andπ,i.e.t/T¼3and3.5inFig.5(b),theLEVsisstilllocatedintherearwardpartof
theairfoil.Thisvortexconfigurationdoesnotcontributemuchtothethrustgeneration.Moreover,atphasesπ/2and3π/2,i.e.
t/T¼3.25and3.75inFig.5(b),theNACA30airfoilwiththecombinedmotionexperienceslargerdragpeakthanthatwiththe
plungingmotion.ThisisduetotheinteractionofthepropagatingLEVsandtheshedtrailingedgevortices(TEVs)duringthe
oscillationstroke.Specifically,forthecombinedmotion,partsoftheTEVsaretrappedonthetrailingedgeasshowninFig.8
(b)and(d),inducingalowpressureregionthereasshowninFig.9(b)and(d).However,fortheplungingmotion,theTEVs
arenottrappedonthetrailingedge,asdisplayedinFig.8(j)and(l).Thustheredoesnotexistalowpressureregioninthe
rearwardpartoftheairfoil.
M.Yuetal./JournalofFluidsandStructures42(2013)166–186 175
Fig. 7. Pressure fields for the NACA0006 airfoil withdifferent motions at k¼1; St¼0.3. (a)–(d) combined pitching and plunging motionat phases
0; π=2;π;3π=2 respectively; (e)–(h) pitching motion at phases 0; π=2;π;3π=2 respectively; (i)–(l) plunging motion at phases 0; π=2;π;3π=2
respectively.
ThebetterperformanceoftheNACA0006airfoilwiththecombinedmotionisattributedtotheLEVsdynamicsaswell.
FromFigs.6and7,itisobservedthatfortheNACA0006airfoilthepitchingmotioncanbettercontroltheLEVsandensure
thatthelowpressureregionislocatedontheforwardpartoftheairfoilatphases0andπ,andnolowpressureregionis
locatedontherearwardpartoftheairfoilatphasesπ/2and3π/2.FromFig.5(a),itisfoundthatincontrasttoNACA0030,the
NACA0006airfoilwiththecombinedmotiondoesnotexperiencelargerdragpeakthanthatwiththeplungingmotion.This
isduetothatnoTEVsaretrappedonthetrailingedgeforboththecombinedandplungingmotions.
3.2. ThicknesseffectsatdifferentreducedfrequenciesandStrouhalnumbers
Asaforementioned,asuitablecombinedmotioncanenhancetheaerodynamicperformancesoftheoscillatingairfoil.In
thissection,thekinematicsoftheairfoilsischosenasacombinedmotionwiththepitchingmotionleadingtheplunging
motionby751.ThicknesseffectsatdifferentreducedfrequenciesforagivenStrouhalnumberarecomparedanddiscussedat
first.TheevolutionofLEVsandtheinteractionbetweenLEVsandTEVsisconfirmedtobevitalfordifferentperformanceof
airfoilswithdifferentthickness.Thephysicsbehindthesephenomenaisunveiledbyanalyzingthevorticityfieldsandthe
correspondinginducedpressuredistribution.Thensimilaranalysisstrategyisadoptedtorevealthicknesseffectsatdifferent
Strouhalnumbersforthesamereducedfrequency.
3.2.1. Thicknesseffectsatdifferentreducedfrequencies
The flow fields, namely the spanwise vorticity fields and the corresponding pressure fields of NACA0020 at k¼1 and
St¼0.3aredisplayedinFig.10.NotethattheflowfieldsofNACA0012aresimilartothoseofNACA0006,andarenotshown
here.BycomparingthevorticityfieldsofNACA0006(Fig.6(a)–(d)),NACA0020(Fig.10(a)–(d))andNACA0030(Fig.8(a)–(d)),
itisobservedthatastheairfoilbecomesthicker,attheendoftheoscillationstrokes,i.e.atphasesπ/2and3π/2,largerparts