Towards a parameter-free method for high Reynolds number turbulent flow simulation based on

15 Pages·2015·0.91 MB·English

by  Johan Hoffman

Availableonlineatwww.sciencedirect.com
ScienceDirect
Comput.MethodsAppl.Mech.Engrg.288(2015)60–74
www.elsevier.com/locate/cma
Towards a parameter-free method for high Reynolds number
turbulent flow simulation based on adaptive finite element
approximation
Johan Hoffmana,b,∗, Johan Janssona,b, Niclas Janssonc, Rodrigo Vilela De Abreua
aDepartmentofHighPerformanceComputingandVisualization,KTHRoyalInstituteofTechnology,10044Stockholm,Sweden
bBasqueCenterforAppliedMathematics,Mazarredo,14.48009Bilbao,BasqueCountry,Spain
cRIKENAdvancedInstituteforComputationalScience,Kobe,Japan
Availableonline29December2014
Highlights
• Wepresentaparameter-freemethodforhighReynoldsnumberturbulentflow.
• Themethodisbasedonadaptivemeshoptimizationusingadjointtechniques.
• Weconnectthenumericalapproximationtodissipativeweaksolutions.
• Turbulentboundarylayersareleftunresolved,noboundarylayermeshisneeded.
• Wehighlightbenchmarkproblemsforwhichthemethodisvalidated.
Abstract
Thisarticleisareviewofourworktowardsaparameter-freemethodforsimulationofturbulentflowathighReynoldsnumbers.
InaseriesofpaperswehavedevelopedamodelforturbulentflowintheformofweaksolutionsoftheNavier–Stokesequations,
approximated by an adaptive finite element method, where: (i) viscous dissipation is assumed to be dominated by turbulent
dissipationproportionaltotheresidualoftheequations,and(ii)skinfrictionatsolidwallsisassumedtobenegligiblecomparedto
inertialeffects.Theresultisacomputationalmodelwithoutempiricaldata,wheretheonlymodelparameteristhelocalsizeofthe
finiteelementmesh.Underadaptiverefinementofthemeshbasedonaposteriorierrorestimation,outputquantitiesofinterestin
theformoffunctionalsofthefiniteelementsolutionconvergetobecomeindependentofthemeshresolution,andthustheresulting
methodhasnoadjustableparameters.Noadhocdesignofthemeshisneeded,insteadthemeshisoptimizedbasedonsolution
features,inparticularnoboundarylayermeshisneeded.Weconnectthecomputationalmethodtothemathematicalconceptofa
dissipativeweaksolutionoftheEulerequations,asamodelofhighReynoldsnumberturbulentflow,andwehighlightanumber
ofbenchmarkproblemsforwhichthemethodisvalidated.Thepurposeofthearticleistopresentthecomputationalframeworkin
aconciseform,toreportonrecentprogress,andtodiscussopenproblemsthataresubjecttoongoingresearch.
⃝c 2014ElsevierB.V.Allrightsreserved.
MSC:65M60;65M50
Keywords:Finiteelementmethod;Adaptivemeshrefinement;Turbulencesimulation
∗Corresponding author at: Department of High Performance Computing and Visualization, KTH Royal Institute of Technology, 10044
Stockholm,Sweden.
E-mailaddresses:[email protected](J.Hoffman),[email protected](J.Jansson),[email protected](N.Jansson),[email protected]
(R.V.DeAbreu).
URL:http://www.csc.kth.se/∼hoffman(J.Hoffman).
http://dx.doi.org/10.1016/j.cma.2014.12.004
0045-7825/⃝c 2014ElsevierB.V.Allrightsreserved.
J.Hoffmanetal./Comput.MethodsAppl.Mech.Engrg.288(2015)60–74 61
1. Introduction
InthispaperwereviewourworktowardsageneralframeworkforsimulationofunsteadyhighReynoldsnumber
turbulent flow using adaptive finite element methods, with the purpose to give a coherent presentation of the com-
putational framework, to report on recent progress, and to highlight challenges and open problems. We can refer to
themethodologyasDirectFiniteElementSimulation(DFS),andwerecallthetheoreticalfoundationandshowhow
DFSofferssolutionstothemainchallengesofturbulencesimulation:howtoassesstheaccuracyinasimulation,how
tooptimizetheuseofcomputationalresources,andhowtominimizethedependencyonapplicationspecificmodel
parameters.
Overtheperiod2010–2013,wehaveparticipatedinanumberofComputationalFluidDynamics(CFD)benchmark
workshops organized by the aerospace community (AIAA, NASA, DLR, Boeing, etc.), namely the BANC series
[1–4](WorkshoponBenchmarkProblemsAirframeNoiseComputation)andHiLiftPW-2[5](2ndAIAACFDHigh
Lift Prediction Workshop). These workshops have gathered leading research groups in the area, as well as major
aerospace companies (e.g. Boeing, Cessna, Aerospace Kawasaki) and the main commercial software developers
(e.g. ANSYS and Exa), with the objective to assess the state of the art computational methods for high Reynolds
number aerodynamics and aeroacoustics in complex geometry. From the results presented in these workshop it is
clear thatour CFDmethodology is uniquein several aspects,in particular: (i)it was theonly method thatprovided
estimatesoftheerrorinthecomputation,(ii)itwastheonlymethodwithautomaticmeshadaptation,and(iii)itwas
theonlymethodwhichwasfreefrommodelparameters.
IntheDFSframework,(i)and(ii)areachievedbygoal-orientedaposteriorierrorestimation,wheretheerrorina
functionalofinterestcanbeexpressedintermsoftheresidualoftheapproximatesolutionweightedbythesolution
of an associated adjoint (dual) problem. A general framework for a posteriori error estimation based on duality for
differentialequationswasdevelopedinthe1990sdescribed,e.g.,in[6–8]andthereferencestherein.Ourownwork
hasbeenfocusedontheextensionofthisframeworktoCFDandturbulentflow[9–11].Themethodisvalidatedfor
arangeofbasicbenchmarkproblemsincludingflowpastcylinders,spheresandsquares[12–15],wheretheadaptive
algorithmisshowntoconvergetotheexperimentalreferencedata.
The parameter-free aspect (iii) of DFS is related to the fact that no subgrid scale stresses or Reynolds stresses
areparametrizedasinLES(LargeEddySimulation)orRANS(ReynoldsAveragedNavier–Stokesequations)[16],
andthatforhighReynoldsnumbersweassumethattheskinfrictionoftheboundarylayersisdominatedbyinertial
effects,sothatafreeslipboundaryconditionisusedatsolidboundaries[17–19].IncontrasttoRANSandLES,no
averagingorfilteringoperatorsareintroducedinDFS,theresolutionoftheturbulentscalesisdeterminedaspartof
thecomputation,withthelocalmeshsizeadaptedbasedontherequiredprecisioninthechosenquantityofinterest.
Themotivationforthefreeslipboundaryconditionisbasedontheuseofamixedboundarycondition,forthenormal
velocityandthetangentialstresses,similartostandardwalllayermodelsofLES[20].
We now proceed to put the DFS method into the context of research on computational methods for turbulence
simulation,andmathematicalresearchontheNavier–Stokesequations.
1.1. Turbulencesimulation
Turbulenceiscentralforourunderstandingoftheworld,includingthegeophysicalflowoftheocean-atmosphere
system, and the aerodynamics of vehicles, aircraft and wind turbines, for example. Simulation of turbulent flow is
anoutstandingproblemthatcutsthroughmathematics,numericalanalysis,scientificcomputingandfluidmechanics.
Thisinterdisciplinarynatureoftheproblemmakesitparticularlychallenging.
ThefundamentalmathematicalmodeloffluidflowisthenonlinearNavier–Stokesequationsin3spacedimensions,
forwhichaproofofexistenceofclassical(smooth)solutionsinthegeneralcaseismissing.Tocomputeapproximate
solutionsoftheNavier–Stokesequationsinthecaseofturbulentflowisextremelyexpensive;anheuristicargument
estimatesthenumberofdegreesoffreedomneededtoresolveallturbulentscalesinaDirectNumericalSimulation
(DNS) to be of the order Re9/4, with Re the Reynolds number, which for many decades to come will make DNS
impossiblee.g.forgeophysicalfloworflightaerodynamicswhere Re≫106.
Areductioninthenumberofdegreesoffreedomispossiblebyinsteadconsideringmodelsforaveragedquantities,
suchasstatisticalmeanfieldsinRANSorspatiallyfilteredfieldsinLES,atthecostofintroducingturbulence/subgrid
62 J.Hoffmanetal./Comput.MethodsAppl.Mech.Engrg.288(2015)60–74
modelstomodelunresolvedstatisticalfluctuations/subgridscales[16].Thisclosureproblemofturbulencemodeling
hasproventobeanoutstandingproblem,wheremodelparametersoftenhavetobetunedtoaparticularapplication.
Similarly,subgridmodelsintroduceparametersthatneedtobechosenbytheuser,inparticularnearsolidwallsitisa
challengetoselectproperparameterstoaccuratelycapturephenomenasuchasflowseparation[20].
In RANS/LES simulations, the relation between model errors and numerical errors is delicate, since the effect
ofthenumericalmethodoftenissimilartotheeffectofturbulence/subgridmodels,e.g.intheformofdissipationof
kineticenergy.Forexample,thewidelyusedSmagorinskysubgridmodel[21]iscloselyrelatedtomethodsofartificial
viscosityusedtostabilizenumericalmethods[22].Totacklethisissue,differentapproacheshavebeenattempted,from
usinghighordernumericalmethodswithminimaldissipation,tointerpretingthenumericalstabilizationasthesubgrid
model[23,11,24,25],oralternativelytryingtobalancethetwosourcesoferrors[26,10].Therelationshipbetweenthe
localresolutionofthecomputationalmeshandtheturbulentlengthscalesofRANSandLESalsoposeschallenges,in
particularnearsolidwallswheresmallturbulentscalesdictateaveryfinemeshresolutionwhichdominatesthetotal
computationalcost[20].
Inthecontextoffiniteelementmethods,VariationalMultiscaleStabilizedmethods(VMS)[24]havebeendevel-
oped where the variational form of the Navier–Stokes equations is split into resolved scales and subgrid scales, by
decomposingtheunderlyingHilbertspaceintoadirectsumofacoarsespaceandafinespace,whereasubgridmodel
isusedtomodelthefinescalecontribution.
1.2. Dissipativeweaksolutions
Intheframeworkpresentedinthispaperwefollowadifferentpath,whereweavoidtheprocessofaveragingand
scaleseparationaltogether,andthusalsotheclosureproblem.ForincompressibleflowthemathematicianJeanLeray
in1934[27]provedtheexistenceofaweaksolution,orturbulentsolution(solutionturbulente)intheterminologyof
Leray, which is a function that satisfies a weak form of the Navier–Stokes equations (in the sense of distributions).
TheworkofLeraylaidthefoundationforthemodernmathematicaltheoryoftheNavier–Stokesequations,andthe
regularized form of the equations used by Leray later inspired a specific class of LES methods [28]. The question
of uniqueness of a weak solution is still an open problem, and whether a weak solution is also a smooth (strong)
solutionisformulatedasa$1millionClayPrizeproblem[29].Aweaksolutiondoesnotsatisfythestandardenergy
balancevalidforasmoothsolutionoftheNavier–Stokesequations,insteadonlyanenergyinequalitycanbeshownto
hold.Partialregularityofaweaksolutionwasprovenundertheassumptionthatageneralizedenergyinequalityholds
[30,31],inwhichcasetheweaksolutionisreferredtoasasuitableweaksolution.Itcanbeshownthatcertaindiscrete
approximations,includingfiniteelementapproximations,convergetosuitablesolutions[32].
Foraweaksolution,dissipationwithoutviscosityispossible.Inparticular,DuchonandRobertin2000introduced
the notion of a dissipative weak solution as a model for high Reynolds number turbulent flow, where an explicit
expressionforthesourceofinvisciddissipationwasderived,whichtooktheformofalocallackofregularityinthe
velocityfield[33].Alreadyin1949,OnsagersuggestedthatweaksolutionsoftheinviscidEulerequations(assuming
suchexist)couldbeusedtomodelhighReynoldsnumberturbulentflow[34]:
“Itisofsomeinteresttonotethatinprinciple,turbulentdissipationasdescribedcouldtakeplacejustasreadily
withoutthefinalassistancebyviscosity.Intheabsenceofviscosity,thestandardproofoftheconservationofenergy
doesnotapply,becausethevelocityfielddoesnotremaindifferentiable!”
A consequence of the idea that turbulent dissipation is independent of viscosity, is that it should be possible to
model turbulence without any empirical parameters. Onsager also conjectured that in such “ideal” turbulence, the
regularity of the velocity field could at most satisfy a Ho¨lder condition with exponent 1/3, otherwise the energy
wouldbeconserved.SeveralproofsofOnsager’sconjecturehavesincethenbeenpresented,andintheir2006review
oftheworkofOnsager,EyinkandSreenivasanendwithaconclusiononthefutureofOnsager’sideas[35]:
“We believe that Onsager’s theoretical vision of an ideal turbulence described by inviscid fluid equations is a
properidealizationforunderstandinghighReynoldsnumberflows.Needlesstosay,inrealphysicalturbulencethere
isviscosity,whichisalwayspositive.[...]Inthesameway,thezero-viscositylimit,whichsupposesaninfinitenumber
of cascade steps, should be a good idealization for turbulence with a large but finite number of cascade steps, that
is, a Reynolds number which is large but finite. The vindication of this belief, if it is true, must come from a set of
calculationaltoolsforthezero-viscositylimit,whichwillmakeit,intheend,atrulypredictivedevice”.
J.Hoffmanetal./Comput.MethodsAppl.Mech.Engrg.288(2015)60–74 63
1.3. Directfiniteelementsimulation
In [11,36,15] we show that approximate dissipative weak solutions can be computed to simulate turbulent flow
usingastabilizedfiniteelementmethodthatsatisfiescertainconditionsonstabilityandconsistency,whichwerefer
to as a General Galerkin (G2) method. We also refer to the computational methodology as DFS, as an extension
to turbulent flow of the general framework for a posteriori error estimation developed in the early 1990s [6–8]. In
particular, for high Reynolds numbers we can set the viscosity to zero, thus realizing a computational model for
turbulentflowfreefromempiricalparameters.
Wearguethatthepropermathematicalobjectsforassessingwell-posednessinDFSarefunctionalsofthecomputed
weaksolutions[37].Torepresenttheerrorinafunctionalofadissipativeweaksolutionweintroduceadual(adjoint)
linearized problem, which opens for a posteriori estimation of the error in outputs of interest, such as the drag and
liftofanairplane.Theaposteriorierrorrepresentationcanbeusedtoconstructefficientadaptivealgorithms,where
the computational mesh and the finite element space are optimized for the computation of a particular output of
interest(goalfunctional)toacertainaccuracyusingminimalcomputationalresources.Thusalsothelocalresolution
of turbulent scales is determined as part of the computation (a posteriori) based on what functional is the goal of
the simulation. For example, approximation of the aeroacoustic sources in a turbulent flow past a landing gear may
demandadifferentmeshresolutionthanapproximationoftheaerodynamicforcesonthesamelandinggear[1,2].
A classical problem of great significance is the prediction of aerodynamic forces acting on a body in a high
Reynolds number flow. The choice of boundary conditions in the simulation is here central. As already noted, for
realistic applications e.g. in flight aerodynamics, full resolution of a turbulent boundary layer in DNS is impossible
duetothecomputationalcost.InLES,typicalmodelsofthenearwallturbulentflowarewallshearstressmodels[20]
orlocalRANSmodelsnearthewall[38],wherebothapproachesdemandhighresolutionboundarylayermeshesand
calibrationofmodelparameters.
In [18] we investigate a simple wall shear stress model parametrized by the skin friction stress generated by the
turbulent boundary layer. We find that for very high Reynolds numbers, corresponding to small skin friction, the
macroscopicflowislargelyindependentofthewallshearstressmodel.Inparticular,theskinfrictionparametercan
be set to zero, corresponding to a free slip boundary condition, and no boundary layer mesh is needed. For high
Reynoldsnumberflow,DFSwithafreeslipboundaryconditionhasbeenshowntobeanefficientandaccuratemodel
forapproximationofaerodynamicforcesandaeroacousticsources[1–3,5].
In this paper the G2/DFS framework is presented in Section 2, and we give a brief description of a detailed
implementationofthemethodinSection3.InSection4wesummarizeasetofvalidationstudiesforthecomputational
model,inSection5wediscusschallengesandopenproblemsrelatedtothedualproblem,andinSection6wegivea
shortsummary.
2. TheG2/DFScomputationalframework
LowMachnumberflowcanbeapproximatedbyincompressibleflow,andthefocusofthispaperisthemodelof
incompressiblehighReynoldsnumberturbulentflowaroundblufforstreamlinedbodies,wheretheReynoldsnumber
is defined as Re = UL, withU a characteristic velocity, L a characteristic length and ν is the kinematic viscosity.
ν
Further,weassumethatturbulentboundarylayersarefullydeveloped,sothattransitiontoturbulenceintheboundary
layersisnotpartofthemodel.
We choose mixed boundary conditions in the form of prescribing the normal component of the velocity and the
tangential components of the viscous shear stresses, with the purpose to avoid resolution of thin boundary layers
at high Reynolds numbers. The normal velocity is set to zero relative the body, corresponding to a no penetration
boundarycondition,whereasthetangentialshearstressesmustbeprovidedasdatatotheproblem.Boundarystress
data is typically not available, but viscous shear stresses, measured in the form of a non-dimensional skin friction
coefficient,isknowntodecreasewithincreasingReynoldsnumber[39].InLES,wall-layermodelstakeasimilarform,
withwallshearstressproportionaltoskinfriction[20].Thus,thisformofboundaryconditionscanbemotivatedboth
bystandardargumentsfromtheLESliterature[16],notingthatbothLESandDFSconcernmodelsoftheunresolved
turbulent flow near the boundary, and by strict mathematical arguments choosing mixed Dirichlet–Neumann type
boundaryconditionsfortheNavier–Stokesequations.
64 J.Hoffmanetal./Comput.MethodsAppl.Mech.Engrg.288(2015)60–74
2.1. Mathematicalmodel
The Navier–Stokes equations for incompressible flow of (small positive) constant kinematic viscosity ν in a
(bounded) volume Ω in R3 around a fixed rigid body with smooth boundary Γ over a time interval I = [0,T],
taketheform:Findthevelocityu =(u ,u ,u )andpressure pdependingon(x,t)∈Ω ∪Γ ×I,suchthat
1 2 3
u˙ +(u·∇)u+∇p−∇·τ =0 inΩ ×I,
∇·u =0 inΩ ×I,
un =0 onΓ ×I, (1)
τ =βu onΓ ×I,
sk tk
u(·,0)=u0 inΩ,
with k = 1,2,u˙ = ∂u,u = u · n the fluid velocity normal to Γ with n a unit outward normal vector, τ =
∂t n
τ(u)=2νϵ(u)theviscousstresswithϵ(u)thestandardvelocitystrain,τ thetangentialstressesandu =u·t the
sk tk k
tangentialvelocitieswitht orthogonaltangentvectors,andu0agiveninitialcondition.
k
We assume suitable far-field inflow/outflow boundary conditions, and we note the Dirichlet–Neumann boundary
conditioninthenormalandtangentialdirections,respectively,withtheeffectofturbulentboundarylayersmodeled
bya(skin)frictionparameterβ.Withzeroviscosityν andboundaryfrictionβ,themodel(1)reducestotheinviscid
Eulerequations.
In [33] it is shown that for any weak solution (u,p) of the Euler equations (assumed to exist), a local energy
equationissatisfied(inthesenseofdistributions)
d 1   1 
|u|2 +∇· u |u|2+ p +D(u)=0, (2)
dt 2 2
with D(u)definedintermsofthelocalsmoothnessofu.If D(u)≥0,then(u,p)isreferredtoasadissipativeweak
solution.
2.2. Finiteelementmethod
InG2/DFSwecomputeapproximatesolutionsoftheNavier–Stokesequations(1)byaweightedresidualstabilized
finiteelementmethod,oftheform:FindUˆ =(U,P)∈ V suchthatforallvˆ =(v,q)∈ V
h h
[R(U;Uˆ),vˆ]+[hR(U;Uˆ),R(U;vˆ)]=0, (3)
whereV isaspace–timefiniteelementspacewithvelocitiesvsatisfyingv·n =0onΓ,[·,·]isan L (Ω ×I)inner
h 2
product, R(U;vˆ) = (v˙ +U ·∇v +∇q,∇ ·v) is the residual, and h is the local mesh size. The first term in (3)
ˆ
establishesU asaweaksolutionof(1)andthesecondtermintroduceskineticenergydissipation
D (Uˆ)=[hR(U;Uˆ),R(U;Uˆ)]=∥h0.5R(U;Uˆ)∥2 (4)
h L2
boundedbydatawith∥·∥ an L (Ω ×I)norm[36].
L2 2
Noticethat here ν and β areset tozero withinsteadthe weightedresidual stabilizationintroducinga dissipative
effectasanautomaticturbulencemodel.ThisisanalogoustothedissipativeweaksolutionsintroducedbyDuchonand
Robert[33],withdissipationcausedbyalackofsmoothnessinthevelocityfield,unrelatedtoviscousdissipation.For
asmoothsolutionUˆ,correspondingtoasmallresidualin L (Ω × I),thedissipativeeffectoftheweightedresidual
2
stabilizationvanishes.InSection4werecallvalidationresultsthatshowthatDFSbasedontheEulerequationsisa
goodapproximationofhighReynoldsnumberturbulentflow.Wealsofindthatwithsufficientmeshresolution,D (Uˆ)
h
in(4)isindependentofthelocalmeshsizeh [15].
InSection3wegiveaprecisedefinitionofaparticularimplementationofG2/DFS,butmanychoicesarepossible
aslongascertainconditionsforstabilityandconsistencyaresatisfied[36,15].
J.Hoffmanetal./Comput.MethodsAppl.Mech.Engrg.288(2015)60–74 65
2.3. Boundaryconditions
Itiswellknownthatforviscousflow,thinboundarylayersdevelopnearsolidwallsoverwhichthevelocitychanges
fromthevelocityofthewalltothefreestreamvelocity.ThethicknessoftheseboundarylayersscalewiththeReynolds
number, and at high Reynolds numbers the boundary layers are so thin that computational resolution is impractical
or impossible. At sufficiently high Reynolds numbers the boundary layers undergo transition to turbulence, where
computationalresolutionofturbulentboundarylayersisevenlessrealisticinanyapplicationofsignificance[20].For
example,fullcomputationalresolutioninDNSoftheturbulentflowaroundacaroranairplane,includingboundary
layers,isimpossibleformanydecadestocome[40].
Ever since the introduction of the boundary layer concept by Prandtl [41], models of the flow decoupling the
boundary layer from the interior flow have been pursued, in which the effect of the boundary layer is modeled as a
boundaryconditionfortheinteriorflow[20].In(1)wemodeltheeffectofturbulentboundarylayersasasmallskin
frictionstressatthesolidboundary,whichdissipateskineticenergyproportionaltoβ.Thesensitivityofthesolution
with respect to perturbations in the skin friction parameter β is investigated in [17,18], where it is found that for
smallβ (correspondingtohighReynoldsnumbers)thesensitivityislow.Inparticular,forexternalflowatReynolds
numbers higher than Re > 106 we typically approximate the small skin friction by zero, which makes DFS into a
parameter-freemodelofhighReynoldsnumberexternalflow.
We have presented evidence in the form of analysis and computational studies, that the macroscopic features of
highReynoldsnumberflow,including3Dflowseparation,ismainlydeterminedbylargescalestabilityaspects,rather
thanboundarylayereffects[19,1–3,5].
2.4. Aposteriorierrorestimation
Thegoalofasimulationcanoftenbeformulatedascomputinganumberofoutputquantities,forexamplethelift
and drag of an airplane. To estimate the error in an output functional of interest we use a posteriori error analysis
basedonsensitivityinformationfromthesolutionofaadjoint(dual)problem.
WedefineatargetfunctionalM(uˆ)=[uˆ,ψˆ]+[pn,ψ ] ,whereψˆ =(ψ,χ)isagivenweightfunctionactingas
Γ Γ
dataforthedualproblem,andψ isboundarydataforthedualvelocity,with[·,·] a L (Γ × I)innerproduct.For
Γ Γ 2
example,tochoosethesumoftheliftanddragofanairplaneasthetargetfunctionalwesetψˆ =0andψ =v +v ,
Γ D L
withv andv unitvectorsindirectionsoppositeandnormaltotheflightdirection,respectively.
D L
Thedifferenceinoutput M(uˆ)− M(Uˆ) = [uˆ,ψˆ]+[pn,ψ ] −[Uˆ,ψˆ]−[Pn,ψ ] oftwoDFSsolutionsuˆ
Γ Γ Γ Γ
ˆ
andU ondifferentmesheswithmaximalmeshsizeh,canberepresentedas
M(uˆ)−M(Uˆ)=[R(u;uˆ)− R(U;Uˆ),ϕˆ] (5)
whereϕˆ =(ϕ,θ)isasolutionoftheadjointlinearizedproblem
−ϕ˙ −(u·∇)ϕ+∇UTϕ+∇θ =ψ inΩ ×I,
∇·ϕ =χ inΩ ×I,
(6)
ϕ =ψ onΓ ×I,
Γ
ϕ(·,T)=0 inΩ,
where(∇UTϕ) =3 ∂U /∂x ϕ .
j i=1 i j i
Ifweassumethatthetimestepisboundedbysomeconstanttimesh,wecanboundtheoutputerrorusingbasic
analysis[36]:
|M(uˆ)−M(Uˆ)|≤C∥h0.5ϕˆ∥ , (7)
H1
with a constant C depending on data and bounds on the magnitude of the discrete velocities u andU. We can thus
guarantee stability in output if ∥h0.5ϕˆ∥ is small, where H1 = H1(Ω × I) is the standard Hilbert space with
H1
associatednorm.Wenote,however,thatforturbulentflowthelossofsharpnessintheboundin(7)maybesignificant,
duetotheuseofaglobalCauchy–Schwarzinequality.Therefore,inadaptivealgorithms,theerrorrepresentationin
(5)isuseddirectly,orwithonlylocalboundsoneachcellinthefiniteelementmesh.
66 J.Hoffmanetal./Comput.MethodsAppl.Mech.Engrg.288(2015)60–74
2.5. Adaptivealgorithms
A posteriori error estimation provides a powerful tool for the development of robust and efficient adaptive algo-
rithms.Thebasicideaofadaptivealgorithmsistooptimizethecomputationalmethodwithrespecttothegoal(output
ofinterest)ofthecomputation.Typicalparametersofanadaptivefiniteelementmethodincludethelocalmeshsize
(h-adaptivity),localdegreeofthefiniteelementapproximation(p-adaptivity),localshapeofthecells(r-adaptivity),
orcombinationsthereof.
ˆ
To construct an adaptive DFS algorithm we compute a solutionU for which we use the error representation (5)
to estimate the output error |M(uˆ)− M(Uˆ)|, with uˆ any DFS solution computed on a finer mesh than Uˆ. We can
eitherchooseuˆ andUˆ tobetwoDFSsolutionscomputedondifferentmeshesinthemeshhierarchygeneratedbythe
adaptivealgorithm,oralternativelytoconsideruˆ tobeanidealsolutioncomputedonafinestmeshpossibleinwhich
caseweusetheapproximationu ≈U inthedualproblem(6).
An adaptive algorithm is based on an estimate of the local contribution to the global error, which is obtained by
splittingtheerrorrepresentation(5)intoasumoferrorindicatorsEn,K overthecells K inthediscretespacemesh
Th,andthetimeintervals In =(tn−1,tn),n =1,...,N,thatis
N N
M(uˆ)−M(Uˆ)=  En,K =  [R(u;uˆ)− R(U;Uˆ),ϕˆ]n,K (8)
n=1K∈Th n=1K∈Th
with[·,·]n,K an L2(K ×In)innerproduct.
2.6. Approximationofthedualsolution
Thereareseveraltechnicalaspectswithrespecttohowtoapplytheerrorrepresentationformula(8)inanadaptive
algorithm,inparticularwithregardstotheapproximationofϕˆ,thesolutiontothedualproblem.
Thedualproblemislinear,butrunsbackwardintimewhichisachallengesincetheDFSvelocitiesuandU actas
coefficientsinthedualproblemandthushavetobestoredasdata.Inpractice,thedualproblemisoftensolvedbyusing
asimilarfiniteelementmethodasusedfortheprimalDFSproblem,andthecoefficientsuandU canbeinterpolatedin
timetominimizedatastorage,withpossiblyarestartschemetoincreaseaccuracy.Also,severaldifferentapproaches
are possible for how to apply the error representation formula (8) in practice [42], to retain sharpness and achieve
maximalrobustness.
When we use the approximation u ≈ U in the dual problem we introduce a linearization error which is hard to
estimate a priori, since worst case estimates based on Gro¨wall’s inequality grossly overestimate the effect [36]. In
practice, computations on a sequence of meshes with different associated approximationsU provide an estimate of
theasymptoticeffectofthislinearizationerrorontheresultingaposteriorierrorestimates.
More specifically, the accuracy in the a posteriori error estimates depends on the accuracy in the approximation
ofthedualproblem,whichisaffectedby(i)perturbationsindataand(ii)errorsfromcomputationalapproximation.
Aspartoftheadaptivealgorithm,approximationsofthedualproblemarecomputedonahierarchyofcomputational
meshes, which can be used as a robustness test with respect to (i)–(ii), improving the reliability of the a posteriori
errorestimates.
3. DFSimplementation
We now give a detailed presentation of one implementation of DFS, in the form of the cG(1)cG(1) method with
piecewiselinearapproximationinspaceandtimeovertetrahedralmeshes.InSection4werecallvalidationstudiesof
thismethod.
3.1. ThecG(1)cG(1)method
The cG(1)cG(1) method is based on the continuous Galerkin method cG(1) in space and time. With cG(1) in
time the trial functions are continuous piecewise linear and the test functions piecewise constant. cG(1) in space
J.Hoffmanetal./Comput.MethodsAppl.Mech.Engrg.288(2015)60–74 67
corresponds to both test functions and trial functions being continuous piecewise linear. Let 0 = t < t < ··· <
0 1
tN = T beasequenceofdiscretetimestepswithassociatedtimeintervals In = (tn−1,tn)oflengthkn = tn −tn−1
andspace–timeslabsS =Ω ×I ,andletW ⊂ H1(Ω)beafiniteelementspaceconsistingofcontinuouspiecewise
n n
linear functions on a tetrahedral mesh Th of mesh size h(x) with W the functions v ∈ W satisfying the Dirichlet
0
boundaryconditionv·n| =0.
Γ
WeseekUˆ =(U,P),continuouspiecewiselinearinspaceandtime:Forn =1,...,N,find(Un,Pn)≡(U(t ),
n
P(t ))withUn ∈ Vn ≡[W ]3and Pn ∈ W,suchthat
n 0 0
((Un −Un−1)k−1+U¯n ·∇U¯n,v)−(Pn,∇·v)
n
+(∇·U¯n,q)+SDn(U¯n,Pn;v,q)=0 ∀vˆ =(v,q)∈ Vn ×Wn, (9)
δ 0
whereU¯n = 1(Un +Un−1)and Pn arepiecewiseconstantintimeover I ,withthestabilizingterm
2 n
SDn(U¯n,Pn;v,q)≡(δ (U¯n ·∇U¯n +∇Pn),U¯n ·∇v+∇q)+(δ ∇·U¯n,∇·v), (10)
δ 1 2
and


(v,w)= v·wdx,
K∈Th K
withthestabilizationparametersδ =κ (k−2+|Un−1|2h−2)−1/2 andδ =κ h|Un−1|,whereκ andκ arepositive
1 1 n 2 2 1 2
constantsofunitsize.Wechooseatimestepsizekn =CCFLminx∈Ω h/|Un−1|,withCCFL aconstantofunitsize.
Remark1. The cG(1) discretization in time does not allow for a fully consistent stabilization of the method, since
the velocity test functions are constant in time. To make the method fully consistent while keeping a least squares
stabilizationoftheresidual,wewouldneedtousetestfunctionsthatarelineardiscontinuousintime,corresponding
to a cG(1)dG(1) method. The main motivation for using cG(1)cG(1) is the reduction of the number of degrees of
freedombyafactor2.
3.2. cG(1)cG(1)energyequation
We recall that a cG(1)cG(1) solution satisfies a local energy equation, with dissipation directly proportional to
the local residual, with the proof given in [15]. We define a family of smooth positive test functions φ (x) for
n
n = 1,...,N, with local compact support supp φ ⊂ Ω, and we define a piecewise constant function in time
n
definedbyφ(x,t)=φ (x)fort ∈ I ,normalizedsuchthat
n n
N 

φ (x)dx k =1.
n n
n=1 Ω
Theorem2. Notingthatδ ≤Ch,fori =1,2,wehavethefollowinglocalenergyestimateforcG(1)cG(1):
i

 N  1(|Un|2−|Un−1|2)k−1+∇·U¯n1|U¯n|2+ Pnφ dx k
 2 n 2 n n
 n=1 Ω

+ N  (δ |R¯ (U¯n,Pn)|2+δ |R¯ (U¯n)|2)φ dx k 
1 1 2 2 n n 
n=1 Ω 
≤Ch1/2
max,φ,n
withhmax,φ,n ≡maxn:suppφn̸=∅(maxx∈suppφnh(x)),andwiththeresiduals R¯1(U¯n,Pn)and R¯2(U¯n)definedby
R¯ (v,q)=v˙+U¯n ·∇v+∇q, (11)
1
R¯ (v,q)=∇·v, (12)
2
fort ∈ I .
n
68 J.Hoffmanetal./Comput.MethodsAppl.Mech.Engrg.288(2015)60–74
Ifwechooseatestfunctionφ withasuitablelocalsupport,andwithφ = 0forall j ̸= n,Theorem2statesthat
j
fort ∈ I ,acG(1)cG(1)solutionuptothefactorCh0.5 weaklysatisfies(integratedagainstasmoothpositivetest
n max,φ,n
functionwithlocalsupport)alocalenergyequationoftheform:
d 1   1 
|Un|2 +∇· U¯n |U¯n|2+ Pn =−Dn(Uˆ), (13)
dt 2 2 h
with Dn(Uˆ)≡δ |R¯ (U¯n,Pn)|2+δ |R¯ (U¯n)|2 ≥0,aresidualbasednumericaldissipation,andwith d (1|Un|2)≡
h 1 1 2 2 dt 2
1(|Un|2−|Un−1|2)k−1.
2 n
This local energy equation connects to a dissipative weak Euler solution [33], with inertial energy dissipation
resultingfromlocalnon-smoothnessofthesolution,andnotfromanyviscosity.Thedissipativeterm Dn(Uˆ)in(13)
h
reflectslocalnon-smoothnessinthesolutionidentifiedbytheresiduals,andisobservedtobeindependentofhunder
sufficientmeshresolution[12,15]
3.3. Aposteriorierroranalysisandanadaptivealgorithm
Usingstandardtechniquesofaposteriorierroranalysis,seee.g.[36],wecanderiveanaposteriorierrorestimate
forcG(1)cG(1).
Theorem3. If Uˆ = (U,P) solves (9), uˆ = (u,p) is a weak solution to (1), and ϕˆ = (ϕ,θ) solves (6), then we
have the following a posteriori error estimate for the output M(Uˆ) = [Uˆ,ψˆ] with respect to the reference output
M(uˆ)=[uˆ,ψˆ]:
N
|M(uˆ)−M(Uˆ)|≤  En,K (14)
n=1K∈Th
where
  
En,K ≡ |R1(Uˆ)|K ·ω1dt + |R2(U)|K ω2dt + |SDδn(Uˆ;ϕˆ)K|dt
In In In
with
R (Uˆ)=U˙ +(U ·∇)U +∇P,
1
R (U)=∇·U,
2
whereSDn(·;·) isalocalversionofthestabilizingform(10),andthestabilityweightsaregivenby
δ K
ω1 =Cnk,Kkn|ϕ˙|K,∞+Cnh,KhK|∇ϕ|K,
ω2 =Cnk,Kkn|θ˙|K,∞+Cnh,KhK|∇θ|K,
where h is the diameter of element K in the mesh Th,Ch ,Ck represent interpolation constants, |w| ≡
K n,K n,K K
(∥w ∥ ,∥w ∥ ,∥w ∥ )with∥w∥ =(w,w)1/2,andthedotdenotesthescalarproductinR3.
1 K 2 K 3 K K K
Basedontheaposteriorierrorestimate(14),wedesignabasicadaptivealgorithmwherethemeshiskeptconstant
intime.Startingfromaninitialmeshtheprimalproblemiscomputedforthefulltimeintervalafterwhichthedual
problemissolvedbackwardsoverthesametimeintervalusingasimilarcG(1)cG(1)discretizationasfortheprimal
problem(sametrialandtestspace,leastsquaresstabilizationoftheresidualexceptthetimederivative,andthesame
stabilizationparametersandtimediscretization).TheerrorindicatorsEK = nN=1En,K arethencomputedforeach
cell K in the mesh. The cells with the highest error indicators E are bisected to generate a new refined mesh, for
K
whichthealgorithmstartsoveruntilconvergencein M(Uˆ)isdetected.
TheprimalsolutionU foreachmeshisstoredtobeusedasdataforthedualproblem,wheretheapproximation
u ≈ U is used for the coefficient u in (6). For demanding large scale problems where it is too expensive to store
ˆ
thecompletetimeseriesoftheprimalsolution,snapshotsofU arestoredwithacertainfrequency,andintermediate
solutionfieldsarereconstructedbyinterpolationintime.
J.Hoffmanetal./Comput.MethodsAppl.Mech.Engrg.288(2015)60–74 69
Adaptive algorithm: given initial coarse meshTh,k =0,
0
(1) compute the solutionUˆ to the primal problem usingTh
k
(2) compute the solutionϕˆ to the dual problem usingTh
k
(3) if convergence in M(Uˆ) is observed, or if the estimated error  E is less than a
K∈Th K
given tolerance, then STOP, else
(4) bisect 10% of cells K ∈Th with highest error indicatorsE
k K
(5) setk =k+1, then goto (1)
3.4. Softwareimplementation
The adaptive algorithm and the cG(1)cG(1) solvers for the primal and dual problems are implemented in the
continuum mechanics solver Unicorn [43] that builds on the high performance (HPC) branch of the finite element
problemsolvingenvironmentDOLFIN,partoftheFEniCSproject[44].DOLFIN-HPCisoptimizedfordistributed
memory architectures using a hybrid MPI+OpenMP approach with efficient parallel I/O (MPI I/O) [45,46], where
ParMETIS [47] is used for mesh partitioning. DOLFIN-HPC supports several parallel linear algebra packages, but
mostlyrelyonPETSc[48].
4. ValidationandapplicationofDFS
DFSisbasedonthecomputationoffiniteelementapproximationsonadaptivelyrefinedcomputationalmeshes.In
particular,noadhocmeshdesignisneeded,sincesolutionsfeaturessuchasflowseparationandturbulentwakesare
automaticallydetectedbytheadaptivemethodthroughaposteriorierrorestimationwithsensitivityinformationfrom
thesolutionofadualproblem.Ahierarchyofadaptivelyrefinedmeshesisgeneratedforwhichconvergencecanbe
observedwithrespecttoanoutputofinterest.Themethodhasbeenappliedforanumberofbluffbodyflowproblems
overarangeofReynoldsnumbers,andwehererecallsomevalidationresultsforDFSintheformofthecG(1)cG(1)
method.
4.1. MediumReynoldsnumberflow
AtmediumReynoldsnumbers,Re < 105,boundarylayersarelaminarandcanberesolvedbythecomputational
mesh.Forthiscase,viscouseffectsarenotnegligiblesothattheviscosityiskeptinthemodel(1),andnoslipboundary
conditionsarechosenwherethevelocityissettozeroonthesolidboundaryΓ.DFSintheformofthecG(1)cG(1)
method has been validated for a number of model problems of simple geometry bluff bodies, including a surface
mountedcubeandarectangularcylinder[12,11],asphere[13]andacircularcylinder[14].Ineachcase,convergence
isobservedforoutputquantitiessuchasdrag,liftandpressurecoefficients,andStrouhalnumbers,andtheadaptive
algorithmleadstoanefficientmethodoftenusingordersofmagnitudefewernumberofdegreesoffreedomcompared
toLESmethodsbasedonadhocdesignofthemesh.
4.2. HighReynoldsnumberflow
InhighReynoldsnumberflow,with Re > 106,boundarylayersareturbulentandareinmostcasestooexpensive
to resolve, and must instead be modeled. To accurately predict e.g. aerodynamic forces it is critical to capture the
correct flow separation, which can be connected to the boundary layer, in case flow separation is not triggered by
sharpfeaturesinthegeometry[15].Dragcrisisforacircularcylinderisanillustrativeexample,whereflowseparation
movesdownstreamthecylinderwallasaconsequenceoftransitiontoturbulenceintheboundarylayer[49].In[18]we
modeldragcrisisbydecreasingtheskinfrictionparameterβinthemodel(1),wherewefindthattheDFSsimulations
reproducetheobserveddragcrisisscenariofromexperiments[50–53].
In particular, we find that for small β the solution is insensitive to the particular value of β, so that we reach
an ultimate regime where the turbulent boundary layer is modeled by a zero friction free slip boundary condition.
Thisfreeslipmodelhasthebenefitofbeingaparameter-freemodelofturbulentboundarylayers,whichrequiresno
boundary layer mesh since no boundary layer is resolved. The basic assumptions underlying the model is that the