A staggered space-time discontinuous Galerkin method for the three-dimensional incompressible Navier-Stokes equations on unstructured tetrahedral meshes

9 MB·

by  Maurizio Tavelli

A staggered space-time discontinuous Galerkin method for the three-dimensional
incompressible Navier-Stokes equations on unstructured tetrahedral meshes
MaurizioTavellia,1,MichaelDumbserb,2,
aDepartmentofMathematics,UniversityofTrento,
ViaSommarive14,I-38050Trento,Italy
bLaboratoryofAppliedMathematics,DepartmentofCivil,EnvironmentalandMechanicalEngineering,UniversityofTrento,ViaMesiano77,
I-38123Trento,Italy
6
1
Abstract
0
2
Inthispaperweproposeanovelarbitraryhighorderaccuratesemi-implicitspace-timediscontinuousGalerkinmethod

n for the solution of the three-dimensional incompressible Navier-Stokes equations on staggered unstructured curved
a tetrahedralmeshes. Astypicalforspace-timeDGschemes,thediscretesolutionisrepresentedintermsofspace-time
J
basisfunctions. Thisallowstoachieveveryhighorderofaccuracyalsointime, whichisnoteasytoobtainforthe
7 incompressibleNavier-Stokesequations. Similartostaggeredfinitedifferenceschemes,inourapproachthediscrete
2
pressure is defined on the primary tetrahedral grid, while the discrete velocity is defined on a face-based staggered

] dual grid. While staggered meshes are state of the art in classical finite difference schemes for the incompressible
A Navier-Stokes equations, their use in high order DG schemes still quite rare. A very simple and efficient Picard
N iterationisusedinordertoderiveaspace-timepressurecorrectionalgorithmthatachievesalsohighorderofaccuracy
. intimeandthatavoidsthedirectsolutionofglobalnonlinearsystems. Formalsubstitutionofthediscretemomentum
h
equationonthedualgridintothediscretecontinuityequationontheprimarygridyieldsaverysparsefive-pointblock
t
a systemforthescalarpressure,whichisconvenientlysolvedwithamatrix-freeGMRESalgorithm. Fromnumerical
m
experimentswefindthatthelinearsystemseemstobereasonablywellconditioned,sinceallsimulationsshowninthis
[ papercouldberunwithouttheuseofanypreconditioner,evenuptoveryhighpolynomialdegrees. Forapiecewise

constantpolynomialapproximationintimeandifpressureboundaryconditionsarespecifiedatleastinonepoint,the
1
v resulting system is, in addition, symmetric and positive definite. This allows us to use even faster iterative solvers,
4 liketheconjugategradientmethod.
8 The flexibility and accuracy of high order space-time DG methods on curved unstructured meshes allows to
3
discretize even complex physical domains with very coarse grids in both, space and time. The proposed method is
7
verifiedforapproximationpolynomialsofdegreeuptofourinspaceandtimebysolvingaseriesoftypical3Dtest
0
. problems and by comparing the obtained numerical results with available exact analytical solutions, or with other
1
numericalorexperimentalreferencedata.
0
6 To the knowledge of the authors, this is the first time that a space-time discontinuous Galerkin finite element
1 method is presented for the three-dimensional incompressible Navier-Stokes equations on staggered unstructured
: tetrahedralgrids.
v
i
X Keywords: highorderschemes,space-timediscontinuousGalerkinfiniteelementschemes,staggeredunstructured
meshes,space-timepressurecorrectionalgorithm,incompressibleNavier-Stokesequationsin3D
r
a
1. Introduction
ThenumericalsolutionofthethreedimensionalincompressibleNavier-Stokesequationsrepresentsaveryimpor-
tant and challenging research topic, both from a numerical and from an application point of view. In the literature,
[email protected] (M. Tavelli)
[email protected] (M. Dumbser)
PreprintsubmittedtoJournalofComputationalPhysics January28,2016
there are many different approaches that have been proposed for the solution of the incompressible Navier-Stokes
equations, for example using classical finite difference methods [1, 2, 3, 4] or continuous finite element schemes
[5,6,7,8,9,10,11]. Veryrecently, alsodifferenthighorderdiscontinuousGalerkin(DG)methodshavebeenpre-
sented for the solution of the incompressible and the compressible Navier-Stokes equations. The first DG schemes
that were able to solve the Navier-Stokes equations were those of Bassi and Rebay [12] and Baumann and Oden
[13, 14]. Many other methods have been presented in the meantime, see for example [12, 13, 14, 15, 16, 17, 18,
19, 20, 21, 22, 23, 24, 25, 26, 27] for a non-exhaustive overview of the ongoing research in this very active field.
In most DG schemes, the DG discretization is only used for space discretization, while the time discretization uses
standardexplicitorimplicittimeintegratorsknownforordinarydifferentialequations,followingtheso-calledmethod
oflinesapproach. ThemethodoflineshasalsobeenusedbyCockburnandShuintheirwell-knownseriesofpapers
[28, 29, 30] on DG schemes for time-dependent nonlinear hyperbolic systems. In contrast to the method of lines
approach,thefamilyofspace-timediscontinuousGalerkinfiniteelementschemes,whichwasintroducedforthefirst
timebyvanderVegtetal. in[31,32,33],treatsspaceandtimeinaunifiedmanner. Thisisachievedbyusingtestand
basisfunctionsthatdependonbothspaceandtime,see[34,35,36,37,38,39,40]foranoverviewofrecentresults.
Foraveryearlyimplementationofcontinuousspace-timefiniteelementschemes,thereaderisalsoreferredto[41].
From an application point of view, it is very important to consider the fully three-dimensional Navier-Stokes
equations,inordertocapturetherelevantflowfeaturesthatareobservedinlaboratoryexperiments,see[42,43,44,
45]. Thismeansthattheuseofatwo-dimensionalalgorithmisinmostcasesinappropriatetoreproducetheresults
ofphysicalexperiments,evenforgeometriesthatcanbeconsideredessentiallytwo-dimensional. Theimportanceof
fully three-dimensional computations has been shown, for example, in [45, 46, 47, 48, 49, 50, 51]. Unfortunately,
themeshgenerationforcomplexandrealistic3Dgeometriesisstillnowadaysquitedifficult,andthecomputational
cost of a fully three-dimensional simulation grows very quickly with increasing mesh resolution. In this context,
it becomes crucial to use unstructured simplex meshes, since they help to simplify the process of mesh generation
significantly compared to unstructured hexahedral meshes. Furthermore, it is at the same time also crucial to use
very high order accurate methods in both space and time, since they allow to reduce the total number of elements
significantly,comparedtolowordermethods,whilekeepingatthesametimeahighlevelofaccuracyofthenumerical
solution. SincethesolutionoftheincompressibleNavier-Stokesequationsrequiresnecessarilythesolutionoflarge
systems of algebraic equations, it is indeed very important to derive a scheme that uses a stencil that is as small as
possible,inordertoimprovethesparsitypatternoftheresultingsystemmatrices.Itisalsodesirabletodesignmethods
thatleadtoreasonablywellconditionedsystemsthatcanbesolvedwithiterativesolvers,liketheconjugategradient
method[52]ortheGMRESalgorithm[53].
Forstructuredgrids,numericalschemescanbeusuallyderivedrathereasilyinmultiplespacedimensions,thanks
totheparticularregularityofthemesh. Onthecontrary,thedevelopmentofnumericalschemesongeneralunstruc-
tured meshes in three space dimensions is not as straightforward and requires some care in the derivation and the
implementationofthemethod. ParticulardifficultiesoftheincompressibleNavier-Stokesequationsarisefromtheir
nonlinearityandfromtheellipticnatureofthePoissonequationforthepressure,thatisalsoobtainedonthediscrete
level when substituting the momentum equation into the discrete continuity equation. A unified analysis of several
variantsoftheDGmethodappliedtoanellipticmodelproblemhasbeenprovidedbyArnoldetal. in[54].
While the use of staggered grids is a very common practice in the finite difference community, its use is not so
widespreadinthecontextofhighorderDGschemes. ThefirststaggeredDGschemes,basedonavertex-baseddual
grid,havebeenproposedin[55,56]. OtherrecenthighorderstaggeredDGschemesthatuseanedge-baseddualgrid
havebeenforwardedin[57,58,59]. Theadvantageinusingedge-basedstaggeredgridsisthattheyallowtoimprove
significantly the sparsity pattern of the final linear system that has to be solved for the pressure. Very recently, a
newfamilyofstaggered semi-implicitDGschemesforthesolutionofthetwodimensionalshallowwaterequations
waspresentedbyDumbser&Casulli[59]andTavelli&Dumbser[60]. Subsequently,thesesemi-implicitstaggered
DG schemes have been successfully extended also to the two-dimensional incompressible Navier-Stokes equations
by Tavelli & Dumbser in [61, 39]. One year later, a staggered DG formulation for the 2D incompressible Navier-
Stokes equations has been reproposed independently also in [62]. Alternative semi-implicit discontinuous Galerkin
schemesoncollocatedgridshavebeenpresented,forexample,in[63,64,65,66,67]. Thesesemi-implicitschemes
trytocombinethesimplicityofexplicitmethodsfornonlinearPDEwiththestabilityandefficiencyofimplicittime
discretizations.
Inthispaperweproposeanew,arbitraryhighorderaccuratestaggered space-timediscontinuousGalerkinfinite
2
elementmethodforthesolutionofthethree-dimensionalincompressibleNavier-Stokesequationsoncurvedunstruc-
turedtetrahedralmeshes,followingsomeoftheideasoutlinedin[39]forthetwo-dimensionalcase. Forthatpurpose
we mimic the philosophy of staggered semi-implicit finite difference schemes, such as discussed and analyzed in
[1,2,3,4,68,69,70,71,72,73,74,75,76,77],wherethediscretepressurefieldisdefinedontheprimarygrid,while
thediscretevelocityfieldisdefinedonanedge-basedstaggereddualgrid.
For the staggered space-time DG scheme proposed in this paper, we use a primal mesh composed of (curved)
tetrahedral elements, and a face-based staggered dual mesh that consists of non-standard five-point hexahedral ele-
ments that are obtained by connecting the three nodes of a face of the primal mesh with the barycenters of the two
tetrahedra that share the common face. The face-based dual grid used here corresponds to the choice made also in
[78,79,80,58].Thesespatialelementsarethenextendedtospace-timecontrolvolumesusingasimpletensorproduct
inthetimedirection.
Since all quantities are readily defined where they are needed, the staggered DG scheme does not require the
useofRiemannsolvers(numericalfluxfunctions),apartfromthenonlinearconvectiveterms,whicharetreatedina
standardway. FortheconvectivepartoftheincompressibleNavier-Stokesequations,weuseastandardDGscheme
for hyperbolic PDE on the main grid, based on the local Lax-Friedrichs (Rusanov) flux [81]. For that purpose, the
velocityfieldisfirstinterpolatedfromthedualgridtothemaingrid,assuggestedin[59]. Thisallowsustousethe
same staggered space-time DG scheme again to discretize also the viscous terms, where now the velocity gradient
thatisneededfortheevaluationoftheviscousfluxesiscomputedontheface-basedstaggereddualgrid. Inthisway,
we can avoid again the use of numerical flux functions for the viscous fluxes, and furthermore, the structure of the
resultinglinearsystemsfortheviscoustermsisverysimilartothepressuresystem.
The discrete momentum equation is then inserted into the discrete continuity equation in order to obtain the
discrete form of the pressure Poisson equation. Thanks to the use of a staggered grid, this leads to a very sparse
five-point block system, with the scalar pressure as the only unknown quantity. Note that the same algorithm on a
collocatedgrid wouldproducea17-pointstencil,sinceitwouldalsoinvolveneighborsofneighbors3. Ontheother
hand,ifonedoesnotsubstitutethemomentumequationintothecontinuityequationonacollocatedgrid,onecould
stillobtainafivepointstencil, butwiththepressureand thevelocityvectorasunknowns, hencethefinalsystemto
solveisfourtimeslargerthanthecorrespondingsystemofourstaggeredDGscheme. Itisthereforeveryclearthat
even in the DG context, the use of a staggered mesh is very beneficial, since it allows to produce a linear system
with the smallest possible stencil and with the smallest number of unknowns, compared to similar approaches on a
collocatedmesh.
Oncethenewpressurefieldisknown,thevelocityvectorfieldcansubsequentlybeupdateddirectly.Averysimple
Picarditerationthatembracestheentireschemeineachtimestepisusedinordertoachievearbitraryhighorderof
accuracyintimealsoforthenonlinearconvectiveandviscousterms,withoutintroducinganonlinearityinthesystem
forthepressure.
The rest of the paper is organized as follows: in Section 2 we derive and present the new numerical method
method. Section2.5containsthedetailsaboutthediscretizationofthenonlinearconvectivetermsonthemaingrid,
while the velocity gradients for the viscous terms are discretized again on the face-based staggered dual mesh. In
Section 2.7 we discuss the important special case of a high order DG discretization in space, while using only a
piecewiseconstantpolynomialapproximationintime,leadingtosymmetricpositivedefinitesystemsforthepressure
andtheviscousterms. Finally, inSection3thenewnumericalschemeproposedinthispaperisrunonasetof3D
benchmarkproblems,comparingthenumericalresultseitherwithexistinganalyticalornumericalreferencesolutions,
orwithavailableexperimentalresults. ThepapercloseswithsomeconcludingremarksprovidedinSection4.
3ThediscretecontinuityequationofaDGschemeonacollocatedgridinvolvesthevelocityintheelementitselfandinitsfourneighbor
elements,duetothenumericalfluxontheelementboundaries.Furthermore,inthediscretemomentumequationthevelocityfieldineachtetrahedral
elementdependsonthepressureinthecellitselfandinitsfourneighbors.Insertingnowthemomentumequationintothecontinuityequationon
thediscretelevelinvolvesatotalof1+4+4·3=17elementsforthepressure!
Onastaggeredmeshinstead,thediscretecontinuityequationinvolvesonlythevelocitiesofthefourdualelementsassociatedwiththefacesofthe
primaryelement. Thediscretemomentumequationwrittenontheface-baseddualgridonlyinvolvesthepressureofthetwotetrahedrathatshare
thecommonface.Hence,substitutingthemomentumequationintothecontinuityequationleadstoa1+4=5pointstencilforthepressure,which
involvesonlytheelementanditsfourneighbors.
3
2. Staggeredspace-timeDGschemeforthe3DincompressibleNavier-Stokesequations
2.1. Governingequations
Thethree-dimensionalincompressibleNavier-Stokesequationscanbewrittenas
∂v
+∇·F +∇p=∇·(ν∇v)+S, (1)
∂t c
∇·v=0, (2)
wherex = (x,y,z)isthevectorofspatialcoordinatesandtdenotesthetime; p = P/ρindicatesthenormalizedfluid
pressure;Pisthephysicalpressureandρistheconstantfluiddensity;ν = µ/ρisthekinematicviscositycoefficient;
v = (u,v,w) is the velocity vector; u, v and w are the velocity components in the x, y and z direction, respectively;
S=S(x,t)isavectorofgivensourceterms;F =v⊗visthefluxtensorofthenonlinearconvectiveterms,namely:
c
 
Fc = uvuu uvvv uvww .
wu wv ww
Theviscositytermcanbegroupedwiththenonlinearconvectiveterm,i.e. themomentumEq. (1)thenreads
∂v
+∇·F+∇p=S, (3)
∂t
whereF=F(v,∇v)=F (v)−ν∇visthenonlinearfluxtensorthatdependsonthevelocityanditsgradient.
c
2.2. Staggeredunstructuredmeshandassociatedspace-timebasisfunctions
Throughoutthispaperweuseamaingridthatiscomposedof(eventuallycurved)tetrahedralsimplexelements,
andastaggeredface-baseddualgrid,consistinginnon-standardfive-pointhexahedralelements. Thesespatialcontrol
volumesarethenextendedtospace-timecontrolvolumesusingatensorproductintimedirection. Inthefollowing,
thestaggeredmeshinspaceisdescribedindetailandissubsequentlyalsoextendedtothetimedirection. Themain
notationistakenastheonepresentedforthetwodimensionalmethodproposedin[39]andissummarizedherefor
thethreedimensionalcase.
2.2.1. Staggeredspace-timecontrolvolumes
The spatial computational domain Ω is covered with a set of N non-overlapping tetrahedral elements T with
e i
i = 1...N . Bydenotingwith N thetotalnumberoffaces,the j−thfacewillbecalledΓ . B(Ω)denotesthesetof
e d j
indices jcorrespondingtoboundaryfaces. TheindicesofthefourfacesofeachtetrahedronT constitutethesetS
i i
definedbyS = {j ∈ [1,N ] | Γ isafaceofT}. Forevery j ∈ [1...N ]−B(Ω)thereexisttwotetrahedrathatshare
i d j i d
acommonfaceΓ . Weassignarbitrarilyaleftandarightelement,calledT andT ,respectively. Thestandard
j (cid:96)(j) r(j)
positivedirectionisassumedtobefromlefttoright. Let(cid:126)n denotetheunitnormalvectordefinedonthefacenumber
j
jandthatisorientedwithrespecttothepositivedirectionfromlefttoright. Foreverytetrahedralelementnumberi
andfacenumber j∈S ,theindexoftheneighbortetrahedronthatsharesthecommonfaceΓ isdenotedby℘(i, j).
i j
For every j ∈ [1,N ]−B(Ω) the dual element (a non-standard5-point hexahedron) associated with Γ is called
d j
H anditisdefinedbythetwocentersofgravityofT andT andthethreeverticesofΓ ,seealso[78,80,60].
j (cid:96)(j) r(j) j
WedenotebyT = H ∩T theintersectionelementforeveryiand j ∈ S . Figures1and2summarizethenotation
i,j j i i
usedonthemaintetrahedralmeshandontheassociateddualgrid. Weexdendourdefinitionsonthemaingridtothe
dualone,namely: N isthetotalamountofsidesofH ;Γ indicatesthel-thside;∀j,thesetofsideslof jisindicated
l j l
withS ;∀l,(cid:96) (l)andr (l)aretheleftandtherighthexahedralelement,respectively;(cid:126)n isthestandardnormalvector
j jl jl l
definedonlandassumedpositivewithrespecttothestandardorientationonl(defined,asforthemaingrid,fromthe
lefttotheright).
Inthetimedirectionwecoverthetimeinterval[0,T]withasequenceoftimes0=t0 <t1 <t2...<tN <tN+1 =T.
Wedenotethetimestepbetweentnandtn+1by∆tn+1 =tn+1−tnandtheassociatedtimeintervalbyTn+1 =[tn,tn+1],for
n = 0...N. Inordertoeasethenotation,sometimeswewillsimplywrite∆t = ∆tn+1. Inthiswaythegenericspace-
timeelementdefinedinthetimeinterval[tn,tn+1]isgivenbyTst =T ×Tn+1 forthemaingridand Hst = H ×Tn+1
i i j j
forthedualgrid.
4
L
Γ
Γj3 j1 Hj i
T
i
T
i
Γ
j2
℘(i, j)
Γ
Γ Γj3 j1 Hj i
j4 r(j)
T
i
(cid:126)n
j T
(cid:96)(j) i
Γ
j2
℘(i, j)
L
Γ
j
Γ
j4 r(j)
Figure1:AtetrahedralelementoftheprimarymeshwithSi={j1,j2,j3,j4}(left)andthestandardorientationusedthroughoutthispaper(right). (cid:96)(j) n(cid:126)j
Γ
H j
j
Γ
Γj3℘(i, j) j1 Hj i
i
T
i
Ti Ti
Γ
j2
℘(i, j)
Γ
j4 r(j)
℘(i, j) Γ
Γj3 j1 Hj i Figure2:Anexampleofadualelement(anon-standard5-pointhexahedron,highlightedinblue)associatedwiththefaceΓj. n(cid:126)
j
(cid:96)(j)
T
i
T
i
2.3. Space-timebasisfunctions
Γ
j2 ℘(i, j) Wefirstconstructthespatialbasisfunctionsandthenweextendthemtothetimedirectionusingasimpletensor Γj
product. Fortetrahedralelements, thebasisfunctionsaregeneratedonastandardreferencetetrahedron, definedby
Ti Γj4 Tref ={(ξ,η,ζ)∈R3 | 0≤ξ+η+ζ ≤1}. Wewritert(hje)basisfunctiononthereferenceelementas
(cid:88)p (cid:88)p−r1p−(cid:88)r2−r1
Γj1 φk(ξ,η,(cid:96)ζ()j)=r1n=(cid:126)0j r2=0 r3=0 αkrξr1ηr2ζr3 :=αkrξr, (4)
for some coefficients α and the multi-index r = (r ,r ,r ). We then define N = (p+1)(p+2)(p+3) nodal points ξ =
kr 1 2 3 φ 6 j
(ξ ,η ,ζ ) = (j /p, j /p, j /p), with the multi-index j = (j , j , j ) and 0 ≤ j + j + j ≤ p, as in standard
cojn1forjm2 inj3gfinitee1leme2nts.W3ethenimposetheclassicalintΓerjpol1atio2nc3onditionforn1odal2finite3elementsφ (ξ )=δ ,
k j kj
withtheusualKroneckersymbolδ . ThismeansthatwehavechosenanodalbasiswhichisdefinedbytheLagrange
kj
interpolationpolynomialsthatpassthroughthenodesgivenbythestandardnodesofconformingfiniteelements.This
leadstothelinearsystemα ξr =δ forthecoefficientsα thatcanbesolvedanalyticallyforeverypolynomialdegree
kr j kj kr
5
ponthereferencetetrahedron. InthiswayweobtainN basisfunctionsonT ,{φ } . Theconnectionbetween
φ ref k k∈[1,Nφ]
the reference coordinates ξ and the physical coordinates x is performed by the map T(·,T) = T : T −→ T for
i i i ref
everyi = 1...N anditsinverse,calledT−1(·,T) = T−1 : T ←− T . Themapsfromthephysicaltothereference
e i i i ref
coordinates can be constructed following a classical sub-parametric or a complete iso-parametric approach and in
generalwewillwrite,foralli=1...N ,φ(i)(x,y,z)=φ (T(x,y,z)).
e k k i
Unfortunately,itisnotsoeasytoconstructasimilarnodalbasisonthedualmesh,duetotheuseofnon-standard
5-pointhexahedralelements. Asdiscussedin[82],thedefinitionofbasisfunctionsbasedonLagrangeinterpolation
polynomialsonthiskindofelementisproblematic, sinceforspecialconfigurationsofthevertexcoordinatesofthe
dual elements, the linear system to be solved for the classical interpolation condition of a nodal basis can become
singular. ThisdoesnotallowtheconstructionofanodalpolynomialbasisforagenericelementH andthereforeone
j
hastopasstorationalfunctionsofpolynomialsinsteadofusingsimplepolynomialfunctionsinthatcase.
Therefore, for the basis functions on the dual grid directly we choose a simple Taylor-type modal basis [83]
directly in the physical space, hence the basis functions will consequently depend on the element j ∈ [1,N ]. The
d
basisfunctionsread
(cid:16)x−x(j)(cid:17)k1(cid:16)y−y(j)(cid:17)k2(cid:16)z−z(j)(cid:17)k3
ψ(j)(x)= 0 0 0 , (5)
k hk1+k2+k3
j
wherexj =(x,y,z)(j)isthecenterofthedualelementandh isacharacteristiclengthofH usedforscalingthebasis.
0 0 j j
Here, 0 ≤ k +k +k ≤ p, i.e. weusetheoptimalnumberofpolynomialsofdegree pinthreespacedimensions,
1 2 3
namelyN = N . Withthischoicewegetonlyamodalbasisforthedualhexaxedralelements,i.e. iftheconvective
ψ φ
termisdirectlycomputedonthedualmeshaccordingtothenaturalextensionofthemethodproposedin[39],thenit
hastobecomputedaccordingtoamodalapproach,whichismoreexpensivethananodalone.
Finally,thetimebasisfunctionsareconstructedonareferenceintervalI =[0,1]forpolynomialsofdegree p . In
γ
thiscasetheresulting N = p +1basisfunctions{γ } aredefinedastheLagrangeinterpolationpolynomials
passing through the Gauγss-Leγgendre quadrature poinktsk∈f[o1,rNγt]he unit interval. For every time interval [tn,tn+1], the
map between the reference interval and the physical one is simply given by t = tn +τ∆tn+1 for every τ ∈ [0,1].
Usingthetensorproductwecanfinallyconstructthebasisfunctionsonthespace-timeelementsTst and Hst suchas
i j
φ˜(ξ,η,ζ,τ) = φ(ξ,η,ζ)·γ(τ) and ψ˜(j)(x,y,z,t) = ψ(j)(x,y,z)·γ(τ(t)). The total number of basis functions becomes
Nst = N ·N andNst = N ·N .
φ φ γ ψ ψ γ
2.4. Staggeredsemi-implicitspace-timeDGscheme
Thediscretepressure p isdefinedonthemaingrid,namely p (x,t)| = p(x,t)= p,whilethediscretevelocity
h h Tst i i
vectorfieldv isdefinedonthedualgrid,namelyv (x,t)| =v (x,t)=vi .
h h Hst j j
j
Thenumericalsolutionof(2)-(3)isrepresentedineachspace-timeelementTst andHst betweentimestnandtn+1
i j
bypiecewisepolynomialsas
Nst
(cid:88)φ
p(x,t)= φ˜(i)(x,t)pˆn+1 =:φ˜(i)(x,t)pˆn+1, (6)
i l l,i i
l=1
Nst
(cid:88)ψ
v (x,t)= ψ˜(j)(x,t)vˆn+1 =:ψ˜(j)(x,t)vˆn+1, (7)
j l l,j j
l=1
wherethevectorofbasisfunctionsφ˜(x,t)isgeneratedfromφ˜(ξ,η,ζ,τ)onT ×[0,1]whileψ˜(j)(x,t)isdefinedfor
std
every j∈[1...N ]−B(Ω)directlyinthephysicalspace.
d
A weak formulation of the continuity equation (2) is obtained by multiplying it with a test function φ˜(i) and
k
integratingoverthespace-timecontrolvolumeTst,foreveryk=1...Nst. Theresultingweakformulationreads
i φ
(cid:90)
φ˜(i)∇·vdxdt=0, (8)
k
Tst
i
6
withdx = dxdydz. Similarly, multiplicationofthemomentumequation(3)bythetestfunctionψ˜(j) andintegrating
k
overacontrolvolumeHst yields
j
(cid:90) (cid:32) (cid:33) (cid:90) (cid:90)
∂v
ψ˜(j) +∇·F dxdt+ ψ˜(j)∇p dxdt= ψ˜(j)Sdxdt, (9)
k ∂t k k
Hst Hst Hst
j j j
forevery j=1...N andk=1...Nst. UsingintegrationbypartsEq. (8)reads
d ψ
(cid:73) (cid:90)
φ˜(i)v·(cid:126)n dSdt− ∇φ˜(i)·vdxdt=0, (10)
k i k
∂Tst Tst
i i
where(cid:126)n indicatestheoutwardpointingunitnormalvector. Duetothediscontinuityof p andv ,equations(9)and
i h h
(10)havetobesplitasfollows:
 
(cid:88)j∈SiΓ(cid:90)st φ˜(ki)vj·(cid:126)nijdSdt−T(cid:90)st ∇φ˜(ki)·vjdxdt=0, (11)
j i,j
and
(cid:90) (cid:32) (cid:33) (cid:90) (cid:90) (cid:90) (cid:90)
∂v (cid:16) (cid:17)
ψ˜(j) j +∇·F dxdt+ ψ˜(j)∇p dxdt+ ψ˜(j)∇p dxdt+ ψ˜(j) p −p (cid:126)n dSdt= ψ˜(j)Sdxdt,
k ∂t k (cid:96)(j) k r(j) k r(j) (cid:96)(j) j k
Hst Tst Tst Γst Hst
j (cid:96)(j),j r(j),j j j
(12)
where(cid:126)nij =(cid:126)ni|Γsjt;Tis,tj =Ti,j×Tn+1;andΓsjt =Γj×Tn+1.NotethatthepressurehasadiscontinuityalongΓsjtinsidethe
hexahedralelementHstandhencethepressuregradientin(9)needstobeinterpretedinthesenseofdistributions,asin
j
path-conservativefinitevolumeschemes[84,85].Thisleadstothejumptermspresentin(12),see[39].Alternatively,
the same jump term can be produced also via forward and backward integration by parts, see e.g. the well-known
workofBassiandRebay[12]. Usingdefinitions(6)and(7),werewritetheaboveequationsas
 
(cid:88)j∈SiΓ(cid:90)st φ˜(ki)ψ˜(lj)(cid:126)nijdSdt ·vˆnl,+j1−T(cid:90)st ∇φ˜(ki)ψ˜(lj)dxdt ·vˆnl,+j1=0, (13)
j i,j
and
(cid:90) (cid:90) (cid:90) (cid:90)
∂v
ψ˜(j) j dxdt+ ψ˜(j)∇·Fdxdt+ ψ˜(j)∇φ˜((cid:96)(j))dxdt pˆn+1 + ψ˜(j)∇φ˜(r(j))dxdt pˆn+1
k ∂t k k l l,(cid:96)(j) k l l,r(j)
Hst Hst Tst Tst
j j (cid:96)(j),j r(j),j
(cid:90) (cid:90) (cid:90)
+ ψ˜(j)φ˜(r(j))(cid:126)n dSdtpˆn+1 − ψ˜(j)φ˜((cid:96)(j))(cid:126)n dSdtpˆn+1 = ψ˜(j)Sdxdt,
k l j l,r(j) k l j l,(cid:96)(j) k
Γst Γst Hst
j j j
(14)
wherewehaveusedthestandardsummationconventionfortherepeatedindexl. Integratingthefirstintegralin(14)
bypartsintimeweobtain
(cid:90) ∂v (cid:90) (cid:90) (cid:90) ∂ψ˜(j)
ψ˜(j) j dxdt= ψ˜(j)(x,tn+1)v (x,tn+1)dx− ψ˜(j)(x,tn)v (x,tn)dx− k v (x,t)dxdt. (15)
k ∂t k j k j ∂t j
Hst H H Hst
j j j j
7
In Eq. (15) we can recognize the fluxes between the current space-time element H ×Tn+1, the future and the past
j
space-timeelements,aswellasaninternalcontributionthatconnectsinanasymmetricwaythedegreesoffreedom
inside the element Hst. Note that the asymmetry appears only in the volume contribution in (15). For the spatial
j
integral at time tn we will insert the boundary-extrapolated numerical solution from the previous time step, which
correspondstoupwindingintimedirectionduetothecausalityprinciple. BysubstitutingEq. (15)into(14)andusing
thecausalityprinciple,weobtainthefollowingweakformulationofthemomentumequation:
 
(cid:90) ψ˜k(j)(x,tn+1)ψ˜(lj)(x,tn+1)dx− (cid:90) ∂ψ∂˜t(kj)ψ˜(lj)dxdtvˆnl,+j1−(cid:90) ψ˜(kj)(x,tn)ψ˜(lj)(x,tn)dxvˆnl,j+ (cid:90) ψ˜(kj)∇·Fdx
H Hst H Hst
j j j j
(cid:90) (cid:90) (cid:90) (cid:90)
+ ψ˜(j)∇φ˜((cid:96)(j))dx pˆn+1 + ψ˜(j)∇φ˜(r(j))dx pˆn+1 + ψ˜(j)φ˜(r(j))(cid:126)n dS pˆn+1 − ψ˜(j)φ˜((cid:96)(j))(cid:126)n dS pˆn+1
k l l,(cid:96)(j) k l l,r(j) k l j l,r(j) k l j l,(cid:96)(j)
Tst Tst Γst Γst
(cid:96)(j),j r(j),j j j
(cid:90)
= ψ˜(j)Sdxdt, (16)
k
Hst
j
Foreveryiand j,Eqs. (13)and(16)canbewritteninacompactmatrixformas
(cid:88)
D vˆn+1 =0, (17)
i,j j
j∈Si
and
(cid:16) (cid:17)
M+− M◦ vˆn+1− M−vˆn+Υ (v,∇v)+R pˆn+1−L pˆn+1 =S , (18)
j j j j j j j r(j) j (cid:96)(j) j
respectively,where:
(cid:90)
M+ = ψ˜(j)(x,t(1))ψ˜(j)(x,t(1))dx, (19)
j k l
H
j
(cid:90)
M− = ψ˜(j)(x,t(0))ψ˜(j)(x,t(1))dx, (20)
j k l
H
j
(cid:90) ∂ψ˜(j)
M◦ = k ψ˜(j)dxdt, (21)
j ∂t l
Hst
j
(cid:90)
Υ = ψ˜(j)∇·Fdxdt (22)
j k
Hst
j
(cid:90) (cid:90)
D = φ˜(i)ψ˜(j)(cid:126)n dSdt− ∇φ˜(i)ψ˜(j)dxdt, (23)
i,j k l ij k l
Γst Tst
j i,j
(cid:90) (cid:90)
R = ψ˜(j)φ˜(r(j))(cid:126)n dSdt+ ψ˜(j)∇φ˜(r(j))dxdt, (24)
j k l j k l
Γst Tst
j r(j),j
(cid:90) (cid:90)
L = ψ˜(j)φ˜((cid:96)(j))(cid:126)n dSdt− ψ˜(j)∇φ˜((cid:96)(j))dxdt, (25)
j k l j k l
Γst Tst
j (cid:96)(j),j
8
(cid:90)
S = ψ˜(j)Sdxdt. (26)
j k
Hst
j
Notehow M◦ introduces,for p > 0,anasymmetriccontributionthatwillleadtoanasymmetryofthemainsystem
j γ
for the discrete pressure. The action of matrices L and R can be generalized by introducing the new matrix Q ,
i,j
definedas
(cid:90) (cid:90)
Q = ψ˜(j)∇φ˜(i)dxdt− ψ˜(j)φ˜(i)σ (cid:126)n dsdt, (27)
i,j k l k l i,j j
Tst Γst
i,j j
whereσ isasignfunctiondefinedby
i,j
r(j)−2i+(cid:96)(j)
σ = . (28)
i,j r(j)−(cid:96)(j)
InthiswayQ =−L andQ =R ,andthenEq. (18)becomesintermsofQ
(cid:96)(j),j j r(j),j j
(cid:16) (cid:17)
M+− M◦ vˆn+1− M−vˆn+Υ (v,∇v)+Q pˆn+1+Q pˆn+1 =S , (29)
j j j j j j r(j),j r(j) (cid:96)(j),j (cid:96)(j) j
or,equivalently,
(cid:16) (cid:17)
M+− M◦ vˆn+1− M−vˆn+Υ (v,∇v)+Q pˆn+1+Q pˆn+1 =S . (30)
j j j j j j i,j i ℘(i,j),j ℘(i,j) j
Inordertoeasethenotationwewilluse M = M+−M◦. Hence,thediscreteequations(17)-(18)readasfollows:
j j j
(cid:88)
D vˆn+1 =0, (31)
i,j j
j∈Si
M vˆn+1− M F(cid:99)v +Q pˆn+1+Q pˆn+1 =0, (32)
j j j j r(j),j r(j) (cid:96)(j),j (cid:96)(j)
where F(cid:99)v is an appropriate discretization of the nonlinear convective, viscous and source terms that will be pre-
j
sentedlater. Formalsubstitutionofthediscretevelocityfieldgivenbythemomentumequation(32)intothediscrete
continuityequation(31),seealso[69,59],yields
(cid:88) (cid:88) (cid:88)
D M−1Q pˆn+1+ D M−1Q pˆn+1 = D F(cid:99)v . (33)
i,j j i,j i i,j j ℘(i,j),j ℘(i,j) i,j j
j∈Si j∈Si j∈Si
Eq. (33)aboverepresentsablockfive-pointsystemforthepressure pˆn+1.
i
2.5. Nonlinearconvectiveandviscousterms
We now have to choose a proper discretization for the nonlinear convective and viscous terms. As discussed in
[39]weintroduceasimplePicarditerationtoupdatetheinformationaboutthepressure,butwithoutintroducingany
nonlinearityintothefinalsystemforthepressure. Hence,fork=1,N ,werewritesystem(33)as
pic
(cid:88)D M−1Q pˆn+1,k+1+(cid:88)D M−1Q pˆn+1,k+1 =(cid:88)D F(cid:99)vn+1,k+21. (34)
i,j j i,j i i,j j ℘(i,j),j ℘(i,j) i,j j
j∈Si j∈Si j∈Si
TherightsideofEq. (34)canbecomputedbyusingthevelocityfieldatthePicarditerationkandincludingthe
viscous effect implicitly, using a fractional step procedure detailed later. Once the new pressure field is known, the
velocityvectorfieldatthenewPicarditerationvˆn+1,k+1 canbereadilyupdatedfromthediscretemomentumequation
(32).
Toclosetheproblemitremainstospecifyhowtoconstructthenonlinearconvective-diffusionoperatorF(cid:99)vn+1,k+21.
j
Atthispointonecantrytoextendtheprocedurealreadyusedin[39]to3D.However,inthiscasetherearesomeissues
thathavetobetakenintoaccount. Inparticular,sinceweareusingamodalbasisonthestaggereddualnon-standard
5-pointhexahedralmesh,wecannotusethesimplenodalapproximationforthenonlinearconvectivetermFˆ =F (vˆ)
c c
thatconsistsinatrivialpoint-wiseevaluationofthenonlinearoperatorF . Inspiredbythegoodpropertiesobtained
c
9
by the use of staggered grids, here we propose a new procedure for the computation of the nonlinear convective
and viscous terms. For that purpose, the velocity field is first interpolated from the dual grid to the main grid. The
nonlinearconvectivetermscanthenbeeasilydiscretizedwithastandard(space-time)DGschemeonthemaingrid.
Then, the staggered mesh is used again in order to define the gradient of the velocity on the dual elements, which
allowsustoproduceaverysimpleandsparsesystemforthediscretizationoftheviscousterms.
Animplicitdiscretizationoftheviscoustermsonthedualgridleadstoalinearsystemforeachvelocitycomponent
thatisaseven-pointnonsymmetricblocksystemwhichis, however, wellconditionedsinceitcanbewrittenasaν
perturbationoftheidentitymatrix,seee.g. [39]. Here,wewilldevelopadiscretizationoftheviscoustermsthatleads
onlytoafive-pointblocksystemand,moreimportantly,issymmetricandpositivedefiniteforν > 0and p = 0,but
γ
isstillbetterconditionedalsointhegeneralcase p >0.
γ
Givenadiscrete velocityfieldv onthedual gridinthetimeinterval[tn,tn+1], wecanproject thevelocityfield
h
fromthedualmeshtothemaingrid(denotedbyv¯)viastandardL projection,
2
(cid:88)
v¯n+1 = M−1 M vˆn+1, ∀i∈[1,N ], (35)
i i i,j j e
j∈Si
wherev¯n+1denotethedegreesoffreedomofthevelocityonthemaingridand
i
(cid:90) (cid:90)
M = φ˜(i)φ˜(i)dxdt, M = φ˜(i)ψ˜(j)dxdt. (36)
i k l i,j k l
Tst Tst
i i,j
Theprojectionbackontothedualgridisgivenby
vˆn+1 = M−1(cid:16)M(cid:62) vn+1+ M(cid:62) vn+1(cid:17). (37)
j j (cid:96)(j),j (cid:96)(j) r(j),j r(j)
with
(cid:90)
M = ψ˜(j)ψ˜(j)dxdt. (38)
j k l
Hst
j
Wecanrewritethenonlinearconvectiveandviscouspartofthemomentumequationbyintroducingtheviscous
stresstensorσ=−ν∇vasauxiliaryvariable. Theconvectiveandviscoussubsystemofthemomentumequationthen
reads
∂v
+∇·F +∇·σ = 0,
∂t c
σ = −ν∇v. (39)
With the averaged velocity v¯n+1 = φ˜(i)v¯n+1 defined on the main grid and the viscous stress tensor σn+1 = ψ˜(j)σn+1
i l l,i j l l,j
definedonthedualgrid,weobtainthefollowingweakformulationof(39):
(cid:90) (cid:90) (cid:90) ∂φ˜(i)
φ˜(i)(x,tn+1)v¯n+1dx− φ˜(i)(x,tn)v¯ndx− k v¯n+1dxdt+
k i k i ∂t i
T T Tst
i i i
 
∂T(cid:90)ist φ˜(ki)FRcS(cid:0)v¯−,v¯+(cid:1)·(cid:126)nidSdt−T(cid:90)ist ∇φ˜(ki)·Fc(v¯ni+1)dxdt+(cid:88)j∈SiΓ(cid:90)sjt φ˜(ki)σnj+1·(cid:126)nijdSdt−T(cid:90)is,tj ∇φ˜k(i)·σnj+1dxdt=0,
 
(cid:90) ψ˜(kj)(x,tn+1)σnj+1dx=−ν (cid:90) ψ˜(kj)∇v¯n(cid:96)(+j1)dxdt+ (cid:90) ψ˜(kj)∇v¯nr(+j1) dxdt+(cid:90) ψ˜(kj)(cid:16)v¯nr(+j1) −v¯n(cid:96)(+j1)(cid:17)⊗(cid:126)njdSdt. (40)
Hst Tst Tst Γst
j (cid:96)(j),j r(j),j j
10