Energy stable and momentum conserving hybrid finite element method for the incompressible Navier-Stokes equations

0.4 MB·English

by  Robert Jan Labeur

ENERGY STABLE AND MOMENTUM CONSERVING HYBRID
FINITE ELEMENT METHOD FOR THE INCOMPRESSIBLE
NAVIER–STOKES EQUATIONS
ROBERT JAN LABEUR∗ AND GARTH N. WELLS†
Abstract. A hybridmethod forthe incompressibleNavier–Stokes equations ispresented. The
method inherits the attractive stabilizing mechanism of upwinded discontinuous Galerkin methods
when momentum advection becomes significant, equal-order interpolations can be used for the ve-
locityandpressurefields,andmasscanbeconservedlocally. Usingcontinuous Lagrangemultiplier
spacestoenforcefluxcontinuityacrosscellfacets,thenumberofglobaldegreesoffreedomisthesame
asforacontinuous Galerkinmethodonthesamemesh. Differentfromourearlierinvestigations on
theapproachfortheNavier–Stokesequations,thepressurefieldinthisworkisdiscontinuousacross
cellboundaries. Itisshownthatthisleads toverygoodlocal massconservation and, foranappro-
2 priatechoiceoffiniteelementspaces,momentumconservation. Also,anewformofthemomentum
1 transportterms forthe method isconstructed such that global energy stabilityisguaranteed, even
0 intheabsence ofapoint-wisesolenoidal velocityfield. Massconservation, momentum conservation
2 andglobalenergystabilityareprovedforthetime-continuous case, andforafullydiscretescheme.
Thepresentedanalysisresultsaresupportedbyarangeofnumericalsimulations.
n
a Key words. Finite element method; hybrid finite element methods; incompressible Navier–
J
Stokes equations.

5
1 AMS subject classifications. 65N12,65N30,76D05,76D07.

]
E 1. Introduction. A method that combines attractive features of discontinuous
C and continuous Galerkinfinite element methods for the incompressible Navier Stokes
. equations was presented in Labeur and Wells [1] and further extended to the case of
s
c movingdomainsandfree-surfaceflowsinLabeurandWells [2]. The methodincorpo-
[ ratednaturallytheevaluationofupwindedadvectivefluxesoncellfacets,inthesame

spirit as discontinuous Galerkin methods, thereby stabilizing flows with significant
3
v momentum advection, and it is possible to use equal-order polynomial bases for the
2 velocityand pressurecomponents. However,the number ofglobaldegreesoffreedom
2 onagivenmeshisthesameasforacontinuousGalerkinmethodusingthesamepoly-
7
nomial orders. The issue regardingthe significantly greater number of globaldegrees
3
of freedom for low- to moderate-order discontinuous Galerkin methods compared to
.
2 continuous Galerkin methods is thus circumvented. However, the method in [1] was
1
restricted to continuous pressure fields, and it could not be proven that the method
0
is globally energy stable. These short-comings are addressed in this paper, with a
1
: formulation presented that permits discontinuous pressure fields, is globally energy
v
stable, conserves momentum and has excellent local mass conservation properties.
i
X The key to the methodology that we present for constructing schemes is the
r postulation of cell-wise balances, subject to weakly enforced boundary conditions.
a
The boundary condition to be satisfied (weakly) is provided by a function that lives
on cell facets only. An equation for this extra field is furnished by insisting on weak
continuity of the associated ‘numerical’ flux. The concept of weak enforcement of
flux continuity across cell facets is central in hybridized finite element methods (for
an overview see [3]). A feature of these methods is that functions on cells are linked
∗FacultyofCivilEngineeringandGeosciences,DelftUniversityofTechnology,Stevinweg1,2628
CNDelft,TheNetherlands ([email protected])
†DepartmentofEngineering,UniversityofCambridge,TrumpingtonStreet,CambridgeCB21PZ,
UnitedKingdom([email protected])
1
2 R.J.LABEURANDG.N.WELLS
to functions on neighboring cells only via functions that live on cell facets, and not
directly via the flux terms. Therefore, functions on cells can be eliminated locally
in favor of the functions that live on cell facets only (via static condensation), thus
reducingthenumberofgloballycoupleddegreesoffreedom. Ifthefunctionsenforcing
the continuity of the fluxes, and which live only on cell facts, are discontinuous, then
point-wise continuous fluxes can be obtained for suitably chosen function spaces. In
contrast, in our method we advocate the use of facet functions that are continuous,
which leads to a significant reduction in the number of globally coupled degrees of
freedom, since the local elimination procedure will lead to a global problem of the
samesizeasacorrespondingcontinuousGalerkinmethod. Yet,itisstraightforwardto
demonstratelocalmomentumandmassconservation,intermsofthenumericalfluxes,
as is typical of discontinuous Galerkin methods. Also the stabilizing mechanism of
thefluxformulation,involvingtheadvectiontermsandthepressure-velocitycoupling,
are directly inherited from the discontinuous Galerkin method and lead to favorable
stability properties.
We are not alone in considering methods that draw on both discontinuous and
continuousGalerkinmethods. Hughesetal.[4]developedamethodfortheadvection–
diffusion equation, and the formulation in [1] for the advection-diffusion equation re-
duces to that of Hughes et al. [4] in the advective limit. In the diffusive case, there
is a subtle difference, with the diffusive flux in Hughes et al. [4] being upwinded,
whereas a centered approach is used in Labeur and Wells [1]. Simulations using the
approach for the advection–diffusion equation exhibited very good stability proper-
ties, minimal dissipation and standard convergence rates. For the case of the linear
advection–diffusion-reaction equation, stability (via an inf–sup condition) and con-
vergence at a rate of k+1 in the diffusive limit and k+1/2 in the advective limit
was later proved[5]. In the context ofhybridized methods, Cockburnand co-workers
have published a number of works (e.g. [6, 7]) that share features with the method-
ology that we consider. A hybrid field on cell interfaces is presented in Egger and
Scho¨berl [8] for the advection–diffusion problem. Gu¨zey et al. [9] present a hybrid
continuous-discontinuous finite element method for elliptic problems, coined embed-
ded discontinuous Galerkin method, which is conceptually related to the method in
Labeur and Wells [1].
The method formulated and analyzed in this work is an extension of the method
presented in Labeur and Wells [1] for the advection–diffusion equation and for the
incompressible Navier–Stokesequations. Unlike in our previous efforts [1, 2], we con-
sider here pressure fields that are discontinuous across cell facets. The impact of this
on the local (cell-wise) mass conservation properties of the method is demonstrated.
Another feature that distinguishes the formulation developed in this work from our
earlier work for the Navier–Stokes equations is the use of a skew-symmetric form of
the advective term. The derivation of the skew-symmetric formulation is not trivial
inthe consideredsetting,butitisshownthatitbringsthe advantageofguaranteeing
stabilityintermsofthe totalkinetic energy,evenwhenthe velocityfieldisnotpoint-
wisesolenoidal. Thecombinationofdiscontinuouspressurefieldsandskew-symmetric
advectiontermsleadstoamethodthatforequal-orderbasisfunctionspreservesmass
and momentum and is also stable in terms of the total kinetic energy. The analy-
sis results that we present are supported by a number of computations for both the
Stokes and incompressibleNavier–Stokesequations. The computer code necessaryto
reproduce all examples presented in this work is available in the supporting material
[10] under a GNU public license.
ENERGYSTABLEMETHODFORTHENAVIER–STOKESEQUATIONS 3
The remainder of this work is structured as follows. We first define concretely
the problem of interest, and then develop a semi-discrete finite element formulation.
Somepropertiesofthe semi-discreteproblemarethenanalyzed. Thisisfollowedbya
particularfully-discreteformulation,anditisshownthattheconsideredpropertiesof
thesemi-discreteproblemareinheritedbythefullydiscreteproblem. Thisisfollowed
by numerical examples, after which conclusions are drawn.
2. Incompressible Navier–Stokes equations. We consider a domain of in-
terest Ω ⊂ Rd, where 1 ≤ d ≤ 3 is the spatial dimension. The boundary ∂Ω is
assumedsufficiently smoothandthe outwardunit normalvectoron∂Ω is denotedby
n. TheboundaryispartitionedsuchthatΓ ∩Γ =∂ΩandΓ ∪Γ =∅. The time
D N D N
interval of interest is I =(0,t ].
N
Thenon-dimensionalincompressibleNavier–StokesproblemonΩ×I reads: given
the viscosityν, theforcingtermf: Ω×I →Rd,the momentumfluxh: Γ ×I →Rd
N
andthesolenoidalinitialconditionu0: Ω→Rd,findthevelocityfieldu: Ω×I →Rd
and the pressure field p: Ω×I →R such that
∂u
+∇·σ =f on Ω×I, (2.1)
∂t
σ =pI−2ν∇su+u⊗u on Ω×I, (2.2)
∇·u=0 on Ω×I, (2.3)
u=0 on Γ ×I, (2.4)
D
σn−max(u·n,0)u=h on Γ ×I, (2.5)
N
u(x,0)=u0(x) on Ω, (2.6)
where σ is the momentum flux, I is the identity tensor, ∇su = ∇u+∇uT /2 is
the symmetric gradient, in which [∇u] = ∂u /∂x , and [u⊗u] = u u . The
ij i j ij i j
(cid:0) (cid:1)
Neumann boundary condition has been formulated such that on portions of Γ on
N
which u·n<0 (inflow boundaries) the total momentum flux is prescribed, while on
portions of Γ on which u·n≥0 (outflow boundaries) only the diffusive part of the
N
momentum flux is prescribed.
3. Finite element method. The hybrid finite element method is defined in
this section. The essence of the method is posing all balance equations cell-wise in a
weak sense, with a suitably constructed numerical flux, and complementing the cell-
wisebalancelawsbyaglobalequationenforcingweakcontinuityofthenumericalflux
across cell facets.
3.1. Definitions. WeconsideratriangulationT ofthedomainΩintoopen,non-
overlappingsub-domainsK (cells). The outwardunitnormalvectoronthe boundary
∂K ofeachcellisdenotedbyn. AdjacentcellsshareacommonfacetF,andF = F
is the union of all facets, including the exterior boundary facets. A measure of the
S
size of a cell K is denoted by h . When evaluated on a shared facet, h is used to
K K
imply the average cell size measure of the adjacent cells.
Consider first the vector finite element spaces V and V¯ :
h h
V := v ∈ L2(T) d,v ∈[P (K)]d ∀ K ∈T , (3.1)
h h h k
V¯h :=nv¯h ∈(cid:2)L2(F)(cid:3)d, v¯h ∈[Pk¯(F)]d ∀ F ∈Fo, v¯h =0 on ΓD , (3.2)
where P (K)ndenote(cid:2)s the s(cid:3)pace of Lagrange polynomials on K of ordoer k > 0, and
k
Pk¯(F) denotes the space of Lagrangepolynomials on F of order k¯ ≥0. The space Vh
4 R.J.LABEURANDG.N.WELLS
contains vector-valued functions that are discontinuous across cell boundaries, while
functions in V¯ are defined on cell facets only. Furthermore, functions in V¯ satisfy
h h
the homogeneous Dirichlet boundary condition on Γ . Scalar finite element spaces
D
Q and Q¯ are defined by:
h h
Q := q ∈L2(T),q ∈P (K) ∀ K ∈T , (3.3)
h h h m
Q¯h :=(cid:8)q¯h ∈L2(F),q¯h ∈Pm¯ (F) ∀ F ∈F (cid:9), (3.4)
where the polynomial ord(cid:8)ers m ≥ 0 and m¯ ≥ 0. Mirroring (cid:9)the vector spaces, Q
h
contains functions that are discontinuous across cell facets, while functions in Q¯ are
h
defined on cell facets only.
For algorithmic reasons, it may be advantageous to compute with the finite ele-
ment spaces
V¯⋆ :=V¯ ∩ H1(F) d, (3.5)
h h
Q¯⋆h :=Q¯h∩(cid:2)H1(F),(cid:3) (3.6)
in place of V¯ and Q¯ , respectively, using polynomial orders k¯ ≥ 1 and m¯ ≥ 1. This
h h
will be discussed in Section 6, and all computational results presented in Section 7
will employ facet functions that are continuous.
3.2. Semi-discrete weak local/global balances. We formulate now a semi-
discrete finite element problem by considering what we will refer to as ‘local’ and
‘global’ equations. The ‘local’ equations solve the problem cell-wise in which the
velocity and pressure boundary conditions are provided by auxiliary fields that live
on cell facets only. To determine the fields that live on cell facets, ‘global’ equations
areformulatedbyrequiringweakcontinuityofthemassandmomentumfluxesacross
element interfaces. The methodology behind the construction of the formulation is
elucidated by presenting a collection of Galerkin problems for the various balances,
afterwhichtheconsideredGalerkinfiniteelementproblemiscompletelyandformally
defined.
3.2.1. Local/global continuity equation. A Galerkin approximation of the
incompressibility constraint(2.3) in a cell-wise fashion requiresthat the approximate
velocity u ∈V satisfies
h h
u ·∇q dx− uˆ ·nq ds=0 ∀ q ∈Q , (3.7)
h h h h h h
K ZK K Z∂K
X X
where uˆ is the ‘numerical’ mass flux on ∂K, and is chosen to be
h
βh
uˆ =u − K (p¯ −p )n, (3.8)
h h h h
ν+1
in which p ∈Q and p¯ ∈Q¯ are pressure fields and β >0 is a parameter required
h h h h
for stability when using equal-order basis functions for the velocity components and
pressurefields. When using alowerpolynomialorderfor the pressurefieldrelativeto
thevelocityfielditispossibletouseβ =0[11]. Penalizationofthepressurejumpwas
used by Hughes and Franca [12] to stabilize equal-order methods with discontinuous
pressure for the Stokes equation, and by other authors for discontinuous Galerkin
methods [13, 14]. However, different from Hughes and Franca [12], we add a non-
dimensional unit viscosity to the term in the denominator to permit consideration
ENERGYSTABLEMETHODFORTHENAVIER–STOKESEQUATIONS 5
of the inviscid limit. Central in (3.8) is that the pressure-stabilizing term involves
the difference between p and p¯ , rather than the jump in p across a facet as in
h h h
other works [12, 13, 14]. Equation (3.7) is local in the sense that there is no direct
interactionbetweenp onneighboring cells. This is a key feature of the method with
h
practical implications that will be elaborated upon in Section 6.
Thenumericalmassfluxin(3.8)isnotuniqueoncellfacets;itcantakeondifferent
values on different sides of a facet. This is in contrast with standard discontinuous
Galerkin methods, in which the numerical mass flux is constructed such that it is
uniquelydefinedonfacets. A‘global’continuityequationisnowfurnishedbyinsisting
that the numerical mass flux uˆ be weakly continuous across cell facets, in that it
h
satisfies
uˆ ·nq¯ ds− u¯ ·nq¯ ds=0 ∀ q¯ ∈Q¯ , (3.9)
h h h h h h
K Z∂K Z∂Ω
X
where u¯ ∈V¯ . Note that (3.9) implies that uˆ ·n=u¯ ·n (weakly) on ∂Ω.
h h h h
3.2.2. Local/global momentum balance in conservative form. At time t
and given the forcing term f ∈ L2(T) d, the viscosity ν, the velocity u¯ ∈ V¯ , and
h h
pressures p ∈ Q and p¯ ∈ Q¯ , consider a Galerkin approximation u ∈ V that
h h h h h h
(cid:2) (cid:3)
satisfies the following weak formulation of the momentum balance (2.1):
∂u
h ·v dx− σ : ∇v dx+ σˆ n·v ds
h h h h h
∂t
ZΩ K ZK K Z∂K
X X
+ 2ν(u¯ −u )·∇sv nds= f ·v dx ∀ v ∈V , (3.10)
h h h h h h
K Z∂K ZΩ
X
where the momentum flux σ on cells is given by
h
σ =σ(u ,p )=p I −2ν∇su +u ⊗u , (3.11)
h h h h h h h
and the ‘numerical’ momentum flux σˆ on cell boundaries is given by
h
σˆ =σˆ +σˆ , (3.12)
h a,h d,h
where the advective flux σˆ is
a,h
σˆ =σˆ (u ,u¯ ,p ,p¯ )=u ⊗uˆ +(u¯ −u )⊗λuˆ , (3.13)
a,h a h h h h h h h h h
in which uˆ is given by (3.8), λ is a function that takes on a value of either one or
h
zero and is defined below, and the diffusive flux σˆ is
d,h
α
σˆ =σˆ (u ,u¯ ,p¯ )=p¯ I−2ν∇su − 2ν(u¯ −u )⊗n. (3.14)
d,h d h h h h h h h
h
K
Thefunctionλtakesonavalueofoneoninflowcellboundaries(whereuˆ·n<0),and
takesonavalueofzeroonoutflowcellboundaries(whereuˆ·n≥0). Theformulation
for the advective interface flux involves upwinding since σˆ = u⊗uˆ on outflow cell
a
boundaries and σˆ = u¯ ⊗uˆ on inflow cell boundaries. In (3.10), the fourth term on
a
the left-hand side ensures symmetry of the diffusion operator (see Arnold et al. [15]).
In (3.14), α is a penalty parameter, and such a term is typical of interior penalty
methods. Just as for standard interior penalty methods, the role of the penalty term
6 R.J.LABEURANDG.N.WELLS
inthis contextis to ensurestability,asdetailedinWells [5]. Equation(3.10)is ‘local’
in the sense that the weak momentum balance equation is posed cell-wise.
Aswiththenumericalmassflux(3.8),thenumericalmomentumflux(3.12)isnot
single-valuedoncellfacets. A‘global’momentumbalanceequationisthereforeformu-
lated by insisting on continuity of the numerical flux σˆ across element facets. This
h
continuity constraint is imposed weakly by requiring that, for a given flux boundary
condition h∈ L2(Γ ) d,
N
(cid:2) (cid:3)
σˆ n·v¯ ds= (1−λ)(u¯ ⊗u¯ )n·v¯ ds
h h h h h
K Z∂K ZΓN
X
+ h·v¯ ds ∀ v¯ ∈V¯ . (3.15)
h h h
ZΓN
Theaboveequationimpliesthatthe numericalmomentumfluxσˆ nandthemomen-
h
tum flux on Γ , given by h+(1−λ)(u¯ ⊗u¯ )n, coincide in a weak sense.
N h h
3.2.3. Local/global momentum balance in advective form. We now re-
phrase the conservative forms of the local and global momentum balance equations
into advective formats with a view to formulating a skew-symmetric version of the
advective terms.
Considering first the local momentum equation, substitution of the fluxes (3.11),
(3.12) and (3.13) into (3.10) yields
∂u
h ·v dx− σ : ∇v dx− (u ⊗u ): ∇v dx
h d,h h h h h
∂t
ZΩ K ZK K ZK
X X
+ σˆ n·v ds+ (uˆ ·n) u ·v ds
d,h h h h h
K Z∂K K Z∂K
X X
+ λ(uˆ ·n) (u¯ −u )·v ds
h h h h
K Z∂K
X
+ 2ν(u¯ −u )·∇sv nds= f ·v dx, (3.16)
h h h h
K Z∂K ZΩ
X
in whichσ =p I−2ν∇su is the diffusive flux oncells. Applying partialintegra-
d,h h h
tion to the advective terms on K in (3.16),
∂u
h ·v dx− σ : ∇v dx+ (∇u u )·v dx
h d,h h h h h
∂t
ZΩ K ZK K ZK
X X
+ (∇·u )u ·v dx+ ((uˆ −u )·n) u ·v ds
h h h h h h h
K ZK K Z∂K
X X
+ σˆ n·v ds+ λ(uˆ ·n) (u¯ −u )·v ds
d,h h h h h h
K Z∂K K Z∂K
X X
+ 2ν(u¯ −u )·∇sv nds= f ·v dx. (3.17)
h h h h
K Z∂K ZΩ
X
Wechoosetodiscardtheintegralsinvolving∇·u and(uˆ −u )·n,whichbyvirtue
h h h
of the continuity equation (3.7) and under an appropriate regularity assumption on
ENERGYSTABLEMETHODFORTHENAVIER–STOKESEQUATIONS 7
the exact solution will not disturb consistency of a Galerkin scheme (this will be
addressed formally in Section 4). A reduced version of (3.17) now reads:
∂u
h ·v dx+ (∇u u )·v dx− σ : ∇v dx
h h h h d,h h
∂t
ZΩ K ZK K ZK
X X
+ λ(uˆ ·n) (u¯ −u )·v ds+ σˆ n·v ds
h h h h d,h h
K Z∂K K Z∂K
X X
+ 2ν(u¯ −u )·∇sv nds= f ·v dx. (3.18)
h h h h
K Z∂K ZΩ
X
Consideringnextthe globalmomentumbalance,inserting the expressionsfor the
numericalflux (3.12) and(3.13)into the globalflux continuityequation(3.15)yields,
σˆ n·v¯ ds+ (uˆ ·n)u ·v¯ ds
d,h h h h h
K Z∂K K Z∂K
X X
+ λ(uˆ ·n)(u¯ −u )·v¯ ds
h h h h
K Z∂K
X
− (1−λ)(u¯ ·n)u¯ ·v¯ ds= h·v¯ ds. (3.19)
h h h h
ZΓN ZΓN
The second integral in (3.19) can be expanded as
(uˆ ·n)u ·v¯ ds= (uˆ ·n)(u −u¯ )·v¯ ds
h h h h h h h
K Z∂K K Z∂K
X X
+ ((uˆ −u¯ )·n)u¯ ·v¯ ds+ (u¯ ·n)u¯ ·v¯ ds, (3.20)
h h h h h h h
K Z∂K Z∂Ω
X
where we have used that u¯ is single-valued on cell facets, by definition. Discarding
h
the term involving (uˆ −u¯ )·n, which is consistent with continuity equation (3.9),
h h
and substituting (3.20) into (3.19) leads to the following advective form of the global
momentum equation:
σˆ n·v¯ ds− (1−λ)(uˆ ·n)(u¯ −u )·v¯ ds
d,h h h h h h
K Z∂K K Z∂K
X X
+ λ(u¯ ·n)u¯ ·v¯ ds= h·v¯ ds, (3.21)
h h h h
ZΓN ZΓN
where we have taken into account that v¯ =0 on ∂Ω\Γ .
h N
3.3. Semi-discrete finite element formulation. We define now a collection
of functionals that together will define a complete finite element problem. For conve-
nience, the notation U :=(u,u¯,p,p¯) will be used.
Based on the local continuity equation (3.7), we define the functional:
Fc(U;q):= u·∇qdx− uˆ ·nqds, (3.22)
K ZK K Z∂K
X X
8 R.J.LABEURANDG.N.WELLS
where Fc is linear in q. Similarly, from the global continuity equation (3.9), the
functional
F¯c(U;q¯):= uˆ ·nq¯ds− u¯ ·nq¯ds, (3.23)
K Z∂K Z∂Ω
X
is defined, where F¯c is linear in q¯. For the momentum equations, we define a local
momentumbalancefunctionalthatisaweightedcombinationofthelocalconservative
balance (3.16) and the local advective balance (3.18):
∂u
Fm(U;v):= ·vdx−χ (u⊗u): ∇vdx
∂t
ZΩ K ZK
X
+(1−χ) (∇uu)·vdx+χ (uˆ ·n) u·vds
K ZK K Z∂K
X X
+ λ(uˆ ·n) (u¯ −u)·vds− σ : ∇vdx
d
K Z∂K K ZK
X X
+ σˆ n·vds+ 2ν(u¯ −u)·∇svnds− f ·vdx, (3.24)
d
K Z∂K K Z∂K ZΩ
X X
where χ ∈ [0,1] and Fm is linear in v. In the same fashion, the global momentum
flux continuity equations (3.19) and (3.21) are weighted and summed, leading to,
F¯m(U;v¯):=χ (uˆ ·n)u·v¯ds−(1−χ) (uˆ ·n)(u¯ −u)·v¯ds
K Z∂K K Z∂K
X X
+ λ(uˆ ·n)(u¯ −u)·v¯ds+ σˆ n·v¯ds
d
K Z∂K K Z∂K
X X
− (χ−λ)(u¯ ·n)u¯ ·v¯ds− h·v¯ds, (3.25)
ZΓN ZΓN
where F¯m is linear in v¯.
Defining now
F(U;W):=Fm(U;v)+F¯m(U;v¯)+Fc(U;q)+F¯c(U;q¯), (3.26)
where W = (v,v¯,q,q¯), a semi-discrete finite element problem at time t involves:
giventhe forcingterm f ∈ L2(Ω) d the boundary conditionh∈ L2(Γ ) d and the
N
viscosity ν, find U ∈V ×V¯ ×Q ×Q¯ such that
h h h h h
(cid:2) (cid:3) (cid:2) (cid:3)
F(U ;W )=0 ∀ W ∈V ×V¯ ×Q ×Q¯ . (3.27)
h h h h h h h
This completes the formulation of the semi-discrete finite element problem.
4. Properties of the semi-discrete formulation. We now demonstrate the
consistency, mass conservation,momentum conservationand energy stability proper-
ties of the method for the semi-discrete formulation in (3.27). The presented results
hold for the spaces V¯ and Q¯ , and deliberately also for the more restrictive case
h h
in which V¯ and Q¯ are replaced by V¯⋆ and Q¯⋆, respectively, which we advocate in
h h h h
practice and will use in numerical examples.
Proposition 4.1 (consistency). If at a given time t, u ∈ H2(Ω) d and p ∈
H1(Ω) solve equations (2.1)–(2.5), and u¯ = γ(u) and p¯= γ(p) on F, where γ is a
(cid:0) (cid:1)
ENERGYSTABLEMETHODFORTHENAVIER–STOKESEQUATIONS 9
trace operator, then
F(U;W )=0 ∀ W ∈V ×V¯ ×Q ×Q¯ , (4.1)
h h h h h h
for any χ∈[0,1].
Proof. Consideringfirstv =0,v¯ =0andq¯ =0,applyingintegrationbyparts
h h h
to (4.1) leads to
βh
(∇·u)q dx− K (p¯−p) q ds=0 ∀ q ∈Q , (4.2)
h h h h
1+ν
K ZK K Z∂K
X X
which holds due to satisfaction of (2.3) and p¯= γ(p). Setting v = 0, v¯ = 0 and
h h
q =0 in (4.1),
h
uˆ ·nq¯ ds+ uˆ ·nq¯ ds− u¯ ·nq¯ ds=0 ∀ q¯ ∈Q¯ , (4.3)
h h h h h
K Z∂K\∂Ω Z∂Ω Z∂Ω
X
which holds due to the regularity of u and because uˆ =γ(u).
Setting v¯ =0, q =0 and q¯ =0 in (4.1) and applying integration by parts,
h h h
∂u
+∇·σ−f ·v dx−(1−χ) (∇·u)u·v dx
h h
∂t
ZΩ(cid:18) (cid:19) ZΩ
+χ ((uˆ −u)·n)u·v ds+ λ(uˆ ·n)(u¯ −u)·v ds
h h
K Z∂K K Z∂K
X X
α
+ (p¯−p)n·v ds− 2ν(u¯ −u)⊗n·v ds
h h
h
K Z∂K K Z∂K K
X X
+ 2ν(u¯ −u)·∇sv nds=0 ∀ v ∈V , (4.4)
h h h
K Z∂K
X
which holds due to u and p satisfying equations (2.1) and (2.3), the regularity of u
and because p¯=γ(p) and u¯ =uˆ =γ(u). Finally, setting v =0, q =0 and q¯ =0
h h h
in (4.1),
α
σn·v¯ ds+ (p¯−p)n·v¯ ds− 2ν(u¯ −u)⊗n·v¯ ds
h h h
h
K Z∂K K Z∂K K Z∂K K
X X X
+ ((uˆ −u)·n)u·v¯ ds+ λ(uˆ ·n)(u¯ −u)·v¯ ds
h h
K Z∂K K Z∂K
X X
−(1−χ) (uˆ ·n)u¯ ·v¯ ds− (χ−λ)(u¯ ·n)u¯ ·v¯ ds
h h
K Z∂K ZΓN
X
= h·v¯ ds ∀ v¯ ∈V¯ , (4.5)
h h
ZΓN
which holds due to the regularity of u, because p¯=γ(p) and u¯ =uˆ =γ(u) and due
to satisfaction of the flux boundary condition (2.5). Equation (4.1) follows from the
summation of (4.2)–(4.5) and the linearity of F in v, v¯, q and q¯.
Keyto the proofofProposition4.1isthe consistentformulationofthe numerical
fluxes, that is σˆ =σ and uˆ =u if u¯ =u and p¯=p.
Proposition 4.2 (mass conservation). If u and u¯ satisfy (3.27), then
h h
uˆ ·nds=0 ∀ K ∈T (4.6)
h
Z∂K
10 R.J.LABEURANDG.N.WELLS
and
u¯ ·nds=0. (4.7)
h
Z∂Ω
Proof. Setting v = v¯ = 0, q¯ = 1, and q = 1 on the cell K and q = 0 on
h h h h h
T \K leads to (4.6). Setting v =v¯ =0 and q =q¯ =1 in (3.27) leads to (4.7).
h h h h
The local conservation property is in terms of the numerical mass flux uˆ, as is
typical for discontinuous Galerkin methods applied to Stokes flow [13, 14]. Classical
local mass conservation would be satisfied if β = 0, but β must be greater than zero
for stability of the equal-order case [12]. If the pressure field is chosen to be one
polynomial order lower than the velocity field, then β can be set equal to zero and
mass is conserved locally and exactly. However, it will be shown that reducing the
size of the pressure space requires compromising on either momentum conservation
or energy stability.
Proposition 4.3 (momentumconservation). If u and u¯ solve (3.27), and the
h h
function spaces V , V¯ , Q and Q¯ are selected such that for a constant but otherwise
h h h h
arbitrary vector c it holds that v ·c∈Q ∀ v ∈V and v¯ ·c∈Q¯ ∀ v¯ ∈V¯ , then
h h h h h h h h
d
u dx= fdx− σˆ nds ∀ K ∈T, (4.8)
h h
dt
ZK ZK Z∂K
and if Γ =∅
D
d
u dx= fdx− (1−λ)(u¯ ·n)u¯ ds− hds. (4.9)
h h h
dt
ZΩ ZΩ Z∂Ω Z∂Ω
Proof. Setting v = e and q = −(1−χ)u ·e on K, where e is a canonical
h j h h j j
unit basis vector, v =0 and q =0 on T\K, v¯ =0 and q¯ =0 in (3.27),
h h h h
d
u ·e dx+ (uˆ ·n)u ·e ds+ λ(uˆ ·n)(u¯ −u )·e ds
h j h h j h h h j
dt
ZK Z∂K Z∂K
+ σˆ n·e ds= f ·e dx, (4.10)
d,h j j
Z∂K ZK
which from the definition of the fluxes in equation (3.13) proves (4.8).
Setting v = e , v¯ = −e , q = −(1−χ)u ·e and q¯ = −(1−χ)u¯ ·e
h j h j h h j h h j
in (3.27) leads to (4.9) directly.
Local momentum conservation is in terms the numerical flux σˆ , as is a typical
h
feature of discontinuous Galerkin methods. Note also the requirement on the size of
the pressure space relative to the components of the velocity space, which would not
be satisfied by methods that use lower-order polynomials for the pressure than for
the velocity,such as Taylor–Hoodelements. Such elements only conserve momentum
when conservative forms of the momentum equation are used which, in the advective
limit, requires compromising on energy stability. Provided that the requirements on
the sizes of the function spaces are met, momentum conservation holds irrespective
of the value of χ, i.e. for advective as well as conservative forms of the advection
operator,and it will be shown that for χ=1/2 the method is also energy stable (see
proposition 4.4).
ForcaseswithDirichletboundaryconditions(Γ 6=∅),demonstrationofmomen-
D
tum conservation is less straightforward since v¯ can not be set equal to e on Γ .
h j D