Grid Convergence for Turbulent Flows (Invited) ∗ BorisDiskin, NationalInstituteofAerospace,Hampton,VA23666 † ‡ JamesL.Thomas, ChristopherL.Rumsey, NASALangleyResearchCenter,Hampton,VA23681 § AxelSchwo¨ppe DLR,Braunschweig,Germany Adetailedgridconvergencestudyhasbeenconductedtoestablishaccuratereferencesolutionscorrespond- ingtotheone-equationlineareddy-viscositySpalart-Allmarasturbulencemodelfortwodimensionalturbu- lentflowsaroundtheNACA0012airfoilandaflatplate. Thestudyinvolvedthreewidelyusedcodes,CFL3D (NASA),FUN3D(NASA),andTAU(DLR),andfamiliesofuniformlyrefinedstructuredgridsthatdifferinthe griddensitypatterns. Solutionscomputedbydifferentcodesondifferentgridfamiliesappeartoconvergeto thesamecontinuouslimit,butexhibitdifferentconvergencecharacteristics.Thegridresolutioninthevicinity of geometric singularities, such as a sharp trailing edge, is found to be the major factor affecting accuracy andconvergenceofdiscretesolutions,moreprominentthandifferencesindiscretizationschemesand/orgrid elements. TheresultsreportedfortheserelativelysimpleturbulentflowsdemonstratethatCFL3D,FUN3D, andTAUsolutionsareveryaccurateonthefinestgridsusedinthestudy,buteventhosegridsarenotsufficient toconclusivelyestablishanasymptoticconvergenceorder. ∗NIAResearchFellowandResearchAssociateProfessor, MAEDepartment, UniversityofVirginia, Charlottesville, VA,AssociateFellow AIAA. †DistinguishedResearchAssociate,ComputationalAeroSciencesBranch,FellowAIAA. ‡SeniorResearchScientist,ComputationalAeroSciencesBranch,FellowAIAA. §ResearchEngineer,InstituteofAerodynamicsandFlowTechnology. 1of50 AmericanInstituteofAeronauticsandAstronautics Nomenclature AR Aspectratioatleadingedge k Meshgradationcoefficient LE C Skinfrictioncoefficientcomponent rˆ Distancefromreentrantcornertip f,x C Pressurecoefficient p Convergenceorder p C Liftcoefficient q ExactsolutionoftheLaplaceequation L C Pitchingmomentcoefficient s Sizeofreentrant-cornercell M C Totaldragcoefficient u,w Velocitycomponents D C Pressuredragcoefficient xˆ,zˆ Coordinatesrelativetoreentrantcornertip Dp C Viscousdragcoefficient x,y,z Orthogonalcoordinatesinstreamwise, Dv M∞ FreestreamMachnumber spanwise,andnormaldirections N Numberofdegreesoffreedom x ,z Reentrantcorner/trailing-edgeposition TE TE P Referencepressure z+ Nondimensionalboundarylayerspacing ref P Totalpressure α Angleofattack t Pr Prandtlnumber αˆ Inverserelativeexteriorangle Pr TurbulentPrandtlnumber β Stretchingfactor t Re Reynoldsnumber (cid:4) Discretizationerror T Referencetemperature θ Polarcoordinateangle ref T Totaltemperature κ MUSCLschemeparameter t a Referencespeedofsound μ Referenceeddyviscosity ref ref c Chordlength μ Eddyviscosity t havg Averagedmeshspacing ν∞ Freestreamkinematicviscosity h,h Characteristicmeshsize ξ Mappingcoordinate eff h Maximummeshsize ω Externalangleforreentrantcorner max I. Introduction Withever-increasingcomputingpowerandrecentadvancementsinsolvertechnology,turbulentflowsareroutinely simulatedonhigh-densitygridswithmanymillionsofdegreesoffreedom.WhileaccurateandreliableComputational Fluid Dynamics (CFD) analysis of attached turbulent flows is now possible, accuracy and robustness of separated turbulent flow simulations in complex geometries are still not adequate. Errors in such simulations are typically attributed to three sources. (1) The modeling error is due to approximations in the continuous formulation (e.g., in differential equations describing turbulent flows, or in geometry definitions, or in boundary conditions) and is defined as the difference between the exact continuous solution of the model formulation and the real flow. (2) The discretizationerrorisduetoapproximationsindiscretizingthecontinuousformulationonaspecificgridandisdefined asthedifferencebetweentheexactdiscreteandcontinuoussolutions. (3)Theiterative(algebraic)errorisduetoan imperfectiterativesolutionprocessforthediscreteformulationandisdefinedasthedifferencebetweentheexactand approximatesolutionstothediscreteformulation. A Turbulence Modeling Resource (TMR) website1,2 has been recently established at NASA Langley Research Centertodescribe,standardize,verify,andvalidateformulationsofcommonturbulencemodels. Thepurposeofthis website is to avoid ambiguity associated with specific implementations of turbulence models in CFD codes. The turbulentflowsconsideredinthispaperaremodeledbyReynolds-averagedNavier-Stokes(RANS)equationsandthe one-equationlineareddy-viscositySpalart-Allmaras(SA)turbulencemodel.3 Theformulationdetailsareavailableat theTMRwebsite. Thispaperreportsonanattempttoeliminate(oratleastminimize)thediscretizationanditerativeerrorsbycon- ductinganextensivegridconvergencestudyforrelativelysimplebenchmarkturbulentflows. Currentguidelinesfor gridconvergencestudies4,5emphasizeaparametricsimilarityofgridsformingafamilyandanasymptoticconvergence order,whichisexpectedtobeobservedonthree-to-fourfinegridsinafamily. Themeshresolutionrequiredforestab- lishingaconvergenceorderissoughtthroughauniformgridrefinement. Forstructuredgrids,afamilyofuniformly refinedgridsistypicallyderivedrecursively, startingfromthefinestgrid. Eachcoarsergridinthefamilyisderived from the preceding finer grid by removing every-other grid plane/line in each dimension. This uniform-refinement approach has rigorous mathematical foundations. However, it is also expensive because it lacks the flexibility of a localrefinement,whichisthebasisforeffectivegridadaptationmethods. Ideally,solutionsobtainedongridsinafamilywouldmonotonicallyconvergeingridrefinementtothecontinuous 2of50 AmericanInstituteofAeronauticsandAstronautics solution with the same asymptotic convergence order for any solution quantity. Then, the entire solution can be extrapolatedtothelimitoftheinfinitegridrefinement. A single asymptotic convergence order characterizing the entire solution is typically observed only for model problems with smooth continuous solutions. In practical computations, an asymptotic convergence order remains elusive. ComputationspresentedattherecentDragPredictionWorkshop(DPW-V)6,7 illustratethedifficultyinreal- izingasymptoticconvergenceorderforturbulentflows. Arelativelycloseagreementwasobservedbetweendifferent computationsinpredictingthedragonathree-dimensional(3D)wing-bodyconfiguration, butthegridconvergence propertieswereverydifferentbetweensolutionscomputedwithdifferentcodesorevenbetweensame-codesolutions computed on different grid families. In anticipation of second-order convergence, the values of drag were plotted versusN−23,whereN isthenumberofdegreesoffreedom,butsecond-orderconvergencewasnotobserved. The study reported in this paper aims at computing highly accurate reference solutions for some benchmark turbulent-flowcasesandatprovidingsomeguidanceonaccuracyvariationforgridfamilieswithdifferentmeshden- sity patterns. The main grid-convergence test case considered is a turbulent flow around the 2D NACA 0012 airfoil at10◦ angleofattack. TheRANSequationsaresolvedonuniformly-refined,structured,high-densitygridsbythree well-establishedCFDcodes: CFL3D(NASA),FUN3D(NASA),andTAU(DLR),whichusedifferentdiscretization anditerationschemes. Advancedturbulent-flowsolvertechnologiesrecentlyimplementedintwoofthesecodes8–10 providemeansforminimizingeffectsofiterativeerrors. FUN3DandTAUconvergeallresiduals,includingtheresid- ualoftheSAturbulence-modelequation,tolevelscomparablewiththemachinezero. CFL3Dconvergesthedensity residualtothelevelof10−13 andtheSAmodelresidualtothelevelof10−7;thecorrespondingaerodynamicforces convergetoatleastfivesignificantdigits;andthepitchingmomentconvergestoatleastfoursignificantdigits. The study began with an attempt to characterize the grid-refined solutions by using grids offered in the “Turbu- lence Model Validation Cases and Grids” section of the TMR website. The grid-refinement study was conducted using FUN3D solutions on the family of grids then available9 — the finest grid had about 1M degrees of freedom. Convergenceoftheliftandpitchingmomentcoefficientsdidnotexhibitanyclearorderproperty. Theliftvaluesin- creasedwithgridrefinementandthendecreasedonthefinestgrid. Thepitchingmomentcoefficientswerecontinually increasingwithgridrefinement,butshowednoconsistentorder. Adetailedinspectionofthesolutionsonthesurface oftheairfoilrevealederraticconvergenceofthepressurecoefficientsnearthetrailingedge. ThisobservationmotivatedthecurrentgridconvergencestudythatinvolvessolutionsobtainedwiththethreeCFD codesonthreeexpandedgridfamilies. Thefinestgridsineachfamilyhaveabout15Mdegreesoffreedom. Thegrid families differ in the trailing-edge resolution and are now available in the “Cases and Grids for Turbulence Model Numerical Analysis” section of the TMR website. Convergence sensitivities to the trailing-edge resolutions as well as to various discretization aspects, such as grid elements and the order of approximation for the turbulence-model convectionterm,havebeenconsidered. BesidestheNACA0012study,anexistingflatplatetestcasehasbeenextendedandusedtostudygridconvergence. FUN3D solutions have been computed on a set of structured grids. The grids are also available in the “Cases and GridsforTurbulenceModelNumericalAnalysis”sectionoftheTMRwebsite. Methodologically, thepaperfollows thecurrentguidelinesforgridconvergencestudy: familiesofuniformlyrefinedgridsareusedandconvergenceorders oflocalandglobalsolutionquantitiesarereported. Additionally,theconvergencedegradationlong-knownforsolutionsofelliptic(purediffusion)equationsonuni- formlyrefinedgridswithgeometricsingularities11 isrevisited. Ellipticequationsdescribediffusionphenomenaand thusapplydirectlytolow-Reynoldsnumber(i.e.,Stokes)flows. Nearsurfaces,turbulent-flowsolutionsareexpected to be similar to Stokes-flow solutions. For a sharp trailing edge, the pure-diffusion solution exhibits a square-root behaviornearthetrailing-edgesingularityandhasunboundedderivatives. Thediscretizationerrorconvergesonuni- formlyrefinedgridswiththefirstorderintheL1normandwithanorderof0.5intheL∞norm.Aseriesofstructured (non-uniformly) refined grids that have a higher refinement rate near the singularity than in the rest of the domain aredevelopedandshowntorecovertheconvergencerateobtainedforsmoothflowsondomainswithoutsingularities. This grid refinement strategy has not been applied to turbulent-flow computations. However, its success in applica- tiontothepure-diffusionequationprovidesanindicationthatimprovedresolutionneargeometricsingularitieswould improvetheaccuracyofturbulent-flowsimulations. The material in the paper is presented in the following order. First, the CFL3D, FUN3D, and TAU codes used inthestudyaredescribedinSectionII, includingdiscretizationdetailsanditerativeconvergencestrategies. Then, a benchmarkturbulentflowaroundtheNACA0012airfoilisintroducedinSectionIIItogetherwithadescriptionofthe gridfamiliesandthecorrespondingnumericalsolutionsobtainedwiththethreecodes. Adetaileddescriptionofthe solutions and grid convergence is provided for future verification of CFD solvers. A study of solution sensitivity to thevariationofdiscretizationmethodsandgridelementsispresented. SectionIVreportsonagridconvergencestudy 3of50 AmericanInstituteofAeronauticsandAstronautics for a turbulent flow around a finite flat-plate configuration. Finally, concluding remarks are offered. The appendix considers effects of geometric singularity on accuracy for solutions of a pure-diffusion equation and presents some strategiestoovercometheconvergencedegradation. II. CFDCodesUsedintheStudy This section introduces the three well-established large-scale practical CFD codes used in the grid-convergence studyfortheNACA0012airfoil. Thecodesrepresentthestateoftheartinaerodynamiccomputations. A. CFL3D CFL3Disastructured-gridmulti-blockcell-centeredfinite-volumecodewidelyappliedforanalysisofcomplexflows. It has been used in many recent workshops involving complex turbulent flows12–15 and for computing benchmark turbulent-flow solutions at the TMR website. It uses second-order, upwind-biased spatial differencing scheme (a MUSCL scheme16,17 corresponding to κ = 1/3 that allows a third-order accuracy in one dimension (1D)) for the convective and pressure terms, and second-order differencing for the viscous terms; it is globally second-order ac- curate. Roe’sfluxdifference-splittingmethodisusedtoobtainfluxesatthecellfaces. Theoptiontomodelthefull Navier-Stokes meanflow equations is exercised for all cases. In distinction from the other two codes that use the SA-neg scheme,18 CFL3D uses the standard SA one-equation turbulence model for this study. The negative values oftheSpalartturbulencevariablearenotallowed; theminimumvaluesareclippedat10−12. Firstandsecond-order approximationsfortheconvectiontermintheturbulence-modelequationareavailable.Asecond-orderapproximation isusedfortheNACA0012caseongridsofFamilyII(seeSectionIIIforgridfamilydefinitions). Initially,thefirst- order approximation was used for all computations; however, the second-order approximation was found to make a significantdifferenceonFamilyIIgrids(seefiguresontheTMRwebsite). Inthisstudy,thefirst-orderapproximation isusedforcomputationsongridsofFamilyIandFamilyIII.Theturbulence-modeldiffusiontermusesthethin-layer approximation.Theiterationschemeislooselycoupled,i.e.,first,themeanflowequationsareadvancedwiththeeddy- viscosity fixed, then the turbulence-model equation is advanced with the meanflow solution fixed. CFL3D employs localtime-stepscaling,gridsequencing,andmultigridtoaccelerateconvergencetosteadystate. B. FUN3D FUN3Disafinite-volume,node-centered,unstructured-gridRANSsolver,whichisalsowidelyusedforhigh-fidelity analysis and adjoint-based design of complex turbulent flows.15,19–25 FUN3D solves governing flow equations on mixed-elementgrids; theelementsaretetrahedra,pyramids,prisms,andhexahedra. Atmedian-dualcontrol-volume faces,theinviscidfluxesarecomputedusinganapproximateRiemannsolver. Roe’sfluxdifferencesplittingisused in the current study. For second-order accuracy, face values are obtained by a MUSCL scheme, with unweighted least-squaresgradientscomputedatthenodes. Ifgridlinesareavailable,e.g.,withinboundarylayersorinthewake, thereisanoptiontouseadirectionalgradientthatexploitsa1Dline-mappingalongthegridlines. Forthisstudy,the MUSCLschemecoefficientissettoκ = 0.5forthemeanflowequationsandtoκ = 0fortheturbulenceconvection term. Theviscousfluxesusefullviscousstresses. Fortetrahedralmeshes,theviscousfluxesarediscretizedusingthe Green-Gausscell-basedgradients;thisisequivalenttoaGalerkintypeapproximation.Fornon-tetrahedralmeshes,the edge-basedgradientsarecombinedwithGreen-Gaussgradients;thisimprovestheh-ellipticityoftheviscousoperator. The diffusion term in the turbulence model is handled in the same fashion as the meanflow viscous terms. FUN3D usestheSA-negvariantoftheSAturbulencemodel18 thatadmitsnegativevaluesfortheSpalartturbulencevariable. Thisvariantwasdesignedforimprovednumericalbehavior. TheSA-negmodelisidenticaltotheoriginalSAmodel forpositivevaluesoftheSpalartturbulencevariable. FUN3Dusesasecond-orderapproximationfortheconvection termintheturbulence-modelequation. Amultigridsolverisusedtoconvergeresiduals. Therelaxationschemeinthismultigridsolverisahierarchical nonlinearscheme. Ontheinnermostlevelitusesapreconditionerbasedonadefect-correctionmethodanditerateson asimplifiedfirst-orderJacobianwithapseudo-timeterm. Onepreconditioneriterationinvolvesanimplicit-linepass through the portionof the domain where implicitgrid lines are defined, followedby a point-implicit sweep through theentiredomain. Thenumberofpreconditioneriterationsmayvaryfordifferentnonlineariterations. Thisvariable preconditioneriswrappedwithaGeneralizedConjugateResidual(GCR)methodtoformaJacobian-freelinearsolver thatusesFrechetderivativestoapproximatelinearresiduals. Anonlinearcontrollerassessesthecorrectioncomputed by the linear solver. The controller is responsible for the CFL adaptation strategy and for deciding when to update the Jacobian. As a result of this assessment, the suggested correction can be applied fully, partially, or completely 4of50 AmericanInstituteofAeronauticsandAstronautics discarded; the current Jacobian may be updated or reused in the next iteration; and the current CFL number may increase, decrease, or stay the same. In the relaxation scheme, the iterations can be tightly or loosely coupled, i.e., applied to the meanflow and turbulence equations collectively or separately. The multigrid iterations are always coupled in the sense that the meanflow and turbulent equations are solved on coarse grids and the meanflow and turbulence variables use a coarse-grid correction. Initially, the CFL number is ramped over a prescribed number of iterations,butthenitautomaticallychangeswithinprescribedbounds. Thecoarse-gridcorrectionsarealsoassessed bythenonlinearcontrollerandcanbeappliedfully,partially,orcompletelydiscarded. C. TAU TAUisafinite-volumenode-centeredunstructured-gridRANSsolverwidelyusedforabroadrangeofaerodynamic andaero-thermodynamicproblems.26 Itofferscouplinginterfacestootherdisciplineslikestructureandflightmechan- ics to allow for multidisciplinary simulations.27 A full derivative is available for adjoint-based shape optimization. TAUsolvesthe3Dcompressibletime-accurateRANSequationsongridswithmixedelements,includingtetrahedra, pyramids,prisms,andhexahedra. Controlvolumesareconstructedbymedian-dualpartition. Thenumericalscheme isbasedonasecond-order,finite-volumeformulation.Atcontrolvolumefaces,theinviscidfluxesarecomputedusing acentraldifferenceschemewithanaddedmatrix-valuedartificialviscosity.10 Todealwithhighlystretchedmeshes, acellstretchingcoefficientisincludedintothescheme. Thefullviscousfluxesofthemeanflowandturbulenceequa- tionsarediscretizedusinganedge-normalgradientformulationasanaugmentedaverageoftheadjacentGreen-Gauss cell gradients.28 Various turbulence models are available, ranging from eddy viscosity to full differential Reynolds stress models,29 including options for Large-Eddy Simulation (LES) and hybrid RANS/LES. The SA-neg model18 isusedastheturbulencemodelinthisstudy, andtheSAmodelconvectiontermisdiscretizedusingasecond-order approximation. Amultigridsolverbasedonagglomeratedcoarsegridsisusedtoconvergetosteadystate. Thebaselinerelaxation scheme of TAU in this multigrid solver is an implicit Lower-Upper Symmetric Gauss-Seidel (LU-SGS) scheme.30 Recently the LU-SGS scheme was embedded in a hierarchy of smoothers derived from a general implicit Runge- KuttamethodtofurtherimprovereliabilityandefficiencyofthesolutionsalgorithmsofTAU.10 Thesmoothers, e.g. a first-order preconditioned Runge-Kutta or Newton-Krylov generalized minimal residual (GMRES) methods, can be considered as simplified Newton methods. The smoothers differ in Jacobian approximations and in the solution methodsusedforthearisinglinearsystems. III. NACA0012Airfoil AgridconvergencestudyforaturbulentflowaroundtheNACA0012airfoilispresentedinthissection. Thistest casecorrespondstotheNACA0012caseinthe“CasesandGridsforTurbulenceModelNumericalAnalysis”section of the TMR website. The goals of this study are (1) to establish an accurate reference solution that can be used for verification of CFD solutions computed with the SA turbulence model, (2) to evaluate the effects of grid resolution nearasharptrailingedgeonconvergenceofturbulent-flowsolutions,and(3)toassesssensitivityofthesolutionsto variationsofdiscretizationmethodsandgridelements. A. FlowParameters,BoundaryConditions,andDiscretizationDetails Aturbulentessentiallyincompressible(M∞ =0.15)flowaroundtheNACA0012airfoilatα=10◦angleofattackis considered. Forthepurposesofthisstudy,thedefinitionoftheNACA0012airfoilisslightlyalteredfromtheoriginal definition, sothattheairfoilclosesatc = 1withasharptrailingedge. TheexactdefinitionisavailableattheTMR website.1 TheReynoldsnumbercomputedperchordlengthisRe = 6M.Thecomputationaldomainandboundary conditionsaresketchedinFig.1.Thefarfieldboundaryconditionsarebasedoninviscidcharacteristicmethods.Ano- slipadiabaticwallconditionisspecifiedattheairfoilsurface. FUN3DandTAUhaveastrongimplementationofthe wallboundaryconditions;CFL3Dhasaweakimplementationofthewallboundaryconditions. T =540◦ Risthe ref freestreamstatictemperature. Althoughafarfieldpointvortexboundaryconditioncorrection31isrecommendedatthe TMRwebsite, theresultsbelowarepresentedwithoutsuchacorrection. Thissimplificationfacilitatescomparisons withemerginghigh-orderandmeshadaptationcapabilities.32,33 ThefarfieldvalueoftheSpalartturbulencevariable isν˜farfield = 3ν∞. ThePrandtlnumberistakentobeconstantatPr = 0.72, andtheturbulentPrandtlnumberis takentobeconstantatPr =0.9. ThemolecularviscosityiscomputedusingSutherland’sLaw.34 t 5of50 AmericanInstituteofAeronauticsandAstronautics Figure1. Domainandboundaryconditions. Table1. NACA0012:summaryofalong-surfacemeshspacingonthefinest7,169×2,049gridsinthreefamilies. Grids x≈0 x≈1 x≈0.5 (Leading-edge) (Trailing-edge) (Middleofthesurface) FamilyI 0.0000125c 0.0001250c 0.00123c FamilyII 0.0000125c 0.0000125c 0.00155c FamilyIII 0.0000125c 0.0000375c 0.00139c B. Grids (a)Farview. (b)Nearview. Figure2. Computationaldomainanda449×129gridforNACA0012airfoil. Three families of grids are generated with a farfield extent of approximately 500 chord lengths. Figures 2 (a) and(b)showtwoviewsofthe449×129gridofFamilyI.FamilyIgridshavethedensitydistributionsimilartothe distributionusedongridsofthefamilyavailableontheTMRwebsitepriortothisstudy. FamilyIIgridsareclustered 6of50 AmericanInstituteofAeronauticsandAstronautics (a)FamilyI. (b)FamilyII. (c)FamilyIII. Figure3. Nearviewoftrailing-edgegridsofFamiliesI,II,andIII. near the trailing edge and Family III grids are intermediate between the two. See Figs. 3 (a), (b), and (c) for near- trailing-edgeviewsofthe449×129gridsforeachfamily.Aseriesofsevennestedgridsaregeneratedforeachfamily, rangingfromthefinest7,169×2,049tothecoarsest113×33grid. Thegridtopologyisaso-called“C-grid.”Each ofthegridswrapsaroundtheairfoilfromthedownstreamfarfield,aroundthelowersurfacetotheupper,thenbackto thedownstreamfarfield;thegridconnectstoitselfinaone-to-onefashioninthewake. Thereare4097pointsonthe airfoilsurfaceonthefinestgrid,with1537pointsalongthewakefromtheairfoiltrailingedgetotheoutflowboundary. Eachfamily’sfinestgridhastheminimumnormalspacingatthewallof10−7. Thespacingalongtheairfoilsurface isdocumentedinTable1. Theleading-edgespacingisthesameforallfamiliesandcorrespondstotheaspectratioof AR = 125. Thetrailing-edgespacingislargestfortheFamilyIgridsandtentimeslargerthantheleading-edge LE spacing.OnFamilyIIgrids,thetrailing-edgespacingisthesameastheleading-edgespacing.OnFamilyIIIgrids,the trailing-edgespacingisbetweenthatofFamilyIandFamilyIIandthreetimeslargerthanthetrailing-edgespacingof thecorrespondingFamilyIIgrids. Thefamilynameconventionisnotconsistentwiththevariationofthetrailing-edge meshspacing. Theauthorschoosetokeepthesamenamesforgridfamiliesasinthe“CasesandGridsforTurbulence ModelNumericalAnalysis”sectionoftheTMRwebsite. Themeshspacinginthemiddleoftheairfoilsurfacechangesbetweenthefamilies. Thealong-surfacespacingsat x≈0.5are0.00123c,0.00139c,and0.00155cforfamiliesI,III,andII,respectively. Thecorrespondingaspectratios are 12300, 13900, and 15500. The relative increase in the middle-of-chord mesh spacing and aspect ratio between familiesIandIIisapproximately25%. Themiddle-of-chordaspectratiosareapproximatelytwoordersofmagnitude higherthanthoseattheleadingedge. C. GridConvergenceofAerodynamicCoefficients ThissectionreportsonconvergenceofaerodynamiccoefficientsongridsoffamiliesI,II,andIII.Figures4–8compare convergenceofthetotaldrag(CD),pressuredrag(CDp),viscousdrag(CDv),lift(CL),andpitchingmomen(cid:2)t(CM) withrespecttothequarter-chordlocation. Thevariationsareshownversusacharacteristicmeshspacingh= 1/N. FUN3DcomputationsareshownonlyonthefourfinestgridsofFamilyIII.Toaccommodateadetailedscaleforthe C andC coefficients,onlyresultsonthethreefinestgridsineachfamilyareshowninFigs.7–8. L M ConvergenceplotsofdragcoefficientsshowninFigs.4–6aresimilarongridsofdifferentfamilies. Convergence plotsofliftandpitchingmomentdiffersignificantlybetweengridfamilies. Note,however,that,relativelyspeaking, the vertical scale for the lift figures is significantly smaller (showing variations in the fourth significant digit) than verticalscalesforthedragandmomentfigures(showingvariationsinthethirdandfirstsignificantdigits,respectively). Relatively large deviations of the CFL3D lift and moment coefficients from the corresponding FUN3D and TAU coefficientsobservedonFamilyIandFamilyIIIgridsarepartiallyexplainedbyvariationsinthediscretizationscheme for the SA model equation. Recall that CFL3D solutions are computed with the first-order approximation for the convectiontermintheSAmodelequationongridsofFamilyIandFamilyIII.Althoughnotshownhere, resultson theTMRwebsitedemonstratetheeffectoftheSAmodeldiscretizationorderwhenusingFamilyIIgrids. Allaerodynamiccoefficientsarepredictedwithasmallvariationbetweenallthethreecodesonthefinestgridsof allfamilies:thedragvariationbetweenthecodesandgridfamiliesislessthan1count(lessthan1%),theliftvariation 7of50 AmericanInstituteofAeronauticsandAstronautics islessthan0.2%,andthepitchingmomentvariationislessthan10%. Themaximumvariationisobservedbetween thethreecodesontheFamilyIgrid;thecorrespondingvariationontheFamilyIIgridisanorderofmagnitudesmaller: 0.1count(lessthan0.1%)forC ,0.005%forC ,and1.5%forC . D L M ThevariationsbetweenthecodesonfinergridswithinthesamefamilyarealsosmallerontheFamilyIIgridsthan ongridsofothertwofamilies. Infact,theliftandpitchingmomentcoefficientscomputedontheFamilyIgridsappear to be converging to values different from those computed on the Family II grids. This discrepancy motivated the introductionofanintermediateFamilyIII.FamilyIIIsolutionsoncoarsergridsappeartobeconvergingtoyetanother limit,butonthefinestgridturntowardthevaluescomputedontheFamilyIIgrids. Thisbehaviorisobservedforall codes and indicates that the solution variations due to differences in the trailing-edge resolution are larger than the variationsduetodifferencesindiscretizationschemes. GreendottedlinesinFigs.7–8showthevaluescorresponding to the infinite grid refinement computed by a linear extrapolation fitting the two finest grids. On grids of Family I, the extrapolated lift coefficients vary between values of 1.0885 and 1.0905 and the extrapolated pitching moment coefficients vary between 0.0069 and 0.0074. On grids of Family II, the extrapolated lift coefficient is 1.0910 and theextrapolatedpitchingmomentcoefficientis0.0068. Notethattheliftandpitchingmomentcoefficientscomputed frompresumablythemostaccuratesolutionsonthefinestFamilyIIgridlieoutsideoftherangespannedbythelift andpitchingmomentvaluesextrapolatedfromsolutionsongridsofFamilyIandFamilyIII. (a)FamilyI. (b)FamilyII. (c)FamilyIII. Figure4. NACA0012:Gridconvergenceofthetotaldragcoefficient(CD). (a)FamilyI. (b)FamilyII. (c)FamilyIII. Figure5. NACA0012:Gridconvergenceofthepressuredragcoefficient(CDp). Figures 9–11 show variations of forces and moment computed on grids of Family II with respect to h2 = 1/N. The results are shown for forces and moment computed over the full airfoil and over the areas near the trailing and leadingedges. ThelocalintegrationareasaredefinedinTable2. Therightendoftheleading-edgeintegrationinterval is selected as the x-coordinate of the surface node on the 897×257 grid nearest to x = 0.1. Analogously, the left endofthetrailing-edgeintegrationintervalisselectedasthex-coordinateofthesurfacenodeonthesame897×257 gridnearesttox = 0.9. Theseendnodesarepresentonfourfinergrids. Thecontributionstotheforcesandmoment 8of50 AmericanInstituteofAeronauticsandAstronautics (a)FamilyI. (b)FamilyII. (c)FamilyIII. Figure6. NACA0012:Gridconvergenceoftheviscousdragcoefficient(CDv). (a)FamilyI. (b)FamilyII. (c)FamilyIII. Figure7. NACA0012:Gridconvergenceoftheliftcoefficient(CL). (a)FamilyI. (b)FamilyII. (c)FamilyIII. Figure8. NACA0012:Gridconvergenceofthepitchingmomentcoefficient(CM). aremuchlargerintheleading-edgeregionthaninthetrailing-edgeregion. TheresultsforCFL3Dcomputationson quadrilateralgridsandforFUN3Dcomputationsontriangulargridsareplottedonlyforthefullairfoil. The convergence plots for lift, moment, and pressure drag are almost linear over the three finest quadrilateral grids for all three codes, indicating apparent second-order convergence. Lift and moment computed by FUN3D on 9of50 AmericanInstituteofAeronauticsandAstronautics Table2. Leadingandtrailing-edgeintegrationareas. LeadingEdge TrailingEdge 0.0≤x≤0.100177952877727 0.899166466843597≤x≤1.0 triangulargridsoverthefullairfoilshowapparentsecond-orderconvergenceoverthethreecoarsergridsandasharp turnonthefinestgrid. ThevaluesofC andC computedbyFUN3Donthefinesttriangulargridareclosetothe L M valuescomputedbyFUN3DandTAUonthefinestquadrilateralgrid. TheC coefficientcomputedbyFUN3Don Dp triangular grids converges with an apparent order higher than second. Convergence plots for the viscous drag show less than second-order convergence for the drag computed by FUN3D and TAU over the full airfoil and over the leading-edgearea. Variationsofdraginthetrailing-edgeareaappearverysmall. Theextrapolated,grid-refinedvalues of aerodynamic coefficients computed with different codes are not the same. CFL3D extrapolates lift and pitching momenttovaluessomewhatdifferentfromthevaluesextrapolatedbyFUN3DandTAU.Thesediscrepanciesmaybe a result of differences in implementation of the SA turbulence model. CFL3D employs a thin-layer approximation forthediffusiontermandastandardSAformulationthatdoesnotallownegativevaluesfortheturbulencevariables; FUN3DandTAUuseafull-diffusionapproximationandtheSA-negvariantoftheSAmodel. Theextrapolatedvalues oftheliftandmomentinthetrailing-edgeareashowsomedifferencesbetweenFUN3DandTAUsolutionsaswell. D. SurfacePressureandSkinFriction ThissectioncomparesthesurfacepressureandskinfrictioncoefficientsfromtheFUN3D,CFL3D,andTAUsolutions onthefinest7,169×2,049gridofFamilyII.Inmoderatelyzoomedviewsfocusedontheleadingandtrailingedges (Fig. 12), the solutions are indistinguishable. Only with a super zoom (Fig. 13) do some differences come to light. Figures13(a)and(b)comparesolutionsclosetotheminimumpressureandthemaximumskinfrictionlocationsnear the leading edge. The CFL3D solution shows a smaller pressure and less skin friction than the FUN3D and TAU solutions, whichareindistinguishablyclosetoeachother, evenonthesuper-zoomview. Thelargestdifferencesare observed in the immediate vicinity of the trailing edge (Figs. 13 (c) and (d)); FUN3D and TAU solutions indicate a smallareaofapositiveload,whiletheCFL3Dsolutionindicatesanegativeloadinthisarea. TheTAUsolutionshows amoreoscillatorysurfacepressure,especiallyonthelowersurface,thanothertwosolutions. Thenear-trailing-edge maximumsofthelower-surfacepressureintheCFL3DandTAUsolutionsarecomparableandlargerthanthatinthe FUN3Dsolution. CFL3DandTAUshowasmallareaofnegativeskinfrictionintheimmediatevicinityofthetrailing edgeindicatingflowseparation;FUN3Dshowsnoflowseparation. Althoughnotshown,FUN3Dsolutionsoncoarse gridsalsohavesomeflowseparation,butitgoesawaywithgridrefinement. 10of50 AmericanInstituteofAeronauticsandAstronautics