# Local conservation laws of continuous Galerkin method for the incompressible
Navier–Stokes equations in EMAC form

Maxim A. Olshanskii Department of Mathematics, University of Houston, Houston,
TX, 77204<EMAIL_ADDRESS>Leo G. Rebholz School of Mathematical and
Statistical Sciences, Clemson University, Clemson, SC, 29364,
<EMAIL_ADDRESS>
###### Abstract

We consider local balances of momentum and angular momentum for the
incompressible Navier-Stokes equations. First, we formulate new weak forms of
the physical balances (conservation laws) of these quantities, and prove they
are equivalent to the usual conservation law formulations. We then show that
continuous Galerkin discretizations of the Navier-Stokes equations using the
EMAC form of the nonlinearity preserve discrete analogues of the weak form
conservation laws, both in the Eulerian formulation and the Lagrangian
formulation (which are not equivalent after discretizations). Numerical tests
illustrate the new theory.

##  1 Introduction

We are interested in conservation properties of continuous Galerkin
discretizations of the incompressible Navier-Stokes equations (NSE), which are
given by

$\left\\{\begin{aligned} \frac{\partial\mathbf{u}}{\partial
t}+(\mathbf{u}\cdot\nabla)\mathbf{u}-\operatorname{div}\boldsymbol{\sigma}&=0\\\
\operatorname{div}\mathbf{u}&=0\end{aligned}\right.\quad\text{in}~{}~{}\Omega.$
(1)

Here, $\boldsymbol{\sigma}$ is the Cauchy stress tensor, and we restrict to
the case of a Newtonian fluid with
$\boldsymbol{\sigma}=2\nu\mathbf{D}(\mathbf{u})-p\mathbf{I}$, where
$\mathbf{D}(\mathbf{u})=\tfrac{1}{2}(\nabla\mathbf{u}+\nabla^{T}\mathbf{u})$
is a rate of deformation tensor.

It is well known that the smooth solution to (1) obeys an array of
conservation laws, including the conservation of momentum, energy, vorticity,
etc., which can be expressed in terms of proper balances for material volumes
of fluid. The development of numerical methods that provide discrete
counterparts for possibly many of these conservation laws is a long-standing
challenge for the computational fluid dynamics community. This challenge has
been addressed by numerous authors and from various perspectives; for example,
see [2, 1, 18, 19, 7, 26, 21, 3, 6] and references therein. Many of these
studies have considered the _global_ conservation properties of numerical
methods, i.e., balances of physical quantities across the entire computational
domain. While properly calibrating these global integral statistics is
necessary for a method to be long-time accurate, it is difficult to see how
this alone can guarantee the quality of a numerical solution.

The proper _local_ balances of momentum, energy, vorticity, etc. represent a
significantly stronger requirement for a numerical solution. Note that
“element-wise conservation” is a common argument used to motivate the
application of discontinuous Galerkin or finite volume discretization
techniques (see, for instance, [5, 17]). At the same time, there is a
widespread belief that continuous (velocity $H^{1}$-conforming) Galerkin
methods inevitably violate local conservation laws; however, see [10, 11] for
a different viewpoint.

Another obstacle in achieving proper discrete counterparts of both local and
global conservation laws for $\mathbf{u}$ is the fact that continuous Galerkin
discretizations of (1) (e.g., conforming finite element methods) typically
enforce the divergence-free constraint only weakly [3]. The purpose of this
paper is to demonstrate that a continuous Galerkin solution, which is only
weakly divergence-free, for (1) does satisfy properly formulated local
conservation laws for momentum and angular momentum when one applies the so-
called EMAC (Energy, Momentum, and Angular Momentum Conserving) formulation of
the NSE.

The EMAC formulation of the discrete NSE was originally developed in [3]. It
re-writes the inertia terms as

$\mathbf{u}\cdot\nabla\mathbf{u}\rightarrow
2\mathbf{D}(\mathbf{u}_{h})\mathbf{u}_{h}+(\operatorname{div}\mathbf{u}_{h})\mathbf{u}_{h},$

along with an altered pressure $p_{h}$ representing
$p-\frac{1}{2}|\mathbf{u}|^{2}$. The motivation for EMAC was that Galerkin
schemes using it can be shown to conserve global energy, momentum and angular
momentum balances when $\operatorname{div}\mathbf{u}_{h}\neq 0$, while schemes
using the common nonlinearity formulations such as convective (CONV:
$\mathbf{u}_{h}\cdot\nabla\mathbf{u}_{h}$), skew-symmetric (SKEW:
$\mathbf{u}_{h}\cdot\nabla\mathbf{u}_{h}+\frac{1}{2}\operatorname{div}\mathbf{u}_{h}$)
and rotational (ROT: $(\nabla\times\mathbf{u}_{h})\times\mathbf{u}_{h}$) do
not preserve some or all of these quantities. Perhaps not surprisingly, the
use of EMAC has become popular for large scale fluid computations in a wide
variety of applications and is shown to give better accuracy especially over
longer time intervals e.g. [22, 15, 23, 8, 16, 14, 27, 3, 4, 20, 12] and is
built into Alya which is a massively parallel multiphysics unstructured finite
element code [28]. In addition to the better discrete physics of EMAC
discussed above, it was proven in [20] that schemes using EMAC are more long
time stable because the Gronwall constant can be shown to be independent of
(explicit) dependence on the Reynolds number; such a result does not hold for
skew-symmetric, convective, or rotation forms for commonly used velocity-
pressure finite elements such as Taylor–Hood elements.

The purpose of this paper is to provide more theoretical justification that
EMAC is superior compared to other discrete nonlinearity formulations, by
proving that continuous Galerkin discretizations using EMAC admit an exact
local balance of momentum and angular momentum. There are very few results for
local conservation properties of continuous finite element methods, with [10,
11] being two fundamental works in this direction. The paper [11] showed that
for NSE, typical Galerkin schemes are not generally conservative, although
this can be ‘fixed’ by multiscale formulation and adding a residual term. One
observation made in this paper is that although local balances written in
different forms – standard Eulerian, Lagrangian, or weak Eulerian and
Lagrangian forms introduced here – represent the same conservation laws of
fluid momentum and angular momentum, after discretization each form can be
different. By considering the weak forms, which we refer to as _diffuse-volume
forms_ , of conservation laws, we can demonstrate that EMAC continuous
Galerkin discretizations exactly preserve properly formulated local momentum
and angular momentum balances. Furthermore, the discrete balances established
here serve as direct analogies to the balances at the partial differential
equation (PDE) level, obviating the need for a multiscale approach and
additional residual terms to establish this connection. We note also that from
the proof construction, it is not possible for SKEW, CONV or ROT to preserve
these local balances of momentum and angular momentum in the same manner that
EMAC does, since they do not preserve them globally.

The rest of the paper is arranged as follows. Section 2 recalls local
conservation laws of momentum and angular momentum. The laws can be
equivalently formulated in Eulerian and Lagrangian forms. Section 3 introduces
a different way to formulate local conservation laws, which is given the name
diffuse-volume form of the conservation laws due to some similarity with
diffuse-interface or phase-field methods in fluid mechanics. We show that this
is another equivalent way to formulate the local balances. Here we also
distinguish between diffuse-volume Eulerian and Lagrangian forms. Section 4
demonstrates how the continuous Galerkin method for the NSE in EMAC form
satisfies discrete counterparts of the (local) diffuse-volume Eulerian and
Lagrangian conservation laws. Section 5 offers a few illustrative numerical
examples.

##  2 Eulerian and Lagrangian forms of momentum and angular momentum
conservation

To formulate local conservation laws satisfied by a _smooth_ solution to (1),
we fix some $t$ and let $\omega\subset\Omega$ be a _fixed_ subdomain of
$\Omega$ with sufficiently smooth boundary $\partial\omega$. For this volume
$\omega$, the balance of momentum and angular momentum take the form:

Moment. $\displaystyle\frac{d}{dt}\int_{\omega}\mathbf{u}\,dx$
$\displaystyle=$
$\displaystyle~{}2\nu\int_{\partial\omega}\mathbf{D}(\mathbf{u})\mathbf{n}\,ds$
$\displaystyle-$ $\displaystyle\quad\int_{\partial\omega}p\mathbf{n}\,ds$
$\displaystyle-$
$\displaystyle\int_{\partial\omega}\mathbf{u}(\mathbf{u}\cdot\mathbf{n})\,ds,$
(2) Ang. Moment.
$\displaystyle\underbrace{\frac{d}{dt}\int_{\omega}\mathbf{u}\times\mathbf{x}\,dx}_{\scriptsize\begin{array}[]{c}\text{momentum}\\\
\text{rate of change}\end{array}}$ $\displaystyle=$
$\displaystyle~{}\underbrace{2\nu\int_{\partial\omega}(\mathbf{D}(\mathbf{u})\mathbf{n})\times\mathbf{x}\,ds}_{\scriptsize\begin{array}[]{c}\text{momentum
change due}\\\ \text{to friction on }\partial\omega\end{array}}$
$\displaystyle-$
$\displaystyle\underbrace{\int_{\partial\omega}p(\mathbf{n}\times\mathbf{x})\,ds}_{\scriptsize\begin{array}[]{c}\text{moment.
change due}\\\ \text{to pressure on }\partial\omega\end{array}}$
$\displaystyle-$
$\displaystyle~{}\underbrace{\int_{\partial\omega}(\mathbf{u}\times\mathbf{x})(\mathbf{u}\cdot\mathbf{n})\,ds.}_{\scriptsize\begin{array}[]{c}\text{the
flux of}\\\ \text{momentum through }\partial\omega\end{array}}$ (11)

Hereafter $\mathbf{n}$ is the outward normal vector on $\partial\omega$.

Local balances (2)–(11) can be interpreted as _Eulerian form_ of the
conservation laws, in contrast to the _Lagrangian form_ formulated for a
material volume below.

We now let $\Omega_{t}\subset\Omega$ be a _material_ volume of the fluid. For
the material volume, the conservation laws for momentum and angular momentum
take the form:

Momentum $\displaystyle\frac{d}{dt}\int_{\Omega_{t}}\mathbf{u}\,dx$
$\displaystyle=$
$\displaystyle~{}2\nu\int_{\partial\Omega_{t}}\mathbf{D}(\mathbf{u})\mathbf{n}\,ds-\int_{\partial\Omega_{t}}p\mathbf{n}\,ds,$
(12) Angular Momentum
$\displaystyle\frac{d}{dt}\int_{\Omega_{t}}\mathbf{u}\times\mathbf{x}\,dx$
$\displaystyle=$
$\displaystyle~{}2\nu\int_{\partial\Omega_{t}}(\mathbf{D}(\mathbf{u})\mathbf{n})\times\mathbf{x}\,ds-\int_{\partial\Omega_{t}}p(\mathbf{n}\times\mathbf{x})\,ds.$
(13)

Of course, for smooth solutions to (1) the Eulerian and Lagrangian forms are
just two different formulations of the same fundamental laws of continuum
mechanics. They both follow from (1), and conversely, together with mass
conservation they imply (1). This equivalence of (1) to the validity of local
conservation laws (specifically, those concerning mass and momentum) is
textbook material. The standard tools used to verify this equivalence include
the divergence theorem, the freedom to choose fluid volumes $\Omega_{t}$ or
$\omega$, and the Reynolds’ transport theorem to handle the Lagrangian form,
which states

$\frac{d}{dt}\int_{\Omega_{t}}f\,dx=\int_{\Omega_{t}}\Big{(}\frac{\partial
f}{\partial t}+\operatorname{div}(f\mathbf{u})\Big{)}\,dx,$

for a smooth scalar function $f$.

Continuous Galerkin methods like the $H^{1}$-conforming finite element method
(FEM) employ finite dimensional subspaces of Sobolev spaces to project (1) and
typically do not offer enough flexibility to verify a direct analogue of
(2)–(11) or (12)–(13). Below we reformulate local conservation laws in a form
more convenient for continuous Galerkin methods.

##  3 Weak form of the conservation laws

The purpose of this section is to derive the conservation laws (2)–(11) and
(12)–(13) in a form more appropriate for a variational formulation. Let
$\omega\subset\Omega$ be an arbitrary subdomain of $\Omega$ with sufficiently
smooth $\partial\omega$. Denote by $\phi$ an _arbitrary_ smooth function such
that $\omega=\mbox{supp}(\phi)$, and set

$\mathbf{\tilde{n}}:=-\frac{\nabla\phi}{|\nabla\phi|}$

for $\mathbf{x}$ such that $\nabla\phi(\mathbf{x})\neq 0$, and let
$\mathbf{\tilde{n}}(\mathbf{x})$ be an arbitrary vector of unit length if
$\nabla\phi(\mathbf{x})=0$. Note that
$\mathbf{\tilde{n}}(\mathbf{x})=\mathbf{n}(\mathbf{x})$ for
$\mathbf{x}\in\partial\omega$. To obtain the weak form of the laws, we
multiply the first equation in (1) by $\phi\mathbf{e}_{i}$ for momentum
conservation and by $\phi\mathbf{e}_{i}\times\mathbf{x}$ for angular momentum
conservation. Doing this for $i=1,\dots,d$, integrating over $\omega$ and by
parts leads after straightforward computations to the following _weak form of
the conservation laws_ :

Moment. $\displaystyle\frac{d}{dt}\int_{\omega}\phi\mathbf{u}\,dx\hskip
17.22217pt=2\nu\int_{\omega}\mathbf{D}(\mathbf{u})\mathbf{\tilde{n}}|\nabla\phi|\,dx-\int_{\omega}p\mathbf{\tilde{n}}|\nabla\phi|\,dx-\int_{\omega}\mathbf{u}(\mathbf{u}\cdot\mathbf{\tilde{n}})|\nabla\phi|\,dx,$
(14) Ang. Moment.
$\displaystyle\frac{d}{dt}\int_{\omega}\phi\mathbf{u}\times\mathbf{x}\,dx=2\nu\int_{\omega}(\mathbf{D}(\mathbf{u})\mathbf{\tilde{n}})\times\mathbf{x}\,|\nabla\phi|\,dx-\int_{\omega}p(\mathbf{\tilde{n}}\times\mathbf{x})|\nabla\phi|\,dx-\int_{\omega}(\mathbf{u}\times\mathbf{x})(\mathbf{u}\cdot\mathbf{\tilde{n}})|\nabla\phi|\,dx.$
(15)

We note that for all calculations below to make sense it is sufficient to
assume $\phi\in W^{1,\infty}(\Omega)$.

Given the freedom in choosing $\omega$ and $\phi$ one can show that (2)–(11)
and (14)–(15) are _equivalent_ if $\mathbf{u}$ is sufficiently smooth and
divergence free. We formulate it as a proposition.

###### Proposition 1

Assume $\mathbf{u}$ and $p$ are smooth and $\operatorname{div}\mathbf{u}=0$,
then (2) (or (11)) holds for any subdomain $\omega\subset\Omega$ iff (14) (or
(15)) holds for any $\phi\in W^{1,\infty}(\Omega)$ with
$\mbox{supp}(\phi)\subset\Omega$.

Proof. We know that (2) and $\operatorname{div}\mathbf{u}=0$ imply (1) by
standard arguments, given that $\omega$ can be taken as an arbitrary subdomain
of $\Omega$ and for any $t$. Similarly, the fact that (14) holds for any
$\phi\in\accentset{\mbox{\large.}}{C}(\Omega)$ leads to

$\int_{\Omega}\Big{(}\frac{\partial\mathbf{u}}{\partial
t}+(\mathbf{u}\cdot\nabla)\mathbf{u}-2\nu\operatorname{div}\mathbf{D}(\mathbf{u})+\nabla
p\Big{)}\phi\,dx=0,\quad\forall\,\phi\in\accentset{\mbox{\large.}}{C}(\Omega),$

which implies (1) due to the density of smooth compactly supported functions
in $L^{2}(\Omega)$. In turn, both (2) and (14) are quick consequences of (1).
Thus (2) implies (14) and vice versa.

The same arguments can be applied to show the equivalence of (11) and (15).
The only difference is that the equivalence of (11) is established not to the
momentum equation in (1), but to the vector product of this equation with
$\mathbf{x}$, and also for (15).

In addition to the above equivalence result, it is easy to see that each
individual term in (2)–(11) can be approximated arbitrarily well by the
corresponding term in (14)–(15). Indeed, fix any $\omega\subset\Omega$ with
smooth $\partial\omega$ and for sufficiently small $\varepsilon>0$ define

$\phi_{\varepsilon}=\begin{cases}\varepsilon^{-1}\mbox{dist}(\mathbf{x},\partial\omega),&\mathbf{x}\in\mathcal{O}_{\varepsilon}(\partial\omega)\cap\omega,\\\
1&\mathbf{x}\in\omega\setminus\mathcal{O}_{\varepsilon}(\partial\omega),\\\
0&\Omega\setminus\omega.\end{cases}$ (16)

We have $\phi_{\varepsilon}\in W^{1,\infty}(\Omega)$ and one easily checks,
letting
$\mathbf{\tilde{n}}=-\nabla\phi_{\varepsilon}/|\nabla\phi_{\varepsilon}|$,
that

$\frac{d}{dt}\int_{\omega}\phi_{\varepsilon}\mathbf{u}\,dx\to\frac{d}{dt}\int_{\omega}\mathbf{u}\,ds,\quad\int_{\omega}\mathbf{D}(\mathbf{u})\mathbf{\tilde{n}}|\nabla\phi_{\varepsilon}|\,dx\to\int_{\partial\omega}\mathbf{D}(\mathbf{u})\mathbf{n}\,ds,~{}\text{for}~{}\varepsilon\to
0,$

and smooth $\mathbf{u}$. Similarly, the limit values of other terms in
(14)–(15) will be their counterparts in (2)–(11). Therefore, (14)–(15) can be
also interpreted as the _diffuse-volume_ version of conservation laws. Eqs.
(14)–(15) imply (2)–(11) term by term without the
$\operatorname{div}\mathbf{u}=0$ assumption or equations (1) being invoked.

Local conservation laws in the Lagrangian form are written for the time-
dependent material volume $\Omega_{t}$. Denote by $\Omega_{0}$ the fluid
volume at a given initial moment $t=t_{0}$ and assume
$\Omega_{t}\subset\Omega$ for $t\in[t_{0},t_{1}]$ for some $t_{1}>t_{0}$. The
evolution of $\Omega_{t}$ is defined by the Lagrangian mapping
$\Phi_{t}:\Omega_{0}\to\Omega_{t}$, i.e. $\mathbf{y}=\Phi_{t}(\mathbf{x})$
solves the Cauchy problem

$\mathbf{y}_{t}=\mathbf{u}(t,\mathbf{y}),\quad
t\in(t_{0},t_{1}],\quad\mathbf{y}(t_{0})=\mathbf{x}.$ (17)

To properly reflect this domain evolution in a weak form of (12)–(13), we want
$\phi$ to be time dependent and such that $\mbox{supp}(\phi)=\Omega_{t}$. To
this end, consider a smooth function $\phi^{0}$ such that
$\mbox{supp}(\phi^{0})=\Omega_{0}$. We define
$\phi=\phi^{0}\circ\Phi_{t}^{-1}.$ The constructed $\phi$ is smooth (since
$\mathbf{u}$ is smooth so is the solution to the Cauchy problem (17)),
$\mbox{supp}(\phi)=\Omega_{t}$, and it satisfies the transport equation

$\frac{\partial\phi}{\partial
t}+(\mathbf{u}\cdot\nabla)\phi=0\quad\text{in}~{}\Omega,~{}t\in(t_{0},t_{1}],\quad\phi(t_{0})=\phi^{0}.$
(18)

Applying the Reynolds’ transport theorem and using (18) and (1) one computes
the following weak Lagrangian form of the local balances:

$\displaystyle\text{Moment.}\hskip
30.1388pt\frac{d}{dt}\int_{\Omega_{t}}\phi\mathbf{u}\,dx$
$\displaystyle~{}=~{}2\nu\int_{\Omega_{t}}\mathbf{D}(\mathbf{u})\mathbf{\tilde{n}}|\nabla\phi|\,dx-\int_{\Omega_{t}}p\mathbf{\tilde{n}}|\nabla\phi|\,dx,$
(19) $\displaystyle\text{Angl.
Moment.}\quad\frac{d}{dt}\int_{\Omega_{t}}\phi\mathbf{u}\times\mathbf{x}\,dx$
$\displaystyle~{}=~{}2\nu\int_{\Omega_{t}}(\mathbf{D}(\mathbf{u})\mathbf{\tilde{n}})\times\mathbf{x}|\nabla\phi|\,dx-\int_{\Omega_{t}}p(\mathbf{\tilde{n}}\times\mathbf{x})|\nabla\phi|\,dx.$
(20)

By the same arguments as we use to prove Proposition 1 we prove the following
proposition.

###### Proposition 2

Assume $\mathbf{u}$ and $p$ are smooth and $\operatorname{div}\mathbf{u}=0$.
Then (12) (or (13)) holds for any material volume $\Omega_{t}$ such that
$\Omega_{t}\subset\Omega$ for $t\in[t_{0},t_{1}]$ iff (19) (or (20)) holds for
any $\phi$ satisfying (18) with $\phi^{0}\in W^{1,\infty}(\Omega_{t_{0}})$,
such that $\mbox{supp}(\phi^{0})=\Omega_{t_{0}}$.

Similar to the Eulerian case, it is easy to see that each individual term in
(12)–(13) can be approximated arbitrarily well by the corresponding term in
(19)–(20). This time $\phi_{\varepsilon}$ is constructed as
$\phi_{\varepsilon}=\phi^{0}_{\varepsilon}\circ\Phi_{t}^{-1}$ with
$\phi^{0}_{\varepsilon}$ defined by the formula in (16) with $\omega$ replaced
by $\Omega_{t_{0}}$. It holds $\phi_{\varepsilon}\in
W^{1,\infty}(\Omega\times[t_{0},t_{1}])$ and one verifies, letting
$\mathbf{\tilde{n}}=-\nabla\phi_{\varepsilon}/|\nabla\phi_{\varepsilon}|$,
that

$\frac{d}{dt}\int_{\Omega_{t}}\phi_{\varepsilon}\mathbf{u}\,dx\to\frac{d}{dt}\int_{\Omega_{t}}\mathbf{u}\,ds,\quad\int_{\Omega_{t}}\mathbf{D}(\mathbf{u})\mathbf{\tilde{n}}|\nabla\phi_{\varepsilon}|\,dx\to\int_{\partial\Omega_{t}}\mathbf{D}(\mathbf{u})\mathbf{n}\,ds,~{}\text{for}~{}\varepsilon\to
0,$

and smooth $\mathbf{u}$. The limit values of other terms in (19)–(20) will be
their counterparts in (12)–(13). Therefore, (19)–(20) can be also interpreted
as the diffuse-volume version of local conservation laws in the Lagrangian
form.

In summary, equations (14)–(15) are equivalent formulations of the fundamental
(local) conservation laws in the Eulerian formulation, while (19)–(20) are
equivalent formulations of the fundamental (local) conservation laws in the
Lagrangian formulation. We will study the ability of a discretization method
to match (14)–(15) and (19)–(20) instead of (2)–(11) and (12)–(13).

##  4 EMAC Galerkin formulation is locally conservative

As an example of a continuous Galerkin method, we consider a conforming finite
element method: Denote by $\mathbf{V}_{h}\subset H^{1}_{0}(\Omega)^{d}$ and
$Q_{h}\subset L^{2}_{0}(\Omega)$ velocity and pressure finite element spaces
with respect to a tessellation $\mathcal{T}_{h}$ of $\Omega$ into elements
(simplexes or more general polygons or polyhedra). We also need the following
auxiliary spaces of continuous finite elements of degree $m+1$ and $m$, with
$m\geq 1$:

$\begin{split}V_{h}&=\\{v\in
H^{1}_{0}(\Omega):\,v\in\mathbb{P}_{m+1}(T)~{}\forall\,T\in\mathcal{T}_{h}\\},\\\
\widetilde{V}_{h}&=\\{v\in
H^{1}_{0}(\Omega):\,v\in\mathbb{P}_{m}(T)~{}~{}~{}~{}\forall\,T\in\mathcal{T}_{h}\\}.\end{split}$
(21)

We only assume that the velocity space contains all piecewise polynomial
continuous functions of degree $m+1$, i.e.

$(V_{h})^{d}\subset\mathbf{V}_{h}.$ (22)

We do not have any further assumptions on finite element spaces, and in
particular, both LBB stable and stabilized finite elements are admitted.

###### Remark 1

Let $\mathcal{T}_{h}$ be a triangulation of $\Omega$ and $m\geq 1$ be a
polynomial degree. The following examples of LBB stable FE pairs satisfy the
assumption: generalized Taylor–Hood $P_{m+1}-P_{m}$, $P_{m+1}-P_{m-1}^{\rm
disc}$ (for $d=2$), $P_{m+1}-P_{m-2}^{\rm disc}$ (for $d=3$, $m>1$),
$P_{m+1}^{\rm bubble}-P_{m-1}^{\rm disc}$ (for $d=3$ with face bubbles),
generalized conforming Crouzeix–Raviart $P_{m+1}^{\rm bubble}-P_{m}^{\rm
disc}$, Scott-Vogelius $P_{m+1}-P_{m}^{\rm disc}$ (SV element is LBB stable
subject to further assumptions on $\mathcal{T}_{h}$ [9]), as well as LBB
unstable equal order $P_{m+1}-P_{m+1}$ elements.

We use $(f,g):=\int_{\Omega}f\cdot g\,dx$ notation for both scalar and vector
functions $f,g$. The EMAC Galerkin formulation of (1) with $\mathbf{u}={\bf
0}$ on $\partial\Omega$ reads: Find $\mathbf{u}_{h}:(0,T)\to\mathbf{V}_{h}$
and $\widehat{p}_{h}:(0,T)\to Q_{h}\cap L^{2}_{0}(\Omega)$

$\left\\{\begin{aligned} \Big{(}\frac{\partial\mathbf{u}_{h}}{\partial
t},\mathbf{v}_{h}\Big{)}+2(\mathbf{D}(\mathbf{u}_{h})\mathbf{u}_{h},\mathbf{v}_{h})+((\operatorname{div}\mathbf{u}_{h})\mathbf{u}_{h},\mathbf{v}_{h})+2\nu(\mathbf{D}(\mathbf{u}_{h}),\mathbf{D}(\mathbf{v}_{h}))+(\widehat{p}_{h},\operatorname{div}\mathbf{v}_{h})&=0\quad\forall\,\mathbf{v}_{h}\in\mathbf{V}_{h},\\\
(\operatorname{div}\mathbf{u}_{h},q_{h})&=0\quad\forall\,q_{h}\in
Q_{h},\end{aligned}\right.$ (23)

where $\widehat{p}_{h}$ approximates the EMAC pressure
$\widehat{p}=p-\tfrac{1}{2}|\mathbf{u}|^{2}$. The EMAC formulation is
equivalent to other commonly used discrete formulations if
$\operatorname{div}\mathbf{u}_{h}=0$ pointwise. However,
$(\operatorname{div}\mathbf{u}_{h},q_{h})=0$ does not imply
$\operatorname{div}\mathbf{u}_{h}=0$ except in special settings. As a
consequence, the discrete solution depends on the form of nonlinear terms used
(i.e. EMAC, SKEW, CONV, ROT, etc.).

Next, we demonstrate that the solution of (23) satisfies discrete counterparts
of local conservation laws in both Eulerian and Lagrangian forms.

###  4.1 Local conservation in Eulerian form

Unlike for the continuous problem, for the discrete case the counterparts of
conservation laws in Eulerian and Lagrangian forms do not follow one from
another and we have to consider them separately. We start with the Eulerian
form.

Conservation of local linear momentum. Consider arbitrary $\phi_{h}\in V_{h}$,
$\phi_{h}|_{\partial\Omega}=0$. Then
$\phi_{h}\mathbf{e}_{i}\in\mathbf{V}_{h}$, for $i=1,\dots,d$, is a legitimate
test function in (23). Letting $\mathbf{v}_{h}=\phi_{h}\mathbf{e}_{i}$ in (23)
we compute for the nonlinear term

$2(\mathbf{D}(\mathbf{u}_{h})\mathbf{u}_{h},\phi_{h}\mathbf{e}_{i})=(\mathbf{u}_{h}\cdot\nabla\mathbf{u}_{h},\phi_{h}\mathbf{e}_{i})+((\phi_{h}\mathbf{e}_{i})\cdot\nabla\mathbf{u}_{h},\mathbf{u}_{h})\\\
=-(\mathbf{u}_{h}\cdot\nabla(\phi_{h}\mathbf{e}_{i}),\mathbf{u}_{h})-((\operatorname{div}\mathbf{u}_{h})\mathbf{u}_{h},\phi_{h}\mathbf{e}_{i})-\tfrac{1}{2}(\operatorname{div}(\phi_{h}\mathbf{e}_{i})\mathbf{u}_{h},\mathbf{u}_{h})\\\
=-(\mathbf{u}_{h}\cdot\nabla\phi_{h},\mathbf{u}_{h}\cdot\mathbf{e}_{i})-((\operatorname{div}\mathbf{u}_{h})\mathbf{u}_{h},\phi_{h}\mathbf{e}_{i})-\tfrac{1}{2}(\mathbf{e}_{i}\cdot\nabla\phi_{h},|\mathbf{u}_{h}|^{2}).$
(24)

Substituting this in the first equation from (23) with
$\mathbf{v}_{h}=\phi_{h}\mathbf{e}_{i}$ we obtain

$\Big{(}\frac{\partial\mathbf{u}_{h}}{\partial
t},\phi_{h}\mathbf{e}_{i}\Big{)}-(\mathbf{u}_{h}\cdot\nabla\phi_{h},\mathbf{u}_{h}\cdot\mathbf{e}_{i})-\tfrac{1}{2}(\mathbf{e}_{i}\cdot\nabla\phi_{h},|\mathbf{u}_{h}|^{2})+2\nu(\mathbf{D}(\mathbf{u}_{h}),\mathbf{D}(\phi_{h}\mathbf{e}_{i}))-(\widehat{p}_{h},\operatorname{div}(\phi_{h}\mathbf{e}_{i}))=0,$
(25)

and after simple re-arrangements,

$\frac{d}{dt}\Big{(}\mathbf{u}_{h}\cdot\mathbf{e}_{i},\phi_{h}\Big{)}-(\mathbf{u}_{h}\cdot\nabla\phi_{h},\mathbf{u}_{h}\cdot\mathbf{e}_{i})+2\nu(\mathbf{D}(\mathbf{u}_{h})\nabla\phi_{h},\mathbf{e}_{i})-(\widehat{p}_{h}+\tfrac{1}{2}|\mathbf{u}_{h}|^{2},\mathbf{e}_{i}\cdot\nabla\phi_{h})=0.$

Let $\mathbf{n}_{h}:=-\nabla\phi_{h}/|\nabla\phi_{h}|$ for
$|\nabla\phi_{h}|\neq 0$ (and arbitrary unit vector otherwise) and define

$\omega_{h}=\mbox{supp}(\phi_{h})\quad\text{and}~{}~{}p_{h}=\widehat{p}_{h}+\tfrac{1}{2}|\mathbf{u}_{h}|^{2},$

then from equation (25) for $i=1,\dots,d$ we get,

$\frac{d}{dt}\int_{\omega_{h}}\phi_{h}\mathbf{u}_{h}\,dx=2\nu\int_{\omega_{h}}\mathbf{D}(\mathbf{u}_{h})\mathbf{n}_{h}|\nabla\phi_{h}|\,dx-\int_{\omega_{h}}p_{h}\mathbf{n}_{h}|\nabla\phi_{h}|\,dx-\int_{\omega_{h}}\mathbf{u}_{h}(\mathbf{u}_{h}\cdot\mathbf{n}_{h})|\nabla\phi_{h}|\,dx,$
(26)

for any $\phi_{h}\in V_{h}$. This is the discrete analogue of the local
momentum conservation in (14).

Conservation of local angular momentum. Consider arbitrary
$\phi_{h}\in\widetilde{V}_{h}$. Then
$\mathbf{x}\times\phi_{h}\mathbf{e}_{i}\in\mathbf{V}_{h}$ for $i=1,\dots,d$ is
a legitimate test function. Letting
$\mathbf{v}_{h}=\mathbf{x}\times\phi_{h}\mathbf{e}_{i}$ in (23) we compute for
the nonlinear term

$2(\mathbf{D}(\mathbf{u}_{h})\mathbf{u}_{h},\mathbf{x}\times\phi_{h}\mathbf{e}_{i})=(\mathbf{u}_{h}\cdot\nabla\mathbf{u}_{h},\mathbf{x}\times\phi_{h}\mathbf{e}_{i})+((\mathbf{x}\times\phi_{h}\mathbf{e}_{i})\cdot\nabla\mathbf{u}_{h},\mathbf{u}_{h})\\\
=-(\mathbf{u}_{h}\cdot\nabla(\mathbf{x}\times\phi_{h}\mathbf{e}_{i}),\mathbf{u}_{h})-((\operatorname{div}\mathbf{u}_{h})\mathbf{u}_{h},\mathbf{x}\times\phi_{h}\mathbf{e}_{i})-\tfrac{1}{2}(\operatorname{div}(\mathbf{x}\times\phi_{h}\mathbf{e}_{i})\mathbf{u}_{h},\mathbf{u}_{h})\\\
=-(\mathbf{u}_{h}\cdot\nabla\phi_{h},(\mathbf{u}_{h}\times\mathbf{x})\cdot\mathbf{e}_{i})-((\operatorname{div}\mathbf{u}_{h})\mathbf{u}_{h},\mathbf{x}\times\phi_{h}\mathbf{e}_{i})-\tfrac{1}{2}(\mathbf{x}\times\nabla\phi_{h},\mathbf{e}_{i}|\mathbf{u}_{h}|^{2}).$
(27)

Substituting this in the first equation from (23) with
$\mathbf{v}_{h}=\mathbf{x}\times\phi_{h}\mathbf{e}_{i}$ we obtain

$\Big{(}\frac{\partial\mathbf{u}_{h}}{\partial
t},\mathbf{x}\times\phi_{h}\mathbf{e}_{i}\Big{)}-(\mathbf{u}_{h}\cdot\nabla\phi_{h},(\mathbf{u}_{h}\times\mathbf{x})\cdot\mathbf{e}_{i})-\tfrac{1}{2}(\mathbf{x}\times\nabla\phi_{h},\mathbf{e}_{i}|\mathbf{u}_{h}|^{2})\\\
+2\nu(\mathbf{D}(\mathbf{u}_{h}),\mathbf{D}(\mathbf{x}\times\phi_{h}\mathbf{e}_{i}))-(\widehat{p}_{h},\operatorname{div}(\mathbf{x}\times\phi_{h}\mathbf{e}_{i}))=0.$
(28)

Simple re-arrangements give

$\frac{d}{dt}\Big{(}\mathbf{u}_{h}\times\mathbf{x},\phi_{h}\mathbf{e}_{i}\Big{)}-(\mathbf{u}_{h}\cdot\nabla\phi_{h},(\mathbf{u}_{h}\times\mathbf{x})\cdot\mathbf{e}_{i})+2\nu(\mathbf{D}(\mathbf{u}_{h})\nabla\phi_{h},\mathbf{x}\times\mathbf{e}_{i})-((\widehat{p}_{h}+\tfrac{1}{2}|\mathbf{u}_{h}|^{2})\mathbf{e}_{i},\nabla\phi_{h}\times\mathbf{x})=0.$

From the above equality for $i=1,\dots,d$ we get

$\frac{d}{dt}\int_{\omega_{h}}\phi_{h}\mathbf{u}_{h}\times\mathbf{x}\,dx=2\nu\int_{\omega_{h}}(\mathbf{D}(\mathbf{u}_{h})\mathbf{n}_{h})\times\mathbf{x}\,|\nabla\phi_{h}|\,dx-\int_{\omega_{h}}p_{h}(\mathbf{n}_{h}\times\mathbf{x})|\nabla\phi_{h}|\,dx-\int_{\omega_{h}}(\mathbf{u}_{h}\times\mathbf{x})(\mathbf{u}_{h}\cdot\mathbf{n}_{h})|\nabla\phi_{h}|\,dx$
(29)

for any $\phi_{h}\in\widetilde{V}_{h}$. This is the discrete analogue of the
local angular momentum conservation from (15).

### Local conservation in Lagrangian form

After discretization, there is no obvious equivalence between the Eulerian and
Lagrangian forms of the local balances. Nevertheless, one can show that EMAC
form also obeys a discrete counterpart of the linear momentum local
conservation in the Lagrangian form. However, we need additional assumption on
$V_{h}$ space. Namely, we assume that the velocity space consists of piecewise
polynomial continuous functions of degree $m+1$:

$(V_{h})^{d}=\mathbf{V}_{h}.$ (30)

Conservation of local linear momentum. Consider $\phi_{h}^{0}\in V_{h}$ and
$\phi_{h}:\,[t_{0},\hat{t}_{1}]\to V_{h}$ solving

$\Big{(}\frac{\partial\phi_{h}}{\partial
t},v_{h}\Big{)}+(\mathbf{u}_{h}\cdot\nabla\phi_{h},v_{h})=0\quad\forall~{}v_{h}\in
V_{h},$ (31)

which is the projection of the transport equation (18) on the finite
dimensional space $V_{h}$ with $\mathbf{u}$ replaced by $\mathbf{u}_{h}$.

Letting $\mathbf{v}_{h}=\phi_{h}\mathbf{e}_{i}$ in (23), we repeat the same
calculations as for the Eulerian case and arrive at (25). Since $\phi_{h}$ is
time dependent, after re-arrangements (25) gives

$\frac{d}{dt}\Big{(}\mathbf{u}_{h}\cdot\mathbf{e}_{i},\phi_{h}\Big{)}-\Big{(}\mathbf{u}_{h}\cdot\mathbf{e}_{i},\frac{\partial\phi_{h}}{\partial
t}\Big{)}-(\mathbf{u}_{h}\cdot\nabla\phi_{h},\mathbf{u}_{h}\cdot\mathbf{e}_{i})+2\nu(\mathbf{D}(\mathbf{u}_{h})\nabla\phi_{h},\mathbf{e}_{i})-(\widehat{p}_{h}+\tfrac{1}{2}|\mathbf{u}_{h}|^{2},\mathbf{e}_{i}\cdot\nabla\phi_{h})=0.$

Thanks to the assumption (30) and equation (31), the second and third terms
add to zero.

Let $\mathbf{n}_{h}:=-\nabla\phi_{h}/|\nabla\phi_{h}|$ for
$|\nabla\phi_{h}|\neq 0$ (and arbitrary unit vector otherwise) and define

$\Omega_{h}(t)=\mbox{supp}(\phi_{h})\quad\text{and}~{}~{}p_{h}=\widehat{p}_{h}+\tfrac{1}{2}|\mathbf{u}_{h}|^{2},$

then from equation (25) for $i=1,\dots,d$ we get,

$\frac{d}{dt}\int_{\Omega_{h}(t)}\phi_{h}\mathbf{u}_{h}\,dx=2\nu\int_{\Omega_{h}(t)}\mathbf{D}(\mathbf{u}_{h})\mathbf{n}_{h}|\nabla\phi_{h}|\,dx-\int_{\Omega_{h}(t)}p_{h}\mathbf{n}_{h}|\nabla\phi_{h}|\,dx,$
(32)

for any $\phi_{h}\in V_{h}$. This is the discrete analogue of the local
momentum conservation in (19).

Conservation of local angular momentum. Consider
$\phi_{h}^{0}\in\widetilde{V}_{h}$ and
$\phi_{h}:\,[t_{0},\hat{t}_{1}]\to\widetilde{V}_{h}$ solving

$\Big{(}\frac{\partial\phi_{h}}{\partial
t},v_{h}\Big{)}+(\mathbf{u}_{h}\cdot\nabla\phi_{h},v_{h})=0\quad\forall~{}v_{h}\in\widetilde{V}_{h}.$

Letting $\mathbf{v}_{h}=\mathbf{x}\times\phi_{h}\mathbf{e}_{i}$ in (23) we
repeat the same calculations as for the Eulerian case and arrive at (28).
Since $\phi_{h}$ is time dependent, after re-arrangements (28) gives

$\frac{d}{dt}\Big{(}\mathbf{u}_{h}\times\mathbf{x},\phi_{h}\mathbf{e}_{i}\Big{)}-\Big{(}\frac{\partial\phi_{h}}{\partial
t},(\mathbf{u}_{h}\times\mathbf{x})\cdot\mathbf{e}_{i}\Big{)}-(\mathbf{u}_{h}\cdot\nabla\phi_{h},(\mathbf{u}_{h}\times\mathbf{x})\cdot\mathbf{e}_{i})\\\
+2\nu(\mathbf{D}(\mathbf{u}_{h})\nabla\phi_{h},\mathbf{x}\times\mathbf{e}_{i})-((\widehat{p}_{h}+\tfrac{1}{2}|\mathbf{u}_{h}|^{2})\mathbf{e}_{i},\nabla\phi_{h}\times\mathbf{x})=0.$

Denote by $I_{m}(\mathbf{u}_{h}\times\mathbf{x})$ a piecewise polynomial of
degree $m$ interpolating $\mathbf{u}_{h}\times\mathbf{x}$, i.e.
$I_{m}(\mathbf{u}_{h}\times\mathbf{x})\in\widetilde{V}_{h}^{3}$. Then
$I_{m}(\mathbf{u}_{h}\times\mathbf{x})\cdot\mathbf{e}_{i}\in\widetilde{V}_{h}$
holds. Therefore,

$\Big{(}\frac{\partial\phi_{h}}{\partial
t},(\mathbf{u}_{h}\times\mathbf{x})\cdot\mathbf{e}_{i}\Big{)}+(\mathbf{u}_{h}\cdot\nabla\phi_{h},(\mathbf{u}_{h}\times\mathbf{x})\cdot\mathbf{e}_{i})=\left(\frac{d\phi_{h}}{dt},(\mathbf{u}_{h}\times\mathbf{x}-I_{m}(\mathbf{u}_{h}\times\mathbf{x}))\cdot\mathbf{e}_{i}\right)=:R_{i},$

where $\frac{d\phi_{h}}{dt}=\frac{\partial\phi_{h}}{\partial
t}+\mathbf{u}_{h}\cdot\nabla\phi_{h}$. Let
$\mathbf{n}_{h}:=-\nabla\phi_{h}/|\nabla\phi_{h}|$ for $|\nabla\phi_{h}|\neq
0$ (and arbitrary unit vector otherwise) and define

$\Omega_{h}(t)=\mbox{supp}(\phi_{h})\quad\text{and}~{}~{}p_{h}=\widehat{p}_{h}+\tfrac{1}{2}|\mathbf{u}_{h}|^{2},$

then from equation (25) for $i=1,\dots,d$ we get,

$\frac{d}{dt}\int_{\Omega_{h}(t)}\phi_{h}\mathbf{u}_{h}\times\mathbf{x}\,dx=2\nu\int_{\Omega_{h}(t)}(\mathbf{D}(\mathbf{u}_{h})\mathbf{n}_{h})\times\mathbf{x}\,|\nabla\phi_{h}|\,dx-\int_{\Omega_{h}(t)}p_{h}(\mathbf{n}_{h}\times\mathbf{x})|\nabla\phi_{h}|\,dx+R$
(33)

for any $\phi_{h}\in\widetilde{V}_{h}$. This is the discrete analogue of the
local momentum conservation in (20) _up to the residual term_
$R=R_{1}+\dots+R_{d}$. If we assume that $\mathbf{u}_{h}$ approximates a
(smooth) solution to the NSE with order $O(h^{r})$, $r\geq m+1$, in some norm
$\|\cdot\|_{\ast}$ then $\|R\|_{\ast}=O(h^{m+1})$ once $\frac{d\phi_{h}}{dt}$
is bounded in the dual norm to $\|\cdot\|_{\ast}$. According to (30) the
optimal approximation order for $\mathbf{u}_{h}$ would be $O(h^{m+2})$ in
$L^{2}(L^{2})$ norm.

##  5 Numerical Tests

We now give numerical examples to illustrate the theory above. For these
tests, the full Navier–Stokes discretization uses BDF temporal discretizations
(and Crank–Nicolson for the initial time steps) and
$(\mathbf{V}_{h},Q_{h})=(P_{2}(\tau_{h})^{d},P_{1}(\tau_{h}))$ Taylor–Hood
elements on a mesh $\tau_{h}$. The nonlinear problem at each time step is
resolved with Newton’s method, and typically it takes just 2 or 3 iterations
to resolve. With temporal discretizations, the precise definitions of the
discrete local balances will change accordingly, and we derive these now
before proceeding to the tests.

Denote by $\omega_{h}$ the approximation of a subdomain $\omega$ whose
boundary consists of element edges from the mesh. Define the functions
$\phi_{h}\in V_{h}$ and $\psi_{h}\in\tilde{V}_{h}$ nodally by

$\displaystyle\phi_{h}(x_{j})=\left\\{\begin{array}[]{c}1\ \mbox{ if $x_{j}$
is a node on $P_{2}(\tau_{h})$ in the interior of $\omega_{h}$}\\\ 0\ \mbox{
otherwise, }\end{array}\right.$ (36)
$\displaystyle\psi_{h}(x_{j})=\left\\{\begin{array}[]{c}1\ \mbox{ if $x_{j}$
is a node on $P_{1}(\tau_{h})$ in the interior of $\omega_{h}$}\\\ 0\ \mbox{
otherwise. }\end{array}\right.$ (39)

In our implementations we apply BDF formulas for the temporal discretization
for the momentum and transport equations. For a sequence
$\\{f^{n}\\}_{n=0,1,\dots}$ of scalar or vector quantities, we use the
shortcut notations

$\displaystyle\left(\frac{df}{dt}\right)^{n}_{\rm bdf3}$
$\displaystyle=\frac{\frac{11}{6}f^{n}-3f^{n-1}+\frac{3}{2}f^{n-2}-\frac{1}{3}f^{n-3}}{\Delta
t},$ $\displaystyle\left(\frac{df}{dt}\right)^{n}_{\rm bdf2}$
$\displaystyle=\frac{3f^{n}-4f^{n-1}+f^{n-2}}{2\Delta t},$
$\displaystyle\left(\frac{df}{dt}\right)^{n}_{\rm bdf1}$
$\displaystyle=\frac{f^{n}-f^{n-1}}{\Delta t}.$

####  5.0.1 Discrete local conservation in Eulerian form

We consider first the discrete Eulerian form of local conservation of momentum
and angular momentum. Choosing $\phi_{h}$ by (36) and repeating the arguments
above that derived (26) but using the BDFk (k=1,2,3) temporal discretization,
we get the following (fully) discrete local momentum balance

$\displaystyle\left(\frac{d\Big{(}\int_{\omega_{h}}\phi_{h}\mathbf{u}_{h}\,dx\Big{)}}{dt}\right)^{n}_{\rm
bdfk}=2\nu\int_{\omega_{h}}\mathbf{D}(\mathbf{u}_{h}^{n})\mathbf{n}_{h}|\nabla\phi_{h}|\,dx-\int_{\omega_{h}}p_{h}^{n}\mathbf{n}_{h}|\nabla\phi_{h}|\,dx-\int_{\omega_{h}}\mathbf{u}_{h}^{n}(\mathbf{u}_{h}^{n}\cdot\mathbf{n}_{h})|\nabla\phi_{h}|\,dx,$

with $\mathbf{n}_{h}=-\nabla\phi_{h}/|\nabla\phi_{h}|$. Similarly, for
discrete local angular momentum conservation we obtain

$\displaystyle\left(\frac{d\Big{(}\int_{\omega_{h}}\psi_{h}\mathbf{u}_{h}\times\mathbf{x}\,dx\Big{)}}{dt}\right)^{n}_{\rm
bdfk}=$ $\displaystyle
2\nu\int_{\omega_{h}}(\mathbf{D}(\mathbf{u}_{h}^{n})\mathbf{n}_{h})\times\mathbf{x}\,|\nabla\psi_{h}|\,dx$
$\displaystyle-\int_{\omega_{h}}p_{h}^{n}(\mathbf{n}_{h}\times\mathbf{x})|\nabla\psi_{h}|\,dx-\int_{\omega_{h}}(\mathbf{u}_{h}^{n}\times\mathbf{x})(\mathbf{u}_{h}^{n}\cdot\mathbf{n}_{h})|\nabla\psi_{h}|\,dx,$

where $\psi_{h}$ is defined by (39), and
$\mathbf{n}_{h}=-\nabla\psi_{h}/|\nabla\psi_{h}|$.

In our tests, we will show plots of discrete local Eulerian momentum error

$e_{E}^{mom}=\left(\frac{d\Big{(}\int_{\omega_{h}}\phi_{h}\mathbf{u}_{h}\,dx\Big{)}}{dt}\right)^{n}_{\rm
bdfk}-2\nu\int_{\omega_{h}}\mathbf{D}(\mathbf{u}_{h}^{n})\mathbf{n}_{h}|\nabla\phi_{h}|\,dx\\\
+\int_{\omega_{h}}p_{h}^{n}\mathbf{n}_{h}|\nabla\phi_{h}|\,dx+\int_{\omega_{h}}\mathbf{u}_{h}^{n}(\mathbf{u}_{h}^{n}\cdot\mathbf{n}_{h})|\nabla\phi_{h}|\,dx,$

and discrete local Eulerian angular momentum error

$e_{E}^{am}=\left(\frac{d\Big{(}\int_{\omega_{h}}\psi_{h}\mathbf{u}_{h}\times\mathbf{x}\,dx\Big{)}}{dt}\right)^{n}_{\rm
bdfk}-2\nu\int_{\omega_{h}}(\mathbf{D}(\mathbf{u}_{h}^{n})\mathbf{n}_{h})\times\mathbf{x}\,|\nabla\psi_{h}|\,dx\\\
+\int_{\omega_{h}}p_{h}^{n}(\mathbf{n}_{h}\times\mathbf{x})|\nabla\psi_{h}|\,dx+\int_{\omega_{h}}(\mathbf{u}_{h}^{n}\times\mathbf{x})(\mathbf{u}_{h}^{n}\cdot\mathbf{n}_{h})|\nabla\psi_{h}|\,dx.$

####  5.0.2 Discrete local conservation in Lagrangian form

Discrete local conservation in Lagrangian form is somewhat more complicated
compared to the Eulerian case due to the $\phi_{h}$ function becoming time
dependent in the momentum and angular momentum balances, as well as the
transport equations involved in these balances being hyperbolic. As our tests
are for illustrative purposes of certain theoretical properties, we
approximate the transport equations in the following way for the purpose of
ease in computations, even though other approaches to solving the transport
equation may be better in practice.

Consider $\phi_{h}^{0}\in V_{h}$ to be defined by (36), and then define
$\phi_{h}^{n}\in V_{h}$ (n=1,2,3,…) via

$\Big{(}\left(\frac{d\phi_{h}}{dt}\right)^{n}_{\rm
bdfj},v_{h}\Big{)}+(\mathbf{u}_{h}^{n}\cdot\nabla\phi_{h}^{n},v_{h})=0\quad\forall~{}v_{h}\in
V_{h},$ (40)

where $j=$1 or 2 in our numerical tests (and if $j=$2 then the first time step
is backward Euler).

Rederiving the discrete local Lagrangian momentum balance (32) but now using
BDFk ($k=$2 or 3) time stepping for Navier-Stokes together with (40), we
obtain the balance

$\int_{\Omega}\left(\phi_{h}\left(\frac{d\mathbf{u}_{h}}{dt}\right)^{n}_{\rm
bdfk}+\left(\frac{d\phi_{h}}{dt}\right)^{n}_{\rm
bdfj}\mathbf{u}_{h}^{n}\right)\
dx=2\nu\int_{\Omega}\mathbf{D}(\mathbf{u}_{h}^{n})\mathbf{n}^{n}_{h}|\nabla\phi^{n}_{h}|\,dx-\int_{\Omega}p_{h}^{n}\mathbf{n}^{n}_{h}|\nabla\phi^{n}_{h}|\,dx,$

and thus define the discrete local Lagrangian momentum error by

$e_{L}^{mom}=\int_{\Omega}\left(\phi_{h}^{n}\left(\frac{d\mathbf{u}_{h}}{dt}\right)^{n}_{\rm
bdfk}+\left(\frac{d\phi_{h}}{dt}\right)^{n}_{\rm
bdfj}\mathbf{u}_{h}^{n}\right)\
dx-\left(2\nu\int_{\Omega}\mathbf{D}(\mathbf{u}_{h}^{n})\mathbf{n}^{n}_{h}|\nabla\phi^{n}_{h}|\,dx-\int_{\Omega}p_{h}^{n}\mathbf{n}^{n}_{h}|\nabla\phi^{n}_{h}|\,dx\right).$

Note that if the transport equations are solved in a different way, then these
definitions of discrete Lagrangian momentum and angular momentum balances need
modified accordingly. For example if an explicit method is used, then the
local balance will be defined with some terms at time $t^{n}$ and others and
time $t^{n-1}$.

For angular momentum, we proceed similarly as for momentum to find a fully
discrete analogue to (33). Let $\psi_{h}^{0}$ be defined by (39) and find
$\psi_{h}^{n}\in\tilde{V}_{h}\cap H^{1}_{0}(\Omega)$ for $n=1,2,3,...$ by

$\Big{(}\left(\frac{d\psi_{h}}{dt}\right)^{n}_{\rm
bdfj},v_{h}\Big{)}+(\mathbf{u}_{h}^{n}\cdot\nabla\psi_{h}^{n},v_{h})=0\quad\forall~{}v_{h}\in\tilde{V}_{h}\cap
H^{1}_{0}(\Omega),$

and using backward Euler for the first time step if $j=$2\. Following similar
steps as the theory above, the discrete Lagrangian local angular momentum
error is then given by

$\displaystyle e_{L}^{am}=\bigg{(}\int_{\Omega}\bigg{(}\psi_{h}^{n}$
$\displaystyle\left(\frac{d\mathbf{u}_{h}}{dt}\right)^{n}_{\rm
bdfk}+\left(\frac{d\psi_{h}}{dt}\right)^{n}_{\rm
bdfj}\mathbf{u}_{h}^{n}\bigg{)}\times\mathbf{x}\,dx$
$\displaystyle-2\nu\int_{\Omega}(\mathbf{D}(\mathbf{u}_{h}^{n})\mathbf{n}_{h}^{n})\times\mathbf{x}\,|\nabla\psi_{h}^{n}|\,dx+\int_{\Omega}p_{h}^{n}(\mathbf{n}_{h}^{n}\times\mathbf{x})|\nabla\psi_{h}^{n}|\,dx\bigg{)}.$

###  5.1 Gresho problem

For our first test we use a slight variation of the classical Gresho problem
on $\Omega=(-0.5,0.5)^{2}$, which consists

Figure 1: Initial velocity for the Gresho problem is shown above, as speed
contours.

of a velocity and pressure

$\displaystyle\small\mathbf{u}=\begin{cases}\begin{bmatrix}-5y~{}5x\end{bmatrix}^{T}&\text{for
}r<2,\\\\[3.87498pt]
\begin{bmatrix}\frac{2y}{r}+5y~{}\frac{2x}{r}-5x\end{bmatrix}^{T}&\text{for
}.2\leq r\leq.4,\\\ \begin{bmatrix}0~{}0\end{bmatrix}^{T}&\text{for
}r>.4,\end{cases},$ $\displaystyle\small
p=\begin{cases}12.5r^{2}+C_{1}&\text{for }r<.2,\\\
12.5r^{2}-20r+4\log(r)+C_{2}&\text{for }.2\leq r\leq.4,\\\ 0&\text{for
}r>.4,\end{cases}$

where $r=\sqrt{x^{2}+y^{2}}$ and

$\displaystyle C_{2}=-12.5(.4)^{2}+20(.4)^{2}-4\log(.4),$ $\displaystyle
C_{1}=C_{2}-20(.2)+4\log(.2).$

This velocity is plotted in figure 1 and is an exact solution of the unforced
steady Euler equations, and hence an accurate solver should preserve the
initial condition in time. It is shown in [3, 20] that a NSE solver with EMAC
nonlinearity and using Crank-Nicolson time stepping together with Taylor-Hood
finite element spatial discretization will preserve pointwise global energy,
momentum and angular momentum for this problem while other common nonlinearity
formulations such as SKEW, ROT and CONV will not preserve these physical
balance laws and moreover will be less accurate in the sense of
$L^{2}(\Omega)$ error.

Figure 2: Shown above is the domain and $\omega$ (left), the mesh (center),
and the mesh zoomed in near $\omega$ for the Gresho problem.

We alter this problem very slightly by changing the viscosity to
$\nu=10^{-10}$ so as not to solve the Euler equations but instead the NSE. We
note this change of viscosity will (very slightly) change the true solution in
time, however this is of no consequence as our interest herein is not the
solution but the local conservation of momentum and angular momentum. We
choose $\omega$ to be the circle of radius $0.05$ centered at $(0.2,0.09)$, as
shown in figure 2 at left. Figure 2 also shows the mesh $\tau_{h}$ used for
the computations below as well as the mesh zoomed in near $\omega$. We define
$\omega_{h}$ to be the approximation of $\omega$ whose boundary consists of
triangle edges from the mesh.

Figure 3: Shown above is error in discrete local Eulerian (left) and
Lagrangian (right) momentum and angular momentum conservation versus time in
the (viscous) Gresho problem.

Computations are done using $T=0.5$, $\Delta t=0.005$ and ${\bf f}={\bf 0}$,
and we show errors in Eulerian conservation laws in figure 3 as absolute
values of errors versus time. We observe these quantities are conserved
pointwise and are stable in time, just as the theory above predicts.

Plots of discrete local Lagrangian momentum and angular momentum versus time
are also shown in figure 3, and we observe as expected that momentum is
preserved pointwise (up to machine roundoff error). We also observe that
Lagrangian angular momentum is not preserved pointwise but instead has values
as large as $O(10^{-4})$, which is consistent with the $O(h^{2})$ residual our
theory above predicts.

###  5.2 Kelvin-Helmholtz flow

Figure 4: Shown above are the absolute vorticity contours of the solution
velocity at $t=$0, 0.5, 1, 2, 3, 4, 4.5 and 5 (left to right, top to bottom).

For our second test we consider a test problem from [25] for simulating 2D
Kelvin-Helmholtz instability. The domain is the unit square, with periodic
boundary conditions at $x=0,1$, representing an infinite extension in the
horizontal direction. At $y=0,1$, we enforce for $t>0$ a no slip condition,
which differs from [25] as they use a no penetration and free slip condition.
However, as these boundaries are far from the physical behavior of interest,
there is little effect on the qualitative behavior of the solution. The
initial condition is set by

$\mathbf{u}_{0}(x,y)=\left(\begin{array}[]{c}u_{\infty}\tanh\left(\frac{2y-1}{\delta_{0}}\right)\\\
0\end{array}\right)+c_{n}\left(\begin{array}[]{c}\partial_{y}\psi(x,y)\\\
-\partial_{x}\psi(x,y)\end{array}\right),$

where $\delta_{0}=\frac{1}{28}$ is the initial vorticity thickness,
$u_{\infty}=1$ is a reference velocity, $c_{n}$ is a noise/scaling factor
taken to be $10^{-3}$, and

$\psi(x,y)=u_{\infty}\exp\left(-\frac{(y-0.5)^{2}}{\delta_{0}^{2}}\right)\left(\cos(8\pi
x)+\cos(20\pi x)\right).$

The Reynolds number is defined by
$Re=\frac{\delta_{0}u_{\infty}}{\nu}=\frac{1}{28\nu}$, and $\nu$ is defined by
selecting $Re$. We use $Re=$100 for our test.

We compute solutions for EMAC discretized with $(P_{2},P_{1})$ Taylor-Hood
elements on a $h=\frac{1}{128}$ uniform mesh, together with BDF2 time stepping
and a time step size of $\Delta t=0.01$. Solutions are computed up to $T=$5,
with plots of vorticity contours shown in figure 4 matching those in [25]
qualitatively well.

The subdomain $\omega$ is defined to be the square
$[\frac{1}{8},\frac{1}{4}]\times[\frac{1}{8},\frac{1}{4}]$. For this domain on
this mesh, we have that $\omega_{h}=\omega$. Plots of discrete Eulerian
momentum and angular momentum are shown in figure 5 at the top, and we observe
that these quantities are conserved exactly, just as the theory above
predicts. For discrete Lagrangian momentum and angular momentum, we solve the
transport equation using backward Euler (BDF1), and plots of momentum and
angular momentum are shown in figure 5 at bottom. We observe exact local
conservation of discrete Lagrangian momentum, and conservation of discrete
Lagrangian angular momentum consistent with the discretization error, as
predicted above.

Figure 5: Shown above is error in discrete local Eulerian (left) and
Lagrangian (right) momentum and angular momentum conservation versus time in
the Kelvin-Helmholtz problem.

###  5.3 2D flow past a cylinder

Figure 6: The domain for the channel flow past a cylinder numerical
experiment.  Figure 7: Shown above is error in discrete local Eulerian
momentum and angular momentum conservation versus time in the 2D channel flow
past a cylinder test.

For our last test we consider the 2D channel flow past a cylinder problem from
[24, 13]. The domain is the rectangle $[0,2.2]\times[0,0.41]$ as shown in
figure 6, with a cylinder centered at $(0.2,0.2)$ with radius $0.05.$ We take
no external forcing ${\bf f}={\bf 0}$, $\nu=0.001$ (which corresponds to
Reynolds number 100), and set inflow/outflow profiles to be

$\displaystyle u_{1}(0,y,t)$
$\displaystyle=u_{1}(2.2,y,t)=\frac{6}{0.41^{2}}y(0.41-y),$ $\displaystyle
u_{2}(0,y,t)$ $\displaystyle=u_{2}(2.2,y,t)=0.$

We define a subdomain $\omega$ to be a circle radius $0.05$ centered at
$(0.35,0.15)$, and $\omega_{h}$ to be its approximation by the mesh. A plot of
$\omega$ is shown in figure 6. We use a mesh that provides approximately 60K
velocity degrees of freedom (dof) and 7K pressure dof when discretized with
Taylor-Hood elements. We compute using BDF3 time stepping to $T=2$ using time
step size $\Delta t=0.001$, and start the flow from rest. Errors in discrete
local Eulerian momentum and angular momentum conservation are shown in figure
7 and we once again observe pointwise local conservation. We do not consider
Lagrangian discrete local conservation for this test, since there is an
outflow and conservation is therefore not expected except for very short times
as the transported quantity will exit the domain through the outflow.

##  6 Future directions

We have shown that continuous Galerkin discretizations of the Navier-Stokes
equations using EMAC nonlinearity form admit (appropriately defined) exact
local balances / conservation laws of momentum and angular momentum. These
discrete local balances are constructed as weak forms of the momentum and
angular momentum conservation laws, and are equivalent to the usual
conservation law definitions before discretization. In the discrete case,
however, these weak formulations are not equivalent to the usual conservation
law definitions and even their Eulerian and Lagrangian constructions are not
equivalent. That the discrete schemes admit any exact local balances at all is
quite rare, and we note that such an analysis is not possible for such common
Navier-Stokes nonlinearity formulations such as convective, skew-symmetric or
rotational. We remark that the ‘conservative’ formulation of nonlinear terms
(referred to as CONS, utilizing
$\operatorname{div}(\mathbf{u}_{h}\mathbf{u}_{h}^{T})$) also maintains the
same conservation properties for momenta as EMAC. However, CONS fails to
achieve a proper global energy balance, when
$\operatorname{div}\mathbf{u}_{h}\neq 0$, unlike EMAC. This deficiency leads
to unstable finite element schemes using CONS; see examples of CONS
underperformance in [3, 4].

Future directions for this work could include an extension of these ideas to
other conservation laws of Navier-Stokes such as energy, helicity, enstrophy
in 2D, vorticity, and others. It is currently unclear to the authors how to
construct appropriate local balances for these quantities for the continuous
Galerkin method.

##  7 Declaration of competing interest

The authors declare the following financial interests/personal relationships
which may be considered as potential competing interests: Both authors report
financial support was provided by National Science Foundation.

##  8 Acknowledgements

M.O. was partially supported by the NSF grant DMS-2309197. L.R. was partially
supported by the NSF grant DMS-2152623. We would like to express our gratitude
to Thomas J.R. Hughes for his insightful discussions that inspired our
research on the local conservation properties of EMAC.

## References

  * [1] R. Abramov and A. Majda, Discrete approximations with additional conserved quantities: deterministic and statistical behavior, Methods Appl. Anal., 10 (2003), pp. 151–190.
  * [2] A. Arakawa, Computational design for long-term numerical integration of the equations of fluid motion: Two dimensional incompressible flow, Part I, J. Comput. Phys., 1 (1966), pp. 119–143.
  * [3] S. Charnyi, T. Heister, M. Olshanskii, and L. Rebholz, On conservation laws of Navier-Stokes Galerkin discretizations, Journal of Computational Physics, 337 (2017), pp. 289–308.
  * [4] S. Charnyi, T. Heister, M. Olshanskii, and L. Rebholz, Efficient discretizations for the EMAC formulation of the incompressible Navier-Stokes equations, Applied Numerical Mathematics, 141 (2019), pp. 220–233.
  * [5] B. Cockburn, Discontinuous Galerkin methods, ZAMM-Journal of Applied Mathematics and Mechanics/Zeitschrift für Angewandte Mathematik und Mechanik: Applied Mathematics and Mechanics, 83 (2003), pp. 731–754.
  * [6] G. Coppola, F. Capuano, and L. de Luca, Discrete energy-conservation properties in the numerical simulation of the Navier–Stokes equations, Applied Mechanics Reviews, 71 (010803) (2019).
  * [7] J. A. Evans and T. J. R. Hughes, Isogeometric divergence-conforming B-splines for the unsteady Navier–Stokes equations, Journal of Computational Physics, 241 (2013), pp. 141–167.
  * [8] O. L. F. Sacco, B. Paun, T. Iles, P. Iaizzo, G. Houzeaux, M. Vzaquex, C. Butakoff, and J. Aguado-Sierra, Evaluating the roles of detailed endocardial structures on right ventricular haemodynamics by means of CFD simulations, International Journal for Numerical Methods in Biomedical Engineering, 34 (2018), pp. 1–14.
  * [9] J. Guzmán and M. Neilan, Inf-sup stable finite elements on barycentric refinements producing divergence–free approximations in arbitrary dimensions, SIAM Journal on Numerical Analysis, 56 (2018), pp. 2826–2844.
  * [10] T. Hughes, G. Engel, L. Mazzei, and M. Larson, The continuous Galerkin method is locally conservative, Journal of Computational Physics, 163 (2000), pp. 467–488.
  * [11] T. Hughes and G. Wells, Conservation properties for the Galerkin and stabilized forms of the advection-diffusion and incompressible Navier-Stokes equations, Computer Methods in Applied Mechanics and Engineering, 194 (2005), pp. 1141–1159.
  * [12] S. Ingimarson, M. Neda, L. Rebholz, J. Reyes, and A. Vu, Improved long time accuracy for projection methods for Navier-Stokes equations using EMAC formulation, International Journal of Numerical Analysis and Modeling, 20 (2023), pp. 176–198.
  * [13] V. John, Reference values for drag and lift of a two dimensional time-dependent flow around a cylinder, International Journal for Numerical Methods in Fluids, 44 (2004), pp. 777–788.
  * [14] V. Kumar, A. Miro, O. Lehmkuhl, and U. Piomelli, Flow separation in airfoils with rough leading edges, AIAA Journal, 61 (2023), pp. 2035–2047.
  * [15] O. Lehmkuhl, G. Houzeaux, H. Owen, G. Chrysokentis, and I. Rodriguez, A low-dissipation finite element scheme for scale resolving simulations of turbulent flows, Journal of Computational Physics, 390 (2019), pp. 51–65.
  * [16] O. Lehmkuhl, U. Piomelli, and G. Houzeaux, On the extension of the integral length-scale approximation model to complex geometries, International Journal of Heat and Fluid Flow, 78 (2019), pp. 1–12.
  * [17] R. J. LeVeque, Finite volume methods for hyperbolic problems, vol. 31, Cambridge university press, 2002.
  * [18] J. Liu and W. Wang, Energy and helicity preserving schemes for hydro and magnetohydro-dynamics flows with symmetry, J. Comput. Phys., 200 (2004), pp. 8–33.
  * [19] M. Olshanskii and L. Rebholz, A note on helicity balance of the Galerkin method for the 3D Navier–Stokes equations, Computer Methods in Applied Mechanics and Engineering, 199 (2010), pp. 1032–1035.
  * [20] M. Olshanskii and L. Rebholz, Longer time accuracy for incompressible Navier-Stokes simulations with the EMAC formulation, Computer Methods in Applied Mechanics and Engineering, 372 (2020), p. 113369.
  * [21] A. Palha and M. Gerritsma, A mass, energy, enstrophy and vorticity conserving (MEEVC) mimetic spectral element discretization for the 2D incompressible Navier–Stokes equations, Journal of Computational Physics, 328 (2017), pp. 200–220.
  * [22] D. Pastrana, J. Cajas, O. Lehmkuhl, I. Rodríguez, and G. Houzeaux, Large-eddy simulations of the vortex-induced vibration of a low mass ratio two-degree-of-freedom circular cylinder at subcritical Reynolds numbers, Computers and Fluids, 173 (2018), pp. 118–132.
  * [23] F. Sacco, B. Paun, O. Lehmkuhl, T. Iles, P. Iaizzo, G. Houzeaux, M. Vzaquex, C. Butakoff, and J. Aguado-Sierra, Left ventricular trabeculations decrease the wall shear stress and increase the intra-ventricular pressure drop in CFD simulations, Frontiers in Physiology, 9 (2018), pp. 1–15.
  * [24] M. Sch$\accentset{\diamond}{\mbox{a}}$fer and S. Turek, The benchmark problem ‘flow around a cylinder’ flow simulation with high performance computers II, in E.H. Hirschel (Ed.), Notes on Numerical Fluid Mechanics, 52, Braunschweig, Vieweg (1996), pp. 547–566.
  * [25] P. Schroeder, V. John, P. Lederer, C. Lehrenfeld, G. Lube, and J. Schoberl, On reference solutions and the sensitivity of the 2D Kelvin-Helmholtz instability problem, Computers and Mathematics with Applications, 77 (2019), pp. 1010–1028.
  * [26] C. Sorgentone, S. L. Cognata, and J. Nordstrom, A new high order energy and enstrophy conserving Arakawa-like Jacobian differential operator, Journal of Computational Physics, 301 (2015), pp. 167–177.
  * [27] P. Varela, P. Suarez, F. Alcantara-Avila, A. Miro, J. Rabault, B. Font, L. Garcia-Cuevas, O. Lehmkuhl, and R. Vinuesa, Deep reinforcement learning for flow control exploits different physics for increasing Reynolds number regimes, Actuators, 11 (2022), p. 359.
  * [28] M. Vazquez, G. Houzeaux, S. Koric, A. Artigues, J. Aguado-Sierra, R. Aris, D. Mira, H. Calmet, F. Cucchietti, H. Owen, A. Taha, E. Burness, J. Cela, and M. Valero, Alya: Multiphysics engineering simulation toward exascale, Journal of Computational Science, 14 (2016), pp. 15–27.