Simulation of Transitional Flow over Airfoils Using the Spectral Difference Method

17 Pages·2010·1.32 MB·English

AIAA 2010-4626
40th Fluid Dynamics Conference and Exhibit
28 June - 1 July 2010, Chicago, Illinois
Simulation of Transitional Flow over Airfoils using the Spectral
Difference Method
PatriceCastonguay∗, ChunleiLiang †, AntonyJameson‡
DepartmentofAeronauticsandAstronautics,StanfordUniversity,Stanford,CA,94305
ThisworkaddressesthesimulationoftransitionalflowoverairfoilsunderlowReynoldsnumberconditions
(Rec ≤60000). TheflowsolutionsareobtainedbymeansofanImplicitLargeEddySimulation(ILES)using
anewlydevelopedunstructured,parallelsolverthatemploysthehigh-orderspectraldifference(SD)method
forspatialdiscretization. ThecalculationsareperformedontheSD7003airfoilsectionatanangleofattack
of 4◦ at Reynolds number of 10000 and 60000. The SD7003 airfoil was selected due to the availability of
experimental and computational results. The use of the SD method without an added subgrid-scale model
appearstobecapableofaccuratelypredictingthelaminarseparation,transitionandreattachmentlocations.
To the authors’s knowledge, the present studyis the first attempt to analyze transitional flow usingthe SD
method.
Nomenclature
Cp Timeandspanwiseaveragecoefficientofpressure,(p −p¯)/(1ρ U2 )
∞ 2 ∞ ∞
C Timeandspanwiseaverageskinfrictioncoefficient,2 µ ∂u
f Re∂n
C Specificheatconstantatconstantpressure
p
F ,G ,H Inviscidfluxvectors
I I I
F ,G ,H Viscousfluxvectors
V V V
J Jacobianmatrixofcoordinatetransformation
M Machnumber
N OrderofaccuracyofSpectralDifferencescheme
Pr Prandtlnumber(0.72forair)
Q Statevector
R Gasconstant
Re FreestreamReynoldsnumber,ρ U C/µ
∞ ∞ ∞
T Non-dimensionaltemperature
U Freestreamvelocity
∞
c Chordlength
k Coefficientofthermalconductivity
p Non-dimensionalstaticpressure
q Componentofheatfluxvector
i
t Non-dimensionaltime
u,v,w Non-dimensionalCartesianvelocitycompoenntsinx,y,zdirections
u′v′ Reynoldsshearstress
x,y,z Non-dimensionalCartesiancoordinates
α Angleofattack
δ Kroneckerdeltafunction
γ Specificheatratio(1.4forair)
∗Ph.D.Candidate,DepartmentofAeronauticsandAstronautics,StanfordUniversity,AIAAStudentMember
†PostdoctoralScholar,DepartmentofAeronauticsandAstronautics,StanfordUniversity,AIAAMember
‡ThomasV.JonesProfessorofEngineering,DepartmentofAeronauticsandAstronautics,StanfordUniversity,AIAAMember
1of17
Copyright © 2010 by P. Castonguay and A. Jameson. Published bAy mtheer AicmanerIincastnit uIntsetiotuftAe eorf oAnearuotnicasutaicnsd aAnsdt rAosntarountiacustics, Inc., with permission.
ρ Non-dimensionaldensity
σ Viscousstresstensor
µ Non-dimensionaldynamicviscosity
ξ,η,ζ ReferenceelementCartesiancoordinates
I. Introduction
AccuratesimulationoflowReynoldsnumberflowsoverairfoilsandwingsisofcrucialimportanceforthedesign
of model airplanes and micro air vehicles (MAV). Due to their size and speed, the Reynolds numbers over these
vehicles are typically of the order of 104 to 105. Accurate performance prediction in this Reynolds number range
is difficult due to the tendencyof the laminar flow to separate and transition overthe wing. At moderate angles of
attack, even the mildest pressure gradientcauses the laminar flow to separate overthe airfoil. After separation, the
flow structures become increasingly irregularand eventually, transition from laminar to turbulentflow occurs. The
turbulentmixingbringshigh-momentumfluidfromthefreestreamtothenear-wallregion,whichcausestheflowto
reattach,thusformingalaminarseparationbubble(LSB).Figure1showsthemainfeaturesofaLSB1.
Figure1.Generalfeaturesofthemainflowfieldinatransitionallaminarseparationbubble(LSB)1
Asthe Reynoldsnumberor angleof attack is increased, the transitionlocation movesfurtherupstream, thereby
diminishing the size of the laminar separation bubble. At some critical Reynolds number, the transition location
coincideswiththeseparationlocationandthebubbledoesnotform.
Oneofthemostpopularapproachtonumericallysimulatetransitionalflowoverairfoilsistocouplea Reynolds
Averaged Navier-Stokes (RANS) solver with a turbulence transition model. The flow is assumed to be laminar up
tothetransitionlocationpredictedbythemodel,thusallowingalaminarseparationbubbletoform. Downstreamof
the transition location, the turbulencemodelis switched on and causesthe flow to reattach. This approachrequires
anefficientandaccuratetransitionmodelinordertoobtainphysicallymeaningfullresults. Oneofthemostpopular
transition modelis the eN method, which is based on linear stability analysisand boundarylayer theory. With this
approach,theOrr-Sommerfeldequationissolvedinordertopredictthelocalgrowthratesofunstablewaves.Recently,
researchers2–4havecoupledaRANSsolverwiththeeN transitionmodelandhavebeenabletoaccuratelycapturethe
time-meanLSBuptostall. However,resultsobtainedwiththeRANS-eN approachareverydependentonthecritical
N-factorusedinthetransitionalmodelandthechoiceofturbulencemodel. Furthermore,thisapproachislimitedby
2of17
AmericanInstituteofAeronauticsandAstronautics
itsassumptionsoftwo-dimensionalparallelsteadyvelocityprofilesandthinboundarylayers.
Recentstudies5–7 haveinvestigatedtheapplicationofhigh-orderspatialdiscretizationschemeswithoutanadded
subgrid-scale(SGS)modelforcompressiblelarge-eddysimulations.Intheseso-calledImplicitLES(ILES)schemes,
dissipation from a filter or from the numerical flux serves to incorporate the additional dissipation associated with
theunder-resolvededdies. Forexample,VisbalandRizetta5 useda6th orderfinitedifferenceschemecoupledwith
a high-order low-pass filter to simulate an isotropic decaying turbulent flow. Results were in very good agreement
with those obtainedfrom a DNS calculation by Spyropoulosand Blaisdell8. ILES approacheshave also been used
by Urangaet al.6 and Galbraith and Visbal7 to simulate transitionalflow overthe SD7003and have producedvery
goodresults. TheILESapproachisveryattractivebecauseitallowstheuseofasingleframeworkformixedlaminar,
transitionaland turbulentthree-dimensionalflows. Furthermore,it doesnotrequireany parametersto be tunedand
isnotlimitedtoflowswithtwo-dimensionalparallelsteadyvelocityprofiles. Thepresentworkaimsatinvestigating
the feasibility of an ILES approach that uses the high-orderSD method for spatial discretization and an artificially
dissipativenumericalflux.
TheSDmethodisanewlydevelopedefficienthigh-orderschemebasedonthedifferentialformofthegoverning
equations.ThefoundationfortheSDschemewasfirstputforwardbyKoprivaandKolias9in1996underthenameof
“staggeredgridChebyshevmultidomain”methods. However,severalyearslaterin2006Liu,WangandVinokular10
presentedamoregeneralformulationforbothtriangularandquadrilateralelements,whichtheytermedtheSDmethod
(anamewhichhasbeenretainedtothepresent).Wangetal.11extendeditto2DEulerequationsontriangulargridsand
Sunetal.12furtherdevelopeditforthethree-dimensionalNavier-Stokesequationsonhexahedralunstructuredmeshes.
TheSDmethodcombineselementsfromfinite-volumeandfinite-differencetechniques,andisparticularlyattractive
becauseitis conservative,hasa simpleformulationandprovidesgeometricflexibility. Similar to the discontinuous
Galerkin method13, the SD scheme achieves high-orderaccuracy by locally approximatingthe solutions as a high-
degreepolynomialinsideeachcell.
In this work, computations are performed on the SD7003 airfoil at different Reynolds numbers. The SD7003
airfoilwasselectedduetotheexistenceofexperimentalandcomputationaldataavailableforcomparison.Resultsare
comparedtoexperimentaldataobtainedbyRadespiel2inawatertunnelattheTechnicalUniversityofBraunschweig
(TU-BS)andbyOletal.14 inawaterchannelattheAirForceResearchLabsHorizontalFree-SurfaceWaterTunnel
(HFWT).Thefreestreamturbulenceintensitieswere0.08%and0.1%respectively.Furthermore,resultsarecompared
to the computationaldata obtained by Galbraith and Visbal7 and Uranga et al.6, who both used an ILES approach
coupledwithahigh-orderspatialdiscretization.GalbraithandVisbaluseda6thorderaccuratecompactschemebased
onthepentadiagonalsystemofLele15, with a10th orderlow-passfilter. Thehigh-orderlow-passfilterwasusedto
stabilizethehighorderschemeandincorporatetheadditionaldissipationassociatedwiththeunder-resolvededdies.
Urangaetal.6utilizedahigh-orderDiscontinuousGalerkinsolver13.
This paper starts by briefly presenting the unfiltered Navier-Stokes equations solved in this work, followed by
a brief introduction of the spectral difference method. Then, the computational methodology used is summarized.
Finally, results are presented for flow over the SD7003 airfoil at various Reynolds number. Comparisons with ex-
perimentalandcomputationalresultsaremadeinordertoassessthefeasibilityofthecurrentapproachtoaccurately
predicttheperformanceofairfoilsandwingsinthetransitionalflowregime.
II. Governing Equations
The computations presented in this work are performed by the newly developed flow solver SD3D, a high-
order accurate, unstructured, parallel solver for the Navier-Stokes equations. The SD3D code solves the unsteady,
three-dimensional,compressible,unfilteredNavier-Stokesequationsusingthespectraldifferencemethod. Thethree-
dimensionalunfilteredNavier-Stokesequationscanbeexpressedasfollows:
∂Q ∂F ∂G ∂H ∂F ∂G ∂H
+ I + I + I = V + V + V (1)
∂t ∂x ∂y ∂z ∂x ∂y ∂z
wherethestatevectorQ,inviscidfluxvectorsF ,G andH ,alongwiththeviscousfluxvectorsF ,G andH are
I I I V V V
describedrespectivelyby
ρ ρu ρv ρw
ρu ρu2+p  ρuv   ρuw 
Q= ρv F = ρuv ,G = ρv2+p ,H = ρuw (2)
ρw I  ρuw  I  ρvw  I ρw2+p
       
ρe ρeu+p ρev+p ρew+p
       
3of17
AmericanInstituteofAeronauticsandAstronautics
0 0 0
 σxx   σyx   σzx 
F = σ ,G = σ ,H = σ . (3)
V  xy  V  yy  V  zy 
 σxz   σyz   σzz 
uiσi1−qx uiσi2−qy uiσi3−qz
     
Inthesedefinitions,ρisthedensity,u isthevelocitycomponentinthex directionandeisthetotalenergyperunit
i i
mass. Thepressureisdeterminedfromtheequationofstate,
1
p=(γ−1)ρ e− (u2+v2+w2) . (4)
(cid:18) 2 (cid:19)
ForaNewtonianfluid,theviscousstressesare
∂u ∂u 2 ∂u
σ =µ i + j − µδ k (5)
ij ij
(cid:18)∂x ∂x (cid:19) 3 ∂x
j i k
andtheheatfluxesare
∂T
q =−k . (6)
i
∂x
i
Thecoefficientofthermalconductivityandthetemperaturearecomputedas
C µ p
k = p , T = (7)
Pr Rρ
wherePr isthePrandtlnumber,C isthespecificheatatconstantpressureandRisthegasconstant. Forthecases
p
consideredinthispaper,γ =1.4andPr=0.72.
ItshouldbenotedthattheSD3DcodesolvestheunfilteredNavier-Stokesequationspresentedabovewithoutchange
inlaminar,transitional,orturbulentflow. Noadditionalsub-gridstressandheatfluxtermsareadded;theunresolved
smalleddiesareaccountedforbymeansofnumericaldissipation.
III. 3DSpectral Difference Scheme on Hexahedral Grids
In the presentwork, the Navier-Stokesequationsare solved using the high-orderspectral differencemethodfor
spatialdiscretization.TheformulationoftheequationsonhexahedralgridsissimilartotheformulationofLiuetal.12,
whichwillbesummarizedbelowforcompleteness. ConsidertheunsteadycompressibleNavier-Stokesequationsin
conservativeformwrittenas
∂Q ∂F ∂G ∂H
+ + + =0 (8)
∂t ∂x ∂y ∂z
whereF = F −F , G = G −G andH = H −H . Toachieveanefficientimplementation,allelementsin
I V I V I V
thephysicaldomain(x,y,z)aretransformedtoastandardcubicelement,0 ≤ ξ ≤ 1,0 ≤ η ≤ 1,0 ≤ ζ ≤ 1. The
transformationcanbewrittenas
x K xi
y= Mi(ξ,η,ζ)yi (9)
z Xi=1 zi
   
whereKisthenumberofpointsusedtodefinethephysicalelements,(x ,y ,z )aretheCartesiancoordinatesatthose
i i i
points,andM (ξ,η,ζ)aretheshapefunctions. Thegoverningequationsinthephysicaldomainarethentransferred
i
intothecomputationaldomain,andtheytaketheform
∂Q˜ ∂F˜ ∂G˜ ∂H˜
+ + + =0 (10)
∂t ∂ξ ∂η ∂ζ
whereQ˜ =|J|Qand
F˜ ξ ξ ξ F
x y z
G˜=|J|ηx ηy ηzG (11)
H˜ ζx ζy ζz H
    
4of17
AmericanInstituteofAeronauticsandAstronautics
TheJacobianmatrixJ isgivenby
x x x
∂(x,y,z) ξ η ζ
J = =yξ yη yζ. (12)
∂(ξ,η,ζ)
z z z
ξ η ζ
 
Inthestandardelement,twosetsofpointsaredefined,namelythesolutionpointsandthefluxpoints,asillustratedin
Figure2fora2Delement. Inordertoconstructadegree(N −1)polynomialineachcoordinatedirection,solutionat
Figure2.Positionofsolution(circles)andflux(squares)pointsonthestandard2Delementfor3rdorderSD
N pointsarerequired.Thesolutionpointsin1DarechosentobetheChebyshevpointsdefinedby
1 2s−1
X = 1−cos π ,s=1,2,...,N. (13)
s
2(cid:20) (cid:18) 2N (cid:19)(cid:21)
The flux points were selected to be the Legendre-Gauss quadrature points plus the two end points 0 and 1. The
Legendre-Gaussquadraturepointsaretherootsoftheequation
2n−1 n−1
Pn(ξ)= (2ξ−1)Pn−1(ξ)− Pn−2(ξ)=0 (14)
n n
wherePn(ξ) is the Legendrepolynomialof ordern, P−1(ξ) = 0 and P0(ξ) = 1. In a recentwork, Jameson16
utilizedafluxreconstruction17formulationtoprovethattheSDmethodisstableinanormofSobolevtype,provided
thattheinteriorfluxcollocationpointsareplacedattheLegendre-Gaussquadraturepoints. UsingthesolutionsatN
solutionpoints,adegree(N −1)polynomialcanbebuiltusingthefollowingLagrangebasisdefinedas
N
X −X
h (X)= s . (15)
i
(cid:18)X −X (cid:19)
s=Y1,s6=i i s
Similarly,usingthefluxesat(N +1)fluxpoints,adegreeN polynomialcanbebuiltforthefluxusingtheLagrange
basis
N X −X
s+1/2
l (X)= (16)
i+1/2 s=Y0,s6=i(cid:18)Xi+1/2−Xs+1/2(cid:19)
Thereconstructedsolutionfortheconservedvariablesinthestandardelementisthetensorproductofthethreeone-
dimensionalpolynomials,
N N N Q˜
Q(ξ,η,ζ)= i,j h (ξ)h (η)h (ζ) (17)
i j k
|J |
Xk=1Xj=1Xi=1 i,j
Similarly,thereconstructedfluxpolynomialstakethefollowingform:
N N N
F˜(ξ,η,ζ)= F˜ ·l (ξ)·h (η)·h (ζ) (18)
i+1/2,j,k i+1/2 j k
Xk=1Xj=1Xi=0
5of17
AmericanInstituteofAeronauticsandAstronautics
N N N
G˜(ξ,η,ζ)= G˜ ·h (ξ)·l (η)·h (ζ) (19)
i,j+1/2,k i j+1/2 k
kX=1Xj=0Xi=1
N N N
H˜(ξ,η,ζ)= H˜ ·h (ξ)·h (η)·l (ζ) (20)
i,j,k+1/2 i j k+1/2
kX=0Xj=1Xi=1
Thereconstructedfluxesareonlyelement-wisecontinuous,butdiscontinuousacrosscellinterfaces. Fortheinviscid
flux,aRiemannsolverisemployedtocomputeacommonfluxatcellinterfacestoensureconservationandstability.
Inthecurrentimplementation,theRusanovsolver18 wasused. TheRusanovschemecomputesthecommonnormal
flux,F∗,using
n
1
F∗ = F (Q+)+F (Q−)−α(Q+−Q−) (21)
n 2 n n
(cid:2) (cid:3)
whereQ+ andQ− representsthesolutiononbothsidesofthesharededgefluxpoint,andα∝|u |+c,whereu is
n n
thevelocitynormaltotheedgeandcisthelocalspeedofsound. Theviscousfluxisafunctionofboththeconserved
variables and their gradients, therefore, the solution gradients have to be calculated at the flux points. The average
approachdescribedinreference12isusedtocomputetheviscousfluxes.
IV. Computational Methodology
The computational mesh used for all test cases is a C-grid generated using the flo103 built-in mesh generator,
whichisextrudedinthespanwisedirection. Itcontains128x24x16cellsintheξ,η,ζ directions,foratotalof49152
cells. Tensorproductsofone-dimensionalsecondandthirdorderpolynomialswereused, to obtainthirdandfourth
orderaccuratesolutionsinspacerespectively. ThenumberofsolutionpointspercellisN3,whereN istheorderof
accuracyoftheschemeandtherefore,thethirdandfourthorderaccuratecomputationshaverespectively1.3million
and3.1milliondegrees-of-freedom. Thewingspantochordratioissetto0.2,thesamevalueusedbyGalbraithand
Visbal7. In their work, the effect of changingthe span to chord ratio was studied and it was concludedthat a span
tochordratioof0.2wouldensurethattheflowisnotartificiallyconstrained,whichwouldpreventthree-dimensional
flowstructurestodevelop.Thefarfieldboundaryislocated30chordsawayupstreamanddownstream,and10chords
awayaboveandbelow,toreduceitsinfluenceonthesolutionnearthe airfoil. TheincomingMachnumberissetto
0.1. Periodicboundaryconditionsareusedalongthespanwisedirectionandano-slip,adiabaticwallconditionwas
usedonthesurfaceoftheairfoil.
(a) (b)
Figure3.C-gridwith49152hexahedralelementsgeneratedwithflo103
In this work, a low-storage three stage, third order TVD Runga-Kutta scheme19 was used as the time-stepping
scheme. ParallelexecutionusingMPIisachievedbypartitioningtheunstructuredmeshusingthegraphpartitioning
softwareMETIS.AnexampleofthepartitioningobtainedforthemeshusedisshowninFigure4.
6of17
AmericanInstituteofAeronauticsandAstronautics
Figure4. PartitionsobtainedfromgraphpartitioningsoftwareMETIS
V. Results
Third and fourth order accurate solutions were obtained using the SD3D solver for the flow over the SD7003
airfoilatanangleofattackof4◦ andReynoldsnumbersof10000and60000. In ordertoreducethecomputational
costassociatedwithinitialtransientsfromanauniforminitialsolution,two-dimensionalsolutionswereusedforthe
initialconditiononthethree-dimensionalmesh.Allmeanvalueswerecalculatedbyaveragingthespanwiseaveraged
timeaccuratesolutionoveranon-dimensionaltimeintervalof40.Anon-dimensionaltimestepof5e−5wasusedand
solutionwasrecordedatevery20stepsforthecomputationofthestatisticalquantities. Furthermore,toeliminatethe
effectofinitialtransients,thetimeaveragingprocesswasinitiatedonlyafterconvergenceoftheaveragedquantities.
A. EffectofSpatialDimensions
1. ReynoldsNumber=10000,α=4◦
In orderto investigatethe necessity of three-dimensionalcomputationsat thisReynoldsnumber,a two-dimensional
solution was obtained by using a computational domain with one cell in the spanwise direction. Figures 5 and 6
show the distribution of the mean coefficient of pressure, Cp, and the mean coefficient of friction, C , along the
f
airfoilforthe2Dand3Dcomputations. Theorderofaccuracyofthescheme,N,wassetto3. Bothfiguresindicate
goodagreementbetweenthe two-dimensionalandthree-dimensionalcalculationssuggestingthattheflow is mostly
two-dimensional, at least over the airfoil. The two-dimensionality of the flow is also seen in figures 7(a) and 7(b)
whichshowscontoursofReynoldsstressforthe2Dand3Dsimulationsrespectively,andfigures8(a)and8(b)which
presentscontourofspanwisevorticity. Iso-surfacesoftheQ-criterionwereusedonthe3Dresultsinordertoidentify
three-dimensionalflowstructures.TheuseoftheQ-criterion,initiallyproposedbyDubeifandDelcayre20,providesa
meanofvisualizingvortexcoresandidentifyturbulentstructures.Itcanbecalculatedfrom
1
Q= (Ω Ω −S S ) (22)
ij ij ij ij
2
whereΩ andS aretheanti-symmetricandsymmetricpartofthevelocitygradienttensor.
ij ij
1 ∂u ∂u 1 ∂u ∂u
Ω = i − j andS = i + j (23)
ij ij
2(cid:18)∂x ∂x (cid:19) 2(cid:18)∂x ∂x (cid:19)
j i j i
TheQ-criterionmeasuresthebalancebetweenrateofvorticityΩ2 =Ω Ω andtherateofstrainS2 =S S .Hence,
ij ij ij ij
regionsofpositiveQ-criterioncorrespondstoregionsoftheflowdominatedbyvorticity,asinthecoreofavortexfor
example. Figure9showsiso-surfacesoftheQ-criterionovertheSD7003airfoilatRe = 10000,α= 4◦. Thefigure
showsnospanwisevariationoverthesurfaceoftheairfoilandconfirmsthattheshearlayerdoesnottransitionover
theairfoil
2. ReynoldsNumber=60000,α=4◦
At a Reynoldsnumberof 60000,a two-dimensionalsolutionwas also obtainedto evaluate the necessity of a three-
dimensionalcomputation.Acomparisonbetweenthetwo-dimensionalandthree-dimensionalmeansurfacepressures
7of17
AmericanInstituteofAeronauticsandAstronautics
andmeanskinfrictioncoefficientsarepresentedinfigures10and11,forafourthorderaccuratesolutioninspace.Fig-
ures12(a)and12(b)showcontoursoftheReynoldsstressforbothcases. TheCpplotofthetwo-dimensionalsolution
agreesreasonablywellwith thethree-dimensionalsolution. However,the coefficientof frictiondistributionandthe
Reynoldsstresscontoursaresignificantlydifferent. Amongotherdiscrepancies,theReynoldsstressreachesamuch
higher value for the two-dimensional solution. Futhermore, the skin friction coefficient from the two-dimensional
solution doesnotrise to the same levelof the three-dimensionalsolution downstreamof the reattachmentpoint. In
orderto furtherillustrate thedifferencesbetweenthe two solutions, instantaneouscontoursofvorticityontheplane
z/c = 0.1areshowninfigures13(a)and13(b)forthe2Dand3Dsimulations,respectively. Inthetwo-dimensional
solution,theshearlayerrollsupintoacoherentvortexthatmaintainsitsshapeasittravelsdownstream.Inthethree-
dimensionalsolutionhowever,thecoherentvortexbreaksdownduetotransitiontoturbulentflow. Figure14shows
iso-surfacesoftheQ-criterionovertheSD7003airfoilatRe = 60000,α= 4◦ forthethree-dimensionalsimulation.
Three-dimensionalvorticalstructuresareeasilyidentified,clearlyindicatingthethree-dimensionalnatureoftheflow
atthisReynoldsnumber.
B. ComparisonwithPreviouslyPublishedComputationalandExperimentalResults
Inthissection,resultsobtainedwiththeSD3Dsolverarecomparedtopreviouslypublishednumericalandexperimen-
talresults.ComputationaldataobtainedbyUrangaetal.6andGalbraithandVisbal7areusedforcomparison.Uranga
etal. used a DiscontinuousGalerkinschemewith a fourthorderaccuratespatialdiscretizationona gridcontaining
52,800tetrahedralelements, fora totalof1milliondegreesoffreedom. GalbraithandVisbalusedan oversetmesh
containing5.7millionspoints,with69%ofthepointslocatedonthesuctionsideoftheairfoil. Thespatialdiscretiza-
tionschemeusedintheirworkisa6th ordercompactfinitedifferenceschemewithahighlydiscrimating10th order
filter,thatwaschosentofilteroutonlytheunder-solvedhighwavenumbers.Thepresentnumericalsolutionswillalso
becomparedtohigh-resolutionvelocityandReynoldsstressexperimentalmeasurementsobtainedbyRadespiel2and
Oletal.14fortheflowovertheSD7003airfoilatRe=60000.
1. ReynoldsNumber=10000,α=4◦
Surface pressure and skin friction coefficients obtained with the spectral difference method with a third and fourth
orderaccurateschemearecomparedtoresultsobtainedbyUrangaetal. andGalbraithandVisbalinfigures15and
16. Itcanbeseen thatallmethodsproducenearlyidenticalresults, exceptforsmalldiscrepanciesinthe coefficient
offrictionnearthetrailingedgeoftheairfoil. Giventhetwo-dimensionalnatureoftheflowatthisReynoldsnumber,
theseresultsareexpected. AtthisReynoldsnumber,thethirdorderaccuratesolutionseemstoprovidethenecessary
mesh resolution to capture all features of the flow since the third and fourth order accurate SD schemes produce
nearly identical results. Spanwise vorticity and Reynolds stress obtained with N = 4 are compared to the results
obtainedbyGalbraithandVisbalinfigures17and18.Again,resultsareingoodagreement.Asshownbycontoursof
Reynoldsstress, theflowremainslaminarovertheairfoil. Turbulenttransitionlocationisdeterminedinaccordance
withRadespiel2andissaidtooccurwhentheReynoldsshearstressreachesavalueof0.1%andexhibitsaclearvisible
rise. Table1showsasummaryoftheresultsobtainedforthecaseRe=10000,α = 4◦. Overall,resultsobtainedwith
theSD3Dcodeagreequitewellwiththeothercomputationalexperiments.
DataSet Separation Transition Reattachment Mean Mean
x /C x /C x /C C C
sep tr r L D
GalbraithandVisbal7(ILES) 0.36 - 0.98 0.36 0.047
Urangaetal6(ILES) 0.34 - - 0.38 0.0504
PresentILES 0.362 - 0.995 0.372 0.0492
Table1.MeasuredandComputedpropertiesofflowoverSD7003atRe=10000,α=4◦
2. ReynoldsNumber=60000,α=4◦
At a Reynolds number of 60000, transition takes place across a laminar separation bubble (LSB) over the airfoil.
Thirdandfourthorderaccuratesolutionswereobtainedusingthespectraldifferencemethodandresultsarecompared
withpreviouslypublishednumericalandexperimentalresults. Averagepressurecoefficientandskinfrictionfriction
8of17
AmericanInstituteofAeronauticsandAstronautics
coefficient distributions on the airfoil are compared with numerical results obtained by Uranga6 and Galbraith and
Visbal7 in figures 19 and 20. Increasing the order of accuracy of the SD scheme yields Cp and C distributions
f
thatapproachthoseobtainedbyGalbraithandVisbal7 whouseda6thorderaccurateschemeonamuchfinermesh.
AlthoughthepressureandskinfrictiongradientsinthetransitionregionarenotassharpasthoseobtainedbyGalbraith
andVisbal,theyoccuratthesamelocationandgoodoverallagreementisfound. Pressureandskinfrictiongradients
obtainedwiththeSD3DcodeinthetransitionregionaresteeperthantheoneobtainedbyUrangaetal6. Thisismost
likelyduetothe finergridresolutionusedwiththe SD3D code. ComputedcontoursofReynoldsstressesobtained
from the TU-BS and HFWT experimentalmeasurementsand the numericalsimulation by Galbraith and Visbal are
showninfigure21.ThoseresultsarecomparedwiththirdandfourthorderaccuratesolutionsobtainedwiththeSD3D
solverin figure22. Similarly, figures23and24 compareaveragespanwise vorticitycontours. ComputedReynolds
stressesobtainedwithN=3andN=4agreewellwithcomputationalresultsfromGalbraithandVisbalandexperimental
resultsintermsofshape,magnitudeandextent.Similarly,generalshapeandextentofthecomputedshearlayershown
withcontoursofspanwisevorticityinfigure24differslittlefromnumericalandexperimentsmeasurements. Table2
comparesseparation,transitionandreattachmentat4◦ angleofattackmeasuredfromthetwoexperimentalfacilities
alongwithsimulationsbyYuanetal.4,Lianetal.3,Uranga,GalbraithandVisbal,andthepresentILEScomputation.
The separation and reattachmentlocationswere recordedby lookingat the streamlines of the mean velocityprofile
showninfigure25. Again,turbulenttransitionisassumedtooccurwhentheReynoldsstressreachesavalueof0.1%
and exhibits a clear visible rise. In the RANS-eN calculations by Yuan et al. and Lian et al., a critical N-factor
of 8 was used, based on the empiricalrelationship between freestream turbulence intensity and critical N-factor by
Mack21. Results obtainedwith the spectraldifferencemethod by meansof a ILES are in excellentagreementwith
computationalresultsofGalbraithandVisbal7.ThetransitionandreattachmentlocationpredictedbythecurrentILES
alsoagreewellwiththeTU-BSmeasurements.Thus,eventhoughthegridresoltionisclearlytoocoarsetocaptureall
scalesoftheflow,andthedissipationfromthenumericalschemeisunabletodiscriminatebetweenresolvedandunder-
resolvedscales,theproposedmethodappearstoaccuratelypredictthelaminarseparationandsubsequenttransitionto
turbulentflow.
DataSet Freestream Separation Transition Reattachment
Turbulence x x /c x /c
sep tr r
TU-BS2 0.08% 0.30 0.53 0.64
HFWT14 0.1% 0.18 0.47 0.58
YuanSGS-LES4 0.1%,N=8 0.21 0.49 0.60
YuanRANS-eN4 0.1%,N=8 0.21 0.49 0.58
LianRANS-eN4 0.1%,N=8 0.21 0.48 -
Visbal(ILES)7 0 0.23 0.55 0.65
Uranga(ILES)6 0 0.23 0.51 0.60
PresentILES,N=3 0 0.23 0.52 0.65
PresentILES,N=4 0 0.23 0.52 0.65
Table2.MeasuredandComputedpropertiesofflowoverSD7003atRe=60000,α=4◦
VI. Conclusion
Thispaperpresentscomputationalresultsforthepredictionoftheformationofalaminarseparationbubble(LSB)
and its subsequent burst to turbulentflow on the SD7003 airfoil. Solutions were obtained by means of an Implicit
LargeEddySimulation(ILES)withanewlydevelopped,parallelNavier-Stokessolverthatusesthespectraldifference
methodforspatialdiscretization.Withafairlycoarsegrid,computedseparation,transitionandreattachmentlocations
were in good agreementwith previouslypublished experimentaland computationalresults. The use of the spectral
differencemethodwhichprovideslow numericaldispersioncombinedwithan artificiallydissipativenumericalflux
seemstopermittheaccuratesimulationoftransitionalflowovertheSD7003airfoilattheReynoldsnumberandangle
of attack considered. As pointedby Visbal and Rizetta5, the artificial dissipation fromthe numericalflux is unable
to discriminatebetweenresolvedandunder-resolvedscalesbutdespitethis, the currentILESschemeachievedvery
goodresultsforthetestcasesconsidered. Futureworkisrequiredtoevaluatethefeasibilityofthecurrentapproach
forotherwelldocumentedtestcasessuchasanisotropicdecayingturbulentfloworaturbulentchannelflow.
9of17
AmericanInstituteofAeronauticsandAstronautics
Acknowledgments
TheauthorswouldliketoacknowledgethesupportforthisworkprovidedbytheStanfordGraduateFellowships
program, the NaturalScience and EngineeringResearch Council(NSERC) of Canada and the FondsQue´be´cois de
la RecherchesurlaNatureetlesTechnologies(FQRNT).Thisresearchwasalso supportedbytheNationalScience
FoundationthroughtheTeraGridresourcesprovidedbytheTexasAdvancedComputingCenter.
References
1Horton,H.,LaminarSeparationBubblesinTwoandThreeDimensionalIncompressibleFlow,Ph.D.thesis,UniversityofLondon,1968.
2R.Radespiel, J.W.andScholz,U.,“NumericalandExperimentalFlowAnalysisofMovingAirfoilswithLaminarSeparationBubbles,”
Aiaapaper2006-501,2006.
3Lian,Y.andShyy,W.,“Laminar-TurbulentTransitionofaLowReynoldsNumberRigidorFlexibleAirfoil,”Aiaapaper2006-3051,2006.
4W.Yuan,M.Khalid,J.W.U.S.andRadespiel,R.,“AnInvestigationofLow-Reynolds-numberFlowspastAirfoils,”Aiaapaper2005-4607,
2005.
5Visbal,M.andd.P.Rizetta,“Large-EddySimulationonCurvilinearGridsUsingCompactDifferencencingandFilteringSchemes,”Journal
ofFluidsEngineering,Vol.124,2002,pp.836–847.
6A.Uranga,P.Persson,M.D.andPeraire,J.,“ImplicitLargeEddySimulationofTransitionalFlowsOverAirfoilsandWings,”Tech.rep.,
16.
7Galbraith, M. andVisbal, M., “Implicit Large EddySimulation ofLow Reynolds NumberFlow Pastthe SD7003 Airfoil,” Aiaa paper
2008-225,46thAerospaceScienceMeetingandExhibit,Reno,Nv,2008.
8Spyropoulos, E.andBlaisdell, G.,“EvaluationoftheDynamicModelforSimulationsofCompressibleDecayingIsotropicTurbulence,”
AIAAJournal,Vol.34,No.5,1996,pp.990–998.
9Kopriva,D.A.andKolias,J.H.,“AConservativeStaggered-GridChebyshevMultidomainMethodforCompressibleFlows,”Journalof
ComputationalPhysics,Vol.125,1996,pp.244–261.
10Liu,Y.,Vinokur,M.,andWang,Z.J.,“SpectraldifferencemethodforunstructuredgridsI:Basicformulation,”JournalofComputational
Physics,Vol.216,2006,pp.780–801.
11Wang,Z.,Liu,Y.,May,G.,andJameson,A.,“SpectralDifferenceMethodforUnstructuredGridsII:ExtensiontotheEulerEquations,”
JournalofScientificComputing,Vol.32,2007,pp.45–71.
12Y.Sun,Z.andLiu,Y.,“High-OrderMultidomainSpectralDifferenceMethodfortheNavier-StokesEquationsonUnstructuredHexahedral
Grids,”CommunicationsinComputationalPhysics,Vol.2,No.2,2007,pp.310–333.
13F.BassiandRebay, S.,“High-orderaccurate discontinuous finiteelementsolutionofthe2Deulerequations,” JournalofComputational
Physics,Vol.138,1997,pp.251–285.
14M.V.Ol,B.R.McAuliffe,E.S.H.U.andKahler,C.,“ComparisonofLaminarSeparationbubbleMeasurementsonaLow-ReynoldsNumber
AirfoilinThreeFacilities,”Aiaapaper2005-5149,2005.
15Lele,S.,“CompactFiniteDifferenceSchemeswithSpectral-LikeResolution,”JournalofComputationalPhysics,Vol.103,1992,pp.16–42.
16Jameson,A.,“Aproofofthestabilityofthespectraldifferencemethodforallordersofaccuracy,”J.Sci.Comput.,2010.
17Huynh,H.,“AFluxReconstructionApproachtoHigh-OrderSchemesIncludingtheDiscontinuousGalerkinMethods,”Aiaapaper2007-
4079,2007.
18Rusanov, V.V., “Calculation ofinteraction ofnon-steady shockwaves with obstacles,” Journal ofComputational MathPhysics USSR,
Vol.1,1961,pp.261–279.
19Gottlieb,S.andShu,C.,“TotalvariationdiminishingRunga-Kuttaschemes,”MathematicsofComputation,Vol.67,No.221,1998,pp.73–
85.
20Dubeif,Y.andDelcayre,F.,“OnCoherent-vortexIdentificationinTurbulence,”JournalofTurbulence,Vol.1,No.11,2000,pp.1–22.
21Mack,L.,“TransitionandLaminarInstability,”Jetpropulsionlaboratorypublication77-15,Pasadena,CA,1977.
10of17
AmericanInstituteofAeronauticsandAstronautics