# A fully conservative and shift-invariant formulation for Galerkin
discretizations of incompressible variable density flow

Lukas Lundgren<EMAIL_ADDRESS>Murtazo Nazarov
<EMAIL_ADDRESS>Department of Information Technology, Division of
Scientific Computing, Uppsala University, Sweden

###### Abstract

This paper introduces a formulation of the variable density incompressible
Navier-Stokes equations by modifying the nonlinear terms in a consistent way.
For Galerkin discretizations, the formulation leads to discrete conservation
of mass, squared density, momentum, angular momentum and kinetic energy
without the divergence-free constraint being strongly enforced. In addition to
favorable conservation properties, the formulation is shown to make the
density field invariant to global shifts. The effect of viscous
regularizations on conservation properties is also investigated. Numerical
tests validate the theory developed in this work. The new formulation shows
superior performance compared to other formulations from the literature, both
in terms of accuracy for smooth problems and in terms of robustness.

###### keywords:

incompressible variable density flow, Navier-Stokes equations, conservation,
EMAC formulation, viscous regularization, structure-preserving discretization

††journal: Journal of Computational Physics

##  1 Introduction.

The simulation of incompressible flow is important for industrial and
scientific applications when modeling subsonic flow. In this work, we consider
incompressible flow when the density field is variable, which arises when
multiple fluids interact. A pair of fluids can remain immiscible as seen in
oil/water, gas/liquids and so on. Alternatively, fluids can diffuse and mix,
as observed in pairs of gases and certain pairs of liquids such as salt
water/fresh water, alcohol/water and more. These phenomena can be modeled
using the incompressible Navier-Stokes equations with variable density, along
with an advection equation for the density field. Complexity is often added to
the model by introducing additional terms to account for factors such as
gravity, chemical reactions and surface tension. At a continuous level, when
viscosity effects are neglected, these equations can be shown to conserve
mass, squared density, momentum, angular momentum and kinetic energy. The
density field can also be shown to be shift-invariant which means that density
remains invariant if it is shifted by a constant. These properties are
fundamental to the flow model and are closely linked to the so-called
divergence-free constraint which is pointwise satisfied. It is favorable for
numerical approximations of variable density flow to fulfill as many of these
physical properties as possible.

Unfortunately, numerical approximations of the incompressible Navier-Stokes
equations are rarely pointwise divergence-free [20] which leads to issues with
conservation. In particular, most Galerkin methods are known to satisfy the
divergence-free condition only weakly and much work has been performed to
improve the conservation properties of these methods. Recently, Charnyi et al.
[5] showed that, by consistently modifying the nonlinear term in the constant
density Navier-Stokes equations, a fully conservative Galerkin method could be
derived. The novelty of the formulation is that a modified pressure is solved
for and the formulation has been named the energy, momentum and angular
momentum conserving (EMAC) formulation. The EMAC formulation has later been
extended to be more efficient by Charnyi et al. [6], further analyzed by
Olshanskii and Rebholz [25] and extended to projection methods by Ingimarson
et al. [16].

Unlike constant density flow, variable density flow has additional conserved
properties and no fully conservative formulation is currently known in the
literature. The current state-of-the-art formulations for variable density
flow are based on skew-symmetric formulations that can be shown to conserve
squared density and kinetic energy [13] or a formulation that can be shown to
conserve mass, kinetic energy, momentum and angular momentum [30, 22]. In this
work, we extend the modified pressure technique [5] to variable density flow
which leads to a new formulation that, when discretized by a Galerkin method,
is shift-invariant and also conserves mass, squared density, kinetic energy,
momentum and angular momentum. The analysis is first performed on the inviscid
problem, i.e., the incompressible Euler equations with variable density.
Later, since viscous regularizations affect the conservation properties of the
model, we also investigate this.

The theoretical findings are summarized in Theorems 3.1 and 4.1, and these are
verified using numerical experiments. The results show that the new
formulation leads to improved accuracy and robustness compared to existing
formulations in the variable density literature. Furthermore, the new
formulation can easily be implemented into existing Galerkin-based variable
density Navier-Stokes solvers.

The paper is organized the following way: In Section 2 we give the governing
equations, present our notation and explain the properties of the governing
equations that we later want our proposed numerical method to mimic. Then, in
Section 3 we introduce an alternative formulation of the governing equations
which, when discretized using a Galerkin method, leads to improved
conservation properties in the semi-discrete form. In Section 4 we present a
family of viscous regularizations that lead to various conserved properties of
the model. In Section 5 a conservative second-order accurate time
discretization is presented which leads to a fully discrete method. In Section
6 we perform numerical validations. In Section 7 we give concluding remarks.

##  2 Preliminaries.

In this section, we introduce the governing equations that model variable
density flow. We start with the inviscid case, i.e., we only consider the
incompressible Euler equations with variable density. We also discuss the
properties of the governing equations and present our notation.

###  2.1 The inviscid model problem.

We first consider the incompressible Euler equations with variable density in
a domain $\Omega\subset\mathbb{R}^{d}$ where $d$ is the dimension and with
finite time interval $[0,T]$

$\displaystyle{\partial}_{t}\rho+{\boldsymbol{u}}{\cdot}\nabla\rho$
$\displaystyle=0,$ (2.1)
$\displaystyle{\partial}_{t}{\boldsymbol{m}}+{\boldsymbol{u}}{\cdot}\nabla{\boldsymbol{m}}+\nabla
p$ $\displaystyle={\boldsymbol{f}},\quad$
$\displaystyle({\boldsymbol{x}},t)\in\Omega{\times}(0,T],$
$\displaystyle\nabla{\cdot\,}{\boldsymbol{u}}$ $\displaystyle=0,$
$\displaystyle{\boldsymbol{u}}({\boldsymbol{x}},0)$
$\displaystyle={\boldsymbol{u}}_{0}({\boldsymbol{x}}),\quad$
$\displaystyle{\boldsymbol{x}}\in\Omega,$
$\displaystyle\rho({\boldsymbol{x}},0)$
$\displaystyle=\rho_{0}({\boldsymbol{x}}),$

where the density $\rho>0$, the velocity field ${\boldsymbol{u}}$ and the
pressure $p$ are the unknowns. We define
${\boldsymbol{m}}\mathrel{\mathop{\mathchar 58\relax}}=\rho{\boldsymbol{u}}$
as the momentum vector, ${\boldsymbol{f}}({\boldsymbol{x}},t)$ represents an
external force, $\rho_{0}({\boldsymbol{x}})$,
${\boldsymbol{u}}_{0}({\boldsymbol{x}})$ are initial conditions for density
and velocity. We assume that the governing equations are supplied with well-
posed boundary conditions. To simplify the analysis we set
${\boldsymbol{f}}=0$ and consider periodic boundary conditions.

######  Remark 2.1.

It is possible to write the governing equations in primitive form, i.e.,
$\rho{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\partial}_{t}}{\boldsymbol{u}}+{\boldsymbol{m}}{\cdot}\nabla{\boldsymbol{u}}$
instead of
${\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\partial}_{t}}{\boldsymbol{m}}+{\boldsymbol{u}}{\cdot}\nabla{\boldsymbol{m}}$.
All of the results presented in this manuscript can straightforwardly be
extended to this form. See A for more details.

###  2.2 Notation and useful identities.

Let $\rho,w\in H^{1}(\Omega)$ and
${\boldsymbol{u}},{\boldsymbol{w}},{\boldsymbol{v}}\in{\boldsymbol{H}}^{1}(\Omega)\mathrel{\mathop{\mathchar
58\relax}}=[H^{1}(\Omega)]^{d}$. We note that
$\rho{\boldsymbol{u}}={\boldsymbol{m}}\in[H^{1}(\Omega)]^{d}$. In this work
the $L^{2}$ inner product is written as $({\cdot},{\cdot})$ with associated
$L^{2}$ norms $\|\cdot\|$ and we define
$(\nabla{\boldsymbol{u}})\mathrel{\mathop{\mathchar
58\relax}}=\partial_{x_{i}}u_{j}$ and
$({\boldsymbol{u}}\otimes{\boldsymbol{v}})_{ij}=u_{i}v_{j}$. We also define a
trilinear form

$b({\boldsymbol{u}},{\boldsymbol{v}},{\boldsymbol{w}})\mathrel{\mathop{\mathchar
58\relax}}=({\boldsymbol{u}}{\cdot}\nabla{\boldsymbol{v}},{\boldsymbol{w}})=\left((\nabla{\boldsymbol{v}})^{\top}{\boldsymbol{u}},{\boldsymbol{w}}\right)=((\nabla{\boldsymbol{v}}){\boldsymbol{w}},{\boldsymbol{u}}).$
(2.2)

The following identities hold:

$\displaystyle
b({\boldsymbol{u}},{\boldsymbol{m}},{\boldsymbol{v}})=b({\boldsymbol{m}},{\boldsymbol{u}},{\boldsymbol{v}})+({\boldsymbol{u}}{\cdot}\nabla\rho,{\boldsymbol{u}}{\cdot}{\boldsymbol{v}}),$
(2.3) $\displaystyle
b({\boldsymbol{u}},{\boldsymbol{u}},{\boldsymbol{m}})=b({\boldsymbol{m}},{\boldsymbol{u}},{\boldsymbol{u}}),$
(2.4)

where (2.3) follows from
${\boldsymbol{u}}{\cdot}\nabla(\rho{\boldsymbol{u}})=\rho{\boldsymbol{u}}{\cdot}\nabla{\boldsymbol{u}}+({\boldsymbol{u}}{\cdot}\nabla\rho){\boldsymbol{u}}$
and (2.4) follows from the definition of the trillinear form (2.2).

Since the domain is periodic the following relations hold due to integration
by parts

$\displaystyle
b({\boldsymbol{u}},{\boldsymbol{v}},{\boldsymbol{w}})=-b({\boldsymbol{u}},{\boldsymbol{w}},{\boldsymbol{v}})-((\nabla{\cdot\,}{\boldsymbol{u}}){\boldsymbol{v}},{\boldsymbol{w}}),$
(2.5) $\displaystyle
b({\boldsymbol{u}},{\boldsymbol{w}},{\boldsymbol{w}})=-\frac{1}{2}((\nabla{\cdot\,}{\boldsymbol{u}}){\boldsymbol{w}},{\boldsymbol{w}}),$
(2.6)
$\displaystyle({\boldsymbol{u}}{\cdot}\nabla\rho,w)=-({\boldsymbol{u}}{\cdot}\nabla
w,\rho)-((\nabla{\cdot\,}{\boldsymbol{u}})\rho,w),$ (2.7)
$\displaystyle({\boldsymbol{u}}{\cdot}\nabla\rho,\rho)=-\frac{1}{2}((\nabla{\cdot\,}{\boldsymbol{u}})\rho,\rho).$
(2.8)

We also define ${\boldsymbol{e}}_{i}$ to be a unit vector pointing in the
$i$th cordinate direction, i.e., in 3D ${\boldsymbol{e}}_{1}=(1,0,0)^{\top}$,
${\boldsymbol{e}}_{2}=(0,1,0)^{\top}$ and
${\boldsymbol{e}}_{3}=(0,0,1)^{\top}$. We define
${\boldsymbol{\phi}}_{i}\coloneqq{\boldsymbol{x}}{\times}{\boldsymbol{e}}_{i}$
and note that ${\boldsymbol{\phi}}_{i}$ has the property that
$\nabla{\cdot\,}{\boldsymbol{\phi}}_{i}=0$ and
$\nabla{\boldsymbol{\phi}}_{i}+(\nabla{\boldsymbol{\phi}}_{i})^{\top}=0$. To
get a correct definition of the cross product in 2D, the last component of all
the vectors is extended by 0. The contraction operator : , is defined as

$A\mathrel{\mathop{\mathchar 58\relax}}B=\sum_{i,j=1}^{d}A_{ij}B_{ij},$ (2.9)

for any two $d\times d$ matrices $A$ and $B$.

####  2.2.1 Finite element preliminaries.

We denote a shape regular computational mesh by $\mathcal{T}_{h}$ which is a
triangulation of $\Omega$ into a finite number of disjoint elements $K$. The
finite element spaces we use for the density, velocity and pressure are
respectively given by

$\begin{split}&{\mathcal{M}}\mathrel{\mathop{\mathchar
58\relax}}=\\{w({\boldsymbol{x}})\mathrel{\mathop{\mathchar
58\relax}}w\in\mathcal{C}^{0}(\Omega),w|_{K}\in{\mathbb{P}}_{k_{\rho}},\forall
K\in\mathcal{T}_{h}\\},\\\ &{{\mathbfcal V}}\mathrel{\mathop{\mathchar
58\relax}}=\Big{[}\\{v({\boldsymbol{x}})\mathrel{\mathop{\mathchar
58\relax}}v\in\mathcal{C}^{0}(\Omega),v|_{K}\in{\mathbb{P}}_{k_{\boldsymbol{u}}},\forall
K\in\mathcal{T}_{h}\\}\Big{]}^{d},\\\ &{\mathcal{Q}}\mathrel{\mathop{\mathchar
58\relax}}=\left\\{q({\boldsymbol{x}})\mathrel{\mathop{\mathchar
58\relax}}q\in\mathcal{C}^{0}(\Omega),q|_{K}\in{\mathbb{P}}_{k_{P}},\forall
K\in\mathcal{T}_{h},\int_{\Omega}q\,\mathrm{d}{\boldsymbol{x}}=0\right\\},\end{split}$
(2.10)

where
${\mathbb{P}}_{k_{\rho}},{\mathbb{P}}_{k_{\boldsymbol{u}}},{\mathbb{P}}_{k_{P}}$
are the set of multivariate polynomials of total degree at most
$k_{\rho},k_{\boldsymbol{u}},k_{P}\geq 1$ defined over $K$. It is well-known
that to satisfy the so-called inf-sup condition [9] we require
$k_{P}<k_{\boldsymbol{u}}$. The corresponding discrete inner products are
defined as

$\begin{split}&(v,w)\mathrel{\mathop{\mathchar
58\relax}}=\sum_{K\in\mathcal{T}_{h}}\int_{K}vw\,\mathrm{d}{\boldsymbol{x}},\quad({\boldsymbol{v}},{\boldsymbol{w}})\mathrel{\mathop{\mathchar
58\relax}}=\sum_{K\in\mathcal{T}_{h}}\int_{K}{\boldsymbol{v}}\cdot{\boldsymbol{w}}\,\mathrm{d}{\boldsymbol{x}},\quad(\nabla{\boldsymbol{v}},\nabla{\boldsymbol{w}})\mathrel{\mathop{\mathchar
58\relax}}=\sum_{K\in\mathcal{T}_{h}}\int_{K}\nabla{\boldsymbol{v}}\mathrel{\mathop{\mathchar
58\relax}}\nabla{\boldsymbol{w}}\,\mathrm{d}{\boldsymbol{x}},\end{split}$
(2.11)

with associated $L^{2}$ norms $\|\cdot\|$. Discrete and continuous inner
products are used interchangeably for brevity. The discrete mesh size function
$h({\boldsymbol{x}})$ is constructed as the following continuous piecewise
linear function: for every nodal point $N_{i}$

$h(N_{i})\mathrel{\mathop{\mathchar
58\relax}}=\min(h_{K}/\max({\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}k_{\boldsymbol{u}}},k_{\rho}),\quad\forall
K\in\mathcal{T}_{h}\text{ such that }N_{i}\in K)$ (2.12)

where $h_{K}$ is the smallest edge of $K$. To facilitate the analysis, we
define the following finite element space

${\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\bar{{\mathcal{M}}}_{s}\mathrel{\mathop{\mathchar
58\relax}}=\\{\overline{w}({\boldsymbol{x}})\mathrel{\mathop{\mathchar
58\relax}}\overline{w}\in
L_{2}(\Omega),\overline{w}|_{K}\in{\mathbb{P}}_{s},\forall
K\in\mathcal{T}_{h}\\},}$ (2.13)

where we note that ${\mathcal{M}}\subset\bar{{\mathcal{M}}}_{s}$ provided that
$k_{\rho}\leq s$.

###  2.3 Properties of the governing equations.

The governing equations can be shown to conserve mass, squared density,
kinetic energy, momentum and angular momentum. These are defined as

Mass $\displaystyle\int_{\Omega}\rho\,\mathrm{d}{\boldsymbol{x}};\quad$
Kinetic energy
$\displaystyle\frac{1}{2}\int_{\Omega}\rho{\boldsymbol{u}}{\cdot}{\boldsymbol{u}}\,\mathrm{d}{\boldsymbol{x}};\quad$
Angular momentum
$\displaystyle\int_{\Omega}{\boldsymbol{m}}{\times}{\boldsymbol{x}}\,\mathrm{d}{\boldsymbol{x}};$
(2.14) Momentum
$\displaystyle\int_{\Omega}{\boldsymbol{m}}\,\mathrm{d}{\boldsymbol{x}};\quad$
Squared density
$\displaystyle\frac{1}{2}\int_{\Omega}\rho^{2}\,\mathrm{d}{\boldsymbol{x}}.$

Another property that is rarely (if at all) mentioned in the variable density
flow literature is that the mass equation is invariant to shifts in the
density field. A more precise definition is given below:

######  Definition 2.1.

A differential operator ${\mathcal{L}}$ is shift invariant if
${\mathcal{L}}(\rho+c)={\mathcal{L}}(\rho)$ for any constant $c$.

######  Example 2.1.

The mass equation is shift-invariant since

${\partial}_{t}\rho+{\boldsymbol{u}}{\cdot}\nabla\rho=0\quad\Leftrightarrow\quad{\partial}_{t}(\rho+c)+{\boldsymbol{u}}{\cdot}\nabla(\rho+c)=0.$
(2.15)

##  3 Properties when $\nabla{\cdot\,}{\boldsymbol{u}}\neq 0.$

For most numerical methods, $\nabla{\cdot\,}{\boldsymbol{u}}=0$ does not hold
pointwise at the discrete level. This is also true for most Galerkin methods
where the numerical solution usually satisfies

$(\nabla{\cdot\,}{\boldsymbol{u}},q)=0,\quad\forall
q\in{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color<EMAIL_ADDRESS>(3.1)

We note that there are some exceptions to this. In particular, some exotic
elements are pointwise divergence-free such as iso-geometric B-splines [7],
Raviart-Thomas elements [27, 8], Scott-Vogelius elements [28] and so on, but
some of these are not included in all FEM software libraries and some
introduce constraints on the computational mesh and polynomial degree.
Additionally, not all divergence-free elements, such as [27, 8], are
$H^{1}$-conforming which is a favorable property to have.

In this work, we analyze many variations of the nonlinear terms. Depending on
how these are chosen, favorable conservation properties can be obtained which
we investigate thoroughly in this work. To this end, we write the governing
equations (2.1) in a more general form

$\displaystyle{\partial}_{t}\rho+{\boldsymbol{u}}{\cdot}\nabla\rho+\alpha_{\rho}(\nabla{\cdot\,}{\boldsymbol{u}})\rho$
$\displaystyle=0,$ (3.2)
$\displaystyle{\partial}_{t}{\boldsymbol{m}}+{\boldsymbol{u}}{\cdot}\nabla{\boldsymbol{m}}+\alpha_{\boldsymbol{m}}(\nabla{\cdot\,}{\boldsymbol{u}}){\boldsymbol{m}}+\nabla
P+\alpha_{P}(\nabla{\boldsymbol{u}}){\boldsymbol{m}}+\alpha_{P}(\nabla{\boldsymbol{m}}){\boldsymbol{u}}$
$\displaystyle=0,$

where $P=p-\alpha_{P}{\boldsymbol{u}}{\cdot}{\boldsymbol{m}}$ is a modified
pressure and
$\alpha_{\rho},\alpha_{\boldsymbol{m}},\alpha_{P}\in{\mathbb{R}}$. The terms
involving $\nabla{\cdot\,}{\boldsymbol{u}}$ are consistent with the governing
equations and have previously been used for variable density flow in different
variations to derive numerical methods, see e.g., [13, 26, 3, 1, 30]. The
modified pressure technique ($\alpha_{P}\neq 0$) has been successfully applied
for constant density flow [5] but has currently not been applied in the
variable density context. Depending on how
$\alpha_{\rho},\alpha_{\boldsymbol{m}},\alpha_{P}$ are chosen different
conservation properties are obtained. By using that
$(\nabla{\boldsymbol{m}}){\boldsymbol{u}}=(\nabla(\rho{\boldsymbol{u}})){\boldsymbol{u}}=(\nabla\rho\otimes{\boldsymbol{u}}){\boldsymbol{u}}+(\nabla{\boldsymbol{u}}){\boldsymbol{u}}\rho=\nabla\rho({\boldsymbol{u}}{\cdot}{\boldsymbol{u}})+(\nabla{\boldsymbol{u}}){\boldsymbol{m}}$,
the system (3.2) can equivalently be expressed as

$\displaystyle{\partial}_{t}\rho+{\boldsymbol{u}}{\cdot}\nabla\rho+\alpha_{\rho}(\nabla{\cdot\,}{\boldsymbol{u}})\rho$
$\displaystyle=0,$ (3.3)
$\displaystyle{\partial}_{t}{\boldsymbol{m}}+{\boldsymbol{u}}{\cdot}\nabla{\boldsymbol{m}}+\alpha_{\boldsymbol{m}}(\nabla{\cdot\,}{\boldsymbol{u}}){\boldsymbol{m}}+\nabla
P+\left(\alpha_{P}+\frac{1}{2}\right)(\nabla{\boldsymbol{u}}){\boldsymbol{m}}$
$\displaystyle+\left(\alpha_{P}-\frac{1}{2}\right)(\nabla{\boldsymbol{m}}){\boldsymbol{u}}+\frac{1}{2}(({\boldsymbol{u}}{\cdot}{\boldsymbol{u}})\nabla\rho-\alpha_{\rho}\nabla(\rho{\boldsymbol{u}}{\cdot}{\boldsymbol{u}}))+\frac{1}{2}\alpha_{\rho}\nabla(\rho{\boldsymbol{u}}{\cdot}{\boldsymbol{u}})$
$\displaystyle=0,$

Lastly, we introduce $\bar{\rho}\in{\mathbb{R}}$ which, if properly chosen,
will make the system below shift-invariant

$\displaystyle{\partial}_{t}\rho+{\boldsymbol{u}}{\cdot}\nabla\rho+\alpha_{\rho}(\nabla{\cdot\,}{\boldsymbol{u}})(\rho-\bar{\rho})$
$\displaystyle=0,$ (3.4)
$\displaystyle{\partial}_{t}{\boldsymbol{m}}+{\boldsymbol{u}}{\cdot}\nabla{\boldsymbol{m}}+\alpha_{\boldsymbol{m}}(\nabla{\cdot\,}{\boldsymbol{u}}){\boldsymbol{m}}+\nabla
P+\left(\alpha_{P}+\frac{1}{2}\right)(\nabla{\boldsymbol{u}}){\boldsymbol{m}}$
$\displaystyle+\left(\alpha_{P}-\frac{1}{2}\right)(\nabla{\boldsymbol{m}}){\boldsymbol{u}}+\frac{1}{2}(\nabla\rho({\boldsymbol{u}}{\cdot}{\boldsymbol{u}})-\alpha_{\rho}\nabla((\rho-\bar{\rho}){\boldsymbol{u}}{\cdot}{\boldsymbol{u}}))+\frac{1}{2}\alpha_{\rho}\nabla(\rho{\boldsymbol{u}}{\cdot}{\boldsymbol{u}})$
$\displaystyle=0,$

where
$P=p-\alpha_{P}\rho{\boldsymbol{u}}{\cdot}{\boldsymbol{u}}-\frac{1}{2}\alpha_{\rho}\bar{\rho}{\boldsymbol{u}}{\cdot}{\boldsymbol{u}}$
is a modified pressure which now includes $\bar{\rho}$. To the best of the
authors’ knowledge, using $\bar{\rho}$ to achieve shift-invariant formulations
is new to the literature for both variable density flow and scalar transport
problems in general.

######  Remark 3.1.

In this work, we provide strong numerical evidence supporting that shift
invariance is a valuable property that enhances both the robustness and
accuracy of variable density flow simulations. If a numerical method lacks
this property, it implies that the method would react artificially differently
when the density is shifted. In simpler terms, such a method would produce
larger discretization errors in its results for an initial condition like
$\rho_{0}({\boldsymbol{x}})+100$ compared to $\rho_{0}({\boldsymbol{x}})$.

###  3.1 Galerkin finite element approximation.

The finite element approximation of (3.4) is derived by testing the mass
equation with $w$, the momentum equations with ${\boldsymbol{v}}$ and the
divergence-free constraint with $q$. The finite element method reads: Find
$(\rho,{\boldsymbol{u}},P)\in{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\mathcal{M}}{\times}\mathbf{{\mathbfcal
V}}{\times}{\mathcal{Q}}}$ such that

$\displaystyle({\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\partial}_{t}}\rho,w)+({\boldsymbol{u}}{\cdot}\nabla\rho,w)+\alpha_{\rho}((\nabla{\cdot\,}{\boldsymbol{u}})(\rho-\bar{\rho}),w)$
$\displaystyle=0,\quad$ $\displaystyle\forall w\in{\mathcal{M}},$ (3.5)
$\displaystyle({\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\partial}_{t}}{\boldsymbol{m}},{\boldsymbol{v}})+b({\boldsymbol{u}},{\boldsymbol{m}},{\boldsymbol{v}})+\alpha_{\boldsymbol{m}}((\nabla{\cdot\,}{\boldsymbol{u}}){\boldsymbol{m}},{\boldsymbol{v}})+(\nabla
P,{\boldsymbol{v}})$
$\displaystyle+\left(\alpha_{P}-\frac{1}{2}\right)b({\boldsymbol{v}},{\boldsymbol{m}},{\boldsymbol{u}})+\left(\alpha_{P}+\frac{1}{2}\right)b({\boldsymbol{v}},{\boldsymbol{u}},{\boldsymbol{m}})+\frac{1}{2}\alpha_{\rho}(\nabla({\boldsymbol{m}}{\cdot}{\boldsymbol{u}}),{\boldsymbol{v}})$
$\displaystyle+\frac{1}{2}\left(\left({\boldsymbol{v}}{\cdot}\nabla\rho,{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\boldsymbol{u}}{\cdot}{\boldsymbol{u}}}\right)-\alpha_{\rho}\left(\nabla\left((\rho-\bar{\rho}){\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\boldsymbol{u}}{\cdot}{\boldsymbol{u}}}\right),{\boldsymbol{v}}\right)\right)$
$\displaystyle=0,\quad$
$\displaystyle\forall{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\boldsymbol{v}}\in{\mathbfcal
V}},$ $\displaystyle(\nabla{\cdot\,}{\boldsymbol{u}},q)$
$\displaystyle=0,\quad$ $\displaystyle\forall q\in{\mathcal{Q}}.$

###  3.2 Semi-discrete properties of the Galerkin method.

We summarize the properties of the semi-discrete method (3.5) in Table 1 and
Theorem 3.1 which is the main result of this work. In the constant density
case, there exists a unique parameter choice such that kinetic energy,
momentum and angular momentum are all conserved which has been named the EMAC
formulation [5]. In the variable density context, we require additional
constraints to achieve all properties in Table 1. All properties are achieved
by setting
$\alpha_{\rho}=\frac{1}{2},\alpha_{\boldsymbol{m}}=1,\alpha_{P}=0.25,\bar{\rho}=\frac{1}{|\Omega|}\int_{\Omega}\rho\,\mathrm{d}{\boldsymbol{x}}$
and by choosing the polynomial degree of density to be less than or equal to
that of pressure ($k_{\rho}\leq k_{P}$). This leads to a formulation that is
shift-invariant and mass, kinetic energy, squared density, momentum and
angular momentum conserving (SI-MEDMAC). If other parameter choices
$\alpha_{\rho},\alpha_{\boldsymbol{m}},\alpha_{P},\bar{\rho},k_{\rho}$ are
considered, we end up with common formulations from the variable density flow
literature and several new ones, see Table 2. We follow a similar naming
convention for the other formulations in Table 2, i.e., MEMAC is mass, kinetic
energy, momentum, angular momentum conserving, EDMAC is kinetic energy,
squared density, momentum, angular momentum conserving and so on. All
formulations with $\alpha_{\rho}=0$ are, for a lack of a better term, named
locally shift-invariant (LSI). If the domain is stretched, these formulations
are invariant to this change but the other formulations are not. The
convective form is the form that is obtained when (2.1) is discretized
directly without any modifications.

Table 1: Semi-discrete properties of the method (3.5) when ${\boldsymbol{f}}=0$ for the inviscid case when $\nabla{\cdot\,}{\boldsymbol{u}}\neq 0$. Property | Condition
---|---
$\partial_{t}\int_{\Omega}\rho{\boldsymbol{u}}{\cdot}{\boldsymbol{u}}\,\mathrm{d}{\boldsymbol{x}}$ = 0 | $\alpha_{\boldsymbol{m}}-\alpha_{P}-\alpha_{\rho}/2=1/2$
$\partial_{t}\int_{\Omega}{\boldsymbol{m}}\,\mathrm{d}{\boldsymbol{x}}$ = 0 | $\alpha_{\boldsymbol{m}}=1$
$\partial_{t}\int_{\Omega}{\boldsymbol{m}}{\times}{\boldsymbol{x}}\,\mathrm{d}{\boldsymbol{x}}=0$ | $\alpha_{\boldsymbol{m}}=1$
$\partial_{t}\frac{1}{2}\int_{\Omega}\rho^{2}\,\mathrm{d}{\boldsymbol{x}}-\int_{\Omega}\bar{\rho}\partial_{t}\rho\,\mathrm{d}{\boldsymbol{x}}=0$ | $\alpha_{\rho}=1/2$
$\partial_{t}\left(\int_{\Omega}\rho^{2}\,\mathrm{d}{\boldsymbol{x}}-\frac{1}{|\Omega|}\left(\int_{\Omega}\rho\,\mathrm{d}{\boldsymbol{x}}\right)^{2}\right)=0$ |  $\alpha_{\rho}=1/2$ and $\bar{\rho}=\frac{1}{|\Omega|}\int_{\Omega}\rho\,\mathrm{d}{\boldsymbol{x}}$
$\partial_{t}\int_{\Omega}\rho^{2}\,\mathrm{d}{\boldsymbol{x}}=0$ | [ $\alpha_{\rho}=1/2$ and $k_{\rho}\leq k_{P}$ ] or $2k_{\rho}\leq k_{P}$
$\partial_{t}\int_{\Omega}\rho\,\mathrm{d}{\boldsymbol{x}}$ = 0 |  $k_{\rho}\leq k_{P}$ or $\alpha_{\rho}=1$
Shift invariance |  $\alpha_{\rho}=0$ or $\bar{\rho}=\frac{1}{|\Omega|}\int_{\Omega}\rho\,\mathrm{d}{\boldsymbol{x}}$
Table 2: Existing formulations in the literature and newly derived formulations for variable density flow. The formulations are later tested and compared in numerical experiments. Name | $\alpha_{\rho}$ | $\alpha_{\boldsymbol{m}}$ | $\alpha_{P}$ | $\bar{\rho}$ | $k_{\rho}$
---|---|---|---|---|---
LSI-EMAC | 0 | 1 | 0.5 | N/A | $\forall k_{\rho}$
MEMAC [30, 22] | 1 | 1 | 0 | 0 | $\forall k_{\rho}$
EDMAC | 0.5 | 1 | 0.25 | 0 | $\forall k_{\rho}$
SI-MEMAC | 1 | 1 | 0 | $\bar{\rho}=\frac{1}{|\Omega|}\int_{\Omega}\rho\,\mathrm{d}{\boldsymbol{x}}$ | $\forall k_{\rho}$
SI-EDMAC | 0.5 | 1 | 0.25 | $\bar{\rho}=\frac{1}{|\Omega|}\int_{\Omega}\rho\,\mathrm{d}{\boldsymbol{x}}$ | $\forall k_{\rho}$
SI-MEDMAC | 0.5 | 1 | 0.25 | $\bar{\rho}=\frac{1}{|\Omega|}\int_{\Omega}\rho\,\mathrm{d}{\boldsymbol{x}}$ | $k_{\rho}\leq k_{P}$
LSI-EC [21] | 0 | 0.5 | 0 | N/A | $\forall k_{\rho}$
Convective [14, Sec 4] | 0 | 0 | 0 | N/A | $\forall k_{\rho}$

######  Theorem 3.1.

The properties of Table 1 hold provided that the conditions inside it are met.

###### Proof.

We divide the proof into several parts. To simplify the proof we assume a
periodic domain and set ${\boldsymbol{f}}=0$. Extending the proofs to other
boundary conditions such as no-slip (${\boldsymbol{u}}=0|_{{\partial}\Omega}$)
can be done by using the proof technique in [5, Sec 3.1.2] or [16, Sec 3.1.2].
Mass: By using the definition we have that (3.5) satisfies

${\partial}_{t}\int_{\Omega}\rho\,\mathrm{d}{\boldsymbol{x}}=\int_{\Omega}{\partial}_{t}\rho\,\mathrm{d}{\boldsymbol{x}}=-({\boldsymbol{u}},\nabla\rho)-\alpha_{\rho}(\nabla{\cdot\,}{\boldsymbol{u}},\rho-\bar{\rho}).$
(3.6)

Since $\nabla\bar{\rho}=0$, (3.6) is equivalent to

$\begin{split}{\partial}_{t}\int_{\Omega}\rho\,\mathrm{d}{\boldsymbol{x}}=&-({\boldsymbol{u}},\nabla(\rho-\bar{\rho}))-\alpha_{\rho}(\nabla{\cdot\,}{\boldsymbol{u}},\rho-\bar{\rho})\\\
=\
&(\nabla{\cdot\,}{\boldsymbol{u}},\rho-\bar{\rho})-\alpha_{\rho}(\nabla{\cdot\,}{\boldsymbol{u}},\rho-\bar{\rho})\\\
=\ &\ (1-\alpha_{\rho})(\nabla{\cdot\,}{\boldsymbol{u}},\rho-\bar{\rho}),\\\
\end{split}$ (3.7)

where integration by parts was used on the first term. Therefore mass is
conserved if $\alpha_{\rho}=1$. If
$\bar{\rho}=\int_{\Omega}\rho\,\mathrm{d}{\boldsymbol{x}}$ then

$\int_{\Omega}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}(}\rho-\bar{\rho}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0})}\,\mathrm{d}{\boldsymbol{x}}=0,$
(3.8)

which means that $(\rho-\bar{\rho})\in{\mathcal{Q}}$ provided that
$k_{\rho}\leq k_{P}$. Invoking the weak divergence-free condition (3.5) shows
that

$(\nabla{\cdot\,}{\boldsymbol{u}},\rho-\bar{\rho})=0,$ (3.9)

which means that mass is conserved if $k_{\rho}\leq k_{P}$.

Shift-invariance: If $\alpha_{\rho}=0$ then the proof is evident from
Definition 2.1. Otherwise, if
$\bar{\rho}=\frac{1}{|\Omega|}\int_{\Omega}\rho\,\mathrm{d}{\boldsymbol{x}}$,
the following is obtained by shifting $\rho$ by a constant $c$ in (3.4)

$\begin{split}{\partial}_{t}(\rho+c)+{\boldsymbol{u}}{\cdot}\nabla(\rho+c)+\alpha_{\rho}(\nabla{\cdot\,}{\boldsymbol{u}})\left(\rho+c-\frac{1}{|\Omega|}\int_{\Omega}(\rho+c)\,\mathrm{d}{\boldsymbol{x}}\right)&=0,\\\
{\partial}_{t}\rho+{\boldsymbol{u}}{\cdot}\nabla\rho+\alpha_{\rho}(\nabla{\cdot\,}{\boldsymbol{u}})\left(\rho+c-\frac{1}{|\Omega|}\int_{\Omega}(\rho+c)\,\mathrm{d}{\boldsymbol{x}}\right)&=0,\\\
{\partial}_{t}\rho+{\boldsymbol{u}}{\cdot}\nabla\rho+\alpha_{\rho}(\nabla{\cdot\,}{\boldsymbol{u}})\left(\rho-\frac{1}{|\Omega|}\int_{\Omega}\rho\,\mathrm{d}{\boldsymbol{x}}\right)&=0,\\\
{\partial}_{t}\rho+{\boldsymbol{u}}{\cdot}\nabla\rho+\alpha_{\rho}(\nabla{\cdot\,}{\boldsymbol{u}})\left(\rho-\bar{\rho}\right)&=0,\\\
\end{split}$ (3.10)

which shows that the formulation is shift-invariant.
Squared density: By using that $\nabla\bar{\rho}=0$, one can show that (3.5)
satisfies

$\partial_{t}\frac{1}{2}\int_{\Omega}\rho^{2}\,\mathrm{d}{\boldsymbol{x}}-\int_{\Omega}\bar{\rho}\partial_{t}\rho\,\mathrm{d}{\boldsymbol{x}}=({\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\partial}_{t}}\rho,\rho-\bar{\rho})=-({\boldsymbol{u}}{\cdot}\nabla(\rho-\bar{\rho}),\rho-\bar{\rho})-\alpha_{\rho}((\nabla{\cdot\,}{\boldsymbol{u}})(\rho-\bar{\rho}),\rho-\bar{\rho}).$
(3.11)

Using (2.8) inside (3.11) yields

$\partial_{t}\frac{1}{2}\int_{\Omega}\rho^{2}\,\mathrm{d}{\boldsymbol{x}}-\int_{\Omega}\bar{\rho}\partial_{t}\rho\,\mathrm{d}{\boldsymbol{x}}=\left(\frac{1}{2}-\alpha_{\rho}\right)((\nabla{\cdot\,}{\boldsymbol{u}})(\rho-\bar{\rho}),\rho-\bar{\rho}).$
(3.12)

If $\alpha_{\rho}=\frac{1}{2}$ then
$\partial_{t}\frac{1}{2}\int_{\Omega}\rho^{2}\,\mathrm{d}{\boldsymbol{x}}-\int_{\Omega}\bar{\rho}\partial_{t}\rho\,\mathrm{d}{\boldsymbol{x}}=0$.
This shows that line 4 in Table 1 is true.

If
$\bar{\rho}=\frac{1}{|\Omega|}\int_{\Omega}\rho\,\mathrm{d}{\boldsymbol{x}}$,
then we can perform the following simplification

$\begin{split}\partial_{t}\frac{1}{2}\int_{\Omega}\rho^{2}\,\mathrm{d}{\boldsymbol{x}}-\int_{\Omega}\bar{\rho}\partial_{t}\rho\,\mathrm{d}{\boldsymbol{x}}&=\partial_{t}\frac{1}{2}\int_{\Omega}\rho^{2}\,\mathrm{d}{\boldsymbol{x}}-\bar{\rho}\int_{\Omega}\partial_{t}\rho\,\mathrm{d}{\boldsymbol{x}}\\\
&=\partial_{t}\frac{1}{2}\int_{\Omega}\rho^{2}\,\mathrm{d}{\boldsymbol{x}}-\frac{1}{|\Omega|}\left(\int_{\Omega}\rho\,\mathrm{d}{\boldsymbol{x}}\right)\int_{\Omega}\partial_{t}\rho\,\mathrm{d}{\boldsymbol{x}}\\\
&=\partial_{t}\frac{1}{2}\int_{\Omega}\rho^{2}\,\mathrm{d}{\boldsymbol{x}}-\frac{1}{|\Omega|}\left(\int_{\Omega}\rho\,\mathrm{d}{\boldsymbol{x}}\right)\partial_{t}\left(\int_{\Omega}\rho\,\mathrm{d}{\boldsymbol{x}}\right)\\\
&=\partial_{t}\frac{1}{2}\left(\int_{\Omega}\rho^{2}\,\mathrm{d}{\boldsymbol{x}}-\frac{1}{|\Omega|}\left(\int_{\Omega}\rho\,\mathrm{d}{\boldsymbol{x}}\right)^{2}\right),\end{split}$
(3.13)

to show that line 5 in Table 1 is true.

Lastly, we consider line 6 in Table 1. By recalling that
$\partial_{t}\int_{\Omega}\rho=0$ if $k_{\rho}\leq k_{P}$, (3.13) simplifies
to

$\partial_{t}\frac{1}{2}\int_{\Omega}\rho^{2}\,\mathrm{d}{\boldsymbol{x}}-\int_{\Omega}\bar{\rho}\partial_{t}\rho\,\mathrm{d}{\boldsymbol{x}}=\partial_{t}\frac{1}{2}\int_{\Omega}\rho^{2}\,\mathrm{d}{\boldsymbol{x}},$
(3.14)

which shows that squared density is conserved if $\alpha_{\rho}=\frac{1}{2}$
and $k_{\rho}\leq k_{P}$. Next, for $A\in{\mathbb{R}}$ we also have

$\begin{split}\partial_{t}\frac{1}{2}\int_{\Omega}\rho^{2}\,\mathrm{d}{\boldsymbol{x}}-\int_{\Omega}A\partial_{t}\rho\,\mathrm{d}{\boldsymbol{x}}&=({\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\partial}_{t}}\rho,\rho-A)=-({\boldsymbol{u}}{\cdot}\nabla(\rho-A),\rho-A)-\alpha_{\rho}((\nabla{\cdot\,}{\boldsymbol{u}})(\rho-\bar{\rho}),\rho-A)\\\
&=\left(\nabla{\cdot\,}{\boldsymbol{u}},\frac{1}{2}(\rho-A)^{2}-\alpha_{\rho}(\rho-\bar{\rho})(\rho-A)\right).\\\
\end{split}$ (3.15)

There exist two roots of $A$ that satisfies

$\int_{\Omega}\left(\frac{1}{2}(\rho-A)^{2}-\alpha_{\rho}(\rho-\bar{\rho})(\rho-A)\right)\,\mathrm{d}{\boldsymbol{x}}=0.$
(3.16)

We only need one $A$ to satisfy (3.16) and one of the roots (found
symbolically using Matlab) are

$\begin{split}A=|\Omega|^{-1}\left(\alpha_{\rho}|\Omega|\bar{\rho}+(1-\alpha_{\rho})\int_{\Omega}\rho\,\mathrm{d}{\boldsymbol{x}}+\left(\alpha_{\rho}^{2}|\Omega|^{2}\bar{\rho}^{2}-2\alpha_{\rho}^{2}|\Omega|\bar{\rho}\int_{\Omega}\rho\,\mathrm{d}{\boldsymbol{x}}+\alpha_{\rho}^{2}\left(\int_{\Omega}\rho\,\mathrm{d}{\boldsymbol{x}}\right)^{2}+2\alpha_{\rho}|\Omega|\int_{\Omega}\rho^{2}\,\mathrm{d}{\boldsymbol{x}}\right.\right.\\\
\left.\left.-2\alpha_{\rho}\left(\int_{\Omega}\rho\,\mathrm{d}{\boldsymbol{x}}\right)^{2}-|\Omega|\int_{\Omega}\rho^{2}\,\mathrm{d}{\boldsymbol{x}}+\left(\int_{\Omega}\rho\,\mathrm{d}{\boldsymbol{x}}\right)^{2}\right)^{1/2}\right).\end{split}$
(3.17)

This means that, if $2k_{\rho}\leq k_{P}$ then there exists an $A$ such that
$\left(\frac{1}{2}(\rho-A)^{2}-\alpha_{\rho}(\rho-\bar{\rho})(\rho-A)\right)\in{\mathcal{Q}}$.
Because of this, we can invoke the weak divergence-free condition (3.5) to
conclude that squared density is conserved.

Kinetic energy: We define the following $L^{2}$-projection. Find
${\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\overline{{\boldsymbol{u}}{\cdot}{\boldsymbol{u}}}}\in{\mathcal{M}}$
such that

${\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\left({\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\overline{{\boldsymbol{u}}{\cdot}{\boldsymbol{u}}}},\overline{w}\right)=\left({\boldsymbol{u}}{\cdot}{\boldsymbol{u}},\overline{w}\right)\quad\forall\overline{w}\in\bar{{\mathcal{M}}}_{k_{\boldsymbol{u}}+k_{\rho}-1}},$
(3.18)

where $\bar{{\mathcal{M}}}_{k_{\boldsymbol{u}}+k_{\rho}-1}$ is defined by
(2.13). Next, using the definition of kinetic energy, the definition of the
$L^{2}$-projection (3.18) and setting
$w={\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\overline{{\boldsymbol{u}}{\cdot}{\boldsymbol{u}}}},{\boldsymbol{v}}={\boldsymbol{u}}$
inside the density and momentum updates (3.5) yields

$\begin{split}\frac{1}{2}{\partial}_{t}\int_{\Omega}\rho{\boldsymbol{u}}{\cdot}{\boldsymbol{u}}\,\mathrm{d}{\boldsymbol{x}}=&\
\frac{1}{2}({\boldsymbol{u}}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\partial}_{t}\rho},{\boldsymbol{u}})+(\rho{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\partial}_{t}}{\boldsymbol{u}},{\boldsymbol{u}})=({\boldsymbol{u}}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\partial}_{t}}\rho,{\boldsymbol{u}})+(\rho{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\partial}_{t}}{\boldsymbol{u}},{\boldsymbol{u}})-\frac{1}{2}({\boldsymbol{u}}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\partial}_{t}\rho},{\boldsymbol{u}})\\\
=&\
({\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\partial}_{t}}{\boldsymbol{m}},{\boldsymbol{u}})-\frac{1}{2}\left({\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\partial}_{t}}\rho,{\boldsymbol{u}}{\cdot}{\boldsymbol{u}}\right)=({\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\partial}_{t}}{\boldsymbol{m}},{\boldsymbol{u}})-\frac{1}{2}\left({\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\partial}_{t}}\rho,{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\overline{{\boldsymbol{u}}{\cdot}{\boldsymbol{u}}}}\right)\\\
=&-b({\boldsymbol{u}},{\boldsymbol{m}},{\boldsymbol{u}})-\alpha_{\boldsymbol{m}}((\nabla{\cdot\,}{\boldsymbol{u}}){\boldsymbol{m}},{\boldsymbol{u}})-(\nabla
P,{\boldsymbol{u}})\\\
&-\left(\alpha_{P}-\frac{1}{2}\right)b({\boldsymbol{u}},{\boldsymbol{m}},{\boldsymbol{u}})-\left(\alpha_{P}+\frac{1}{2}\right)b({\boldsymbol{u}},{\boldsymbol{u}},{\boldsymbol{m}})-\frac{1}{2}\alpha_{\rho}(\nabla({\boldsymbol{m}}{\cdot}{\boldsymbol{u}}),{\boldsymbol{u}})\\\
&-\frac{1}{2}\left(\left({\boldsymbol{u}}{\cdot}\nabla\rho,{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\boldsymbol{u}}{\cdot}{\boldsymbol{u}}}\right)-\alpha_{\rho}\left(\nabla\left((\rho-\bar{\rho}){\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\boldsymbol{u}}{\cdot}{\boldsymbol{u}}}\right),{\boldsymbol{u}}\right)\right)\\\
&+\frac{1}{2}\left({\boldsymbol{u}}{\cdot}\nabla\rho+\alpha_{\rho}(\nabla{\cdot\,}{\boldsymbol{u}})(\rho-\bar{\rho}),{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\overline{{\boldsymbol{u}}{\cdot}{\boldsymbol{u}}}}\right).\end{split}$
(3.19)

Integration by parts on the pressure term and invoking the weak divergence-
free constraint inside (3.5), using (2.5) on
$\left(\alpha_{P}-\frac{1}{2}\right)b({\boldsymbol{u}},{\boldsymbol{m}},{\boldsymbol{u}})$
and integration by parts on
$\alpha_{\rho}\left(\nabla\left((\rho-\bar{\rho}){\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\boldsymbol{u}}{\cdot}{\boldsymbol{u}}}\right),{\boldsymbol{u}}\right)$
yields

$\begin{split}\frac{1}{2}{\partial}_{t}\int_{\Omega}\rho{\boldsymbol{u}}{\cdot}{\boldsymbol{u}}\,\mathrm{d}{\boldsymbol{x}}=&-b({\boldsymbol{u}},{\boldsymbol{m}},{\boldsymbol{u}})+\left(\alpha_{P}-\frac{1}{2}-\alpha_{\boldsymbol{m}}\right)((\nabla{\cdot\,}{\boldsymbol{u}}){\boldsymbol{m}},{\boldsymbol{u}})\\\
&+\left(\alpha_{P}-\frac{1}{2}\right)b({\boldsymbol{u}},{\boldsymbol{u}},{\boldsymbol{m}})-\left(\alpha_{P}+\frac{1}{2}\right)b({\boldsymbol{u}},{\boldsymbol{u}},{\boldsymbol{m}})-\frac{1}{2}\alpha_{\rho}(\nabla({\boldsymbol{m}}{\cdot}{\boldsymbol{u}}),{\boldsymbol{u}})\\\
&-\frac{1}{2}\left(\left({\boldsymbol{u}}{\cdot}\nabla\rho,{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\boldsymbol{u}}{\cdot}{\boldsymbol{u}}}\right)+\alpha_{\rho}\left((\nabla{\cdot\,}{\boldsymbol{u}}),(\rho-\bar{\rho}){\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\boldsymbol{u}}{\cdot}{\boldsymbol{u}}}\right)\right)\\\
&+\frac{1}{2}\left({\boldsymbol{u}}{\cdot}\nabla\rho+\alpha_{\rho}(\nabla{\cdot\,}{\boldsymbol{u}})(\rho-\bar{\rho}),{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\overline{{\boldsymbol{u}}{\cdot}{\boldsymbol{u}}}}\right)\\\
\\\
=&-b({\boldsymbol{u}},{\boldsymbol{m}},{\boldsymbol{u}})+\left(\alpha_{P}-\frac{1}{2}-\alpha_{\boldsymbol{m}}\right)((\nabla{\cdot\,}{\boldsymbol{u}}){\boldsymbol{m}},{\boldsymbol{u}})\\\
&-b({\boldsymbol{u}},{\boldsymbol{u}},{\boldsymbol{m}})-\frac{1}{2}\alpha_{\rho}(\nabla({\boldsymbol{m}}{\cdot}{\boldsymbol{u}}),{\boldsymbol{u}})\\\
&{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}-\frac{1}{2}\left({\boldsymbol{u}}{\cdot}\nabla\rho+\alpha_{\rho}(\nabla{\cdot\,}{\boldsymbol{u}})(\rho-\bar{\rho}),{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\boldsymbol{u}}{\cdot}{\boldsymbol{u}}}\right)}\\\
&{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}+\frac{1}{2}\left({\boldsymbol{u}}{\cdot}\nabla\rho+\alpha_{\rho}(\nabla{\cdot\,}{\boldsymbol{u}})(\rho-\bar{\rho}),{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\overline{{\boldsymbol{u}}{\cdot}{\boldsymbol{u}}}}\right).}\end{split}$
(3.20)

Since
$\left({\boldsymbol{u}}{\cdot}\nabla\rho+\alpha_{\rho}(\nabla{\cdot\,}{\boldsymbol{u}})(\rho-\bar{\rho})\right)\in\bar{{\mathcal{M}}}_{k_{\boldsymbol{u}}+k_{\rho}-1}$,
we can use (3.18) to show that

$\begin{split}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\frac{1}{2}{\partial}_{t}\int_{\Omega}\rho{\boldsymbol{u}}{\cdot}{\boldsymbol{u}}\,\mathrm{d}{\boldsymbol{x}}=}&{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}-b({\boldsymbol{u}},{\boldsymbol{m}},{\boldsymbol{u}})+\left(\alpha_{P}-\frac{1}{2}-\alpha_{\boldsymbol{m}}\right)((\nabla{\cdot\,}{\boldsymbol{u}}){\boldsymbol{m}},{\boldsymbol{u}})}\\\
&{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}-b({\boldsymbol{u}},{\boldsymbol{u}},{\boldsymbol{m}})-\frac{1}{2}\alpha_{\rho}(\nabla({\boldsymbol{m}}{\cdot}{\boldsymbol{u}}),{\boldsymbol{u}}).}\end{split}$
(3.21)

Next, using (2.5) on $b({\boldsymbol{u}},{\boldsymbol{m}},{\boldsymbol{u}})$
and integration by parts on
$\frac{1}{2}\alpha_{\rho}(\nabla({\boldsymbol{m}}{\cdot}{\boldsymbol{u}}),{\boldsymbol{u}})$
yields

$\begin{split}\frac{1}{2}{\partial}_{t}\int_{\Omega}\rho{\boldsymbol{u}}{\cdot}{\boldsymbol{u}}\,\mathrm{d}{\boldsymbol{x}}=&\left(\alpha_{P}+\frac{1}{2}-\alpha_{\boldsymbol{m}}+\frac{1}{2}\alpha_{\rho}\right)((\nabla{\cdot\,}{\boldsymbol{u}}){\boldsymbol{m}},{\boldsymbol{u}}),\\\
\end{split}$ (3.22)

which means that kinetic energy is conserved if
$\alpha_{\boldsymbol{m}}-\alpha_{P}-\alpha_{\rho}/2=1/2$.
Momentum: We set ${\boldsymbol{v}}={\boldsymbol{e}}_{i}$ inside the momentum
equations (3.5) and obtain

$\begin{split}{\partial}_{t}\int_{\Omega}{\boldsymbol{m}}_{i}\,\mathrm{d}{\boldsymbol{x}}=&\
({\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\partial}_{t}}{\boldsymbol{m}},{\boldsymbol{e}}_{i})=-b({\boldsymbol{u}},{\boldsymbol{m}},{\boldsymbol{e}}_{i})-\alpha_{\boldsymbol{m}}((\nabla{\cdot\,}{\boldsymbol{u}}){\boldsymbol{m}},{\boldsymbol{e}}_{i})-(\nabla
P,{\boldsymbol{e}}_{i})\\\
&-\left(\alpha_{P}-\frac{1}{2}\right)b({\boldsymbol{e}}_{i},{\boldsymbol{m}},{\boldsymbol{u}})-\left(\alpha_{P}+\frac{1}{2}\right)b({\boldsymbol{e}}_{i},{\boldsymbol{u}},{\boldsymbol{m}})+\frac{1}{2}\alpha_{\rho}((\nabla{\cdot\,}{\boldsymbol{e}}_{i}){\boldsymbol{m}},{\boldsymbol{u}})\\\
&-\frac{1}{2}\left(\left({\boldsymbol{e}}_{i}{\cdot}\nabla\rho,{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\boldsymbol{u}}{\cdot}{\boldsymbol{u}}}\right)+\alpha_{\rho}\left((\nabla{\cdot\,}{\boldsymbol{e}}_{i})(\rho-\bar{\rho}),{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\boldsymbol{u}}{\cdot}{\boldsymbol{u}}}\right)\right)\\\
\\\
=&-b({\boldsymbol{u}},{\boldsymbol{m}},{\boldsymbol{e}}_{i})-\alpha_{\boldsymbol{m}}((\nabla{\cdot\,}{\boldsymbol{u}}){\boldsymbol{m}},{\boldsymbol{e}}_{i})-(\nabla
P,{\boldsymbol{e}}_{i})\\\
&-\left(\alpha_{P}-\frac{1}{2}\right)b({\boldsymbol{e}}_{i},{\boldsymbol{m}},{\boldsymbol{u}})-\left(\alpha_{P}+\frac{1}{2}\right)b({\boldsymbol{e}}_{i},{\boldsymbol{u}},{\boldsymbol{m}})-\frac{1}{2}\left({\boldsymbol{e}}_{i}{\cdot}\nabla\rho,{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\boldsymbol{u}}{\cdot}{\boldsymbol{u}}}\right),\end{split}$
(3.23)

since ${\boldsymbol{e}}_{i}$ is a constant. Using (2.5) on
$\left(\alpha_{P}-\frac{1}{2}\right)b({\boldsymbol{e}}_{i},{\boldsymbol{m}},{\boldsymbol{u}})$
and $b({\boldsymbol{u}},{\boldsymbol{m}},{\boldsymbol{e}}_{i})$ and
integration by parts on the pressure term yield

$\begin{split}{\partial}_{t}\int_{\Omega}{\boldsymbol{m}}_{i}\,\mathrm{d}{\boldsymbol{x}}=&\
b({\boldsymbol{u}},{\boldsymbol{e}}_{i},{\boldsymbol{m}})+((\nabla{\cdot\,}{\boldsymbol{u}}){\boldsymbol{m}},{\boldsymbol{e}}_{i})-\alpha_{\boldsymbol{m}}((\nabla{\cdot\,}{\boldsymbol{u}}){\boldsymbol{m}},{\boldsymbol{e}}_{i})+(P,\nabla{\cdot\,}{\boldsymbol{e}}_{i})\\\
&+\left(\alpha_{P}-\frac{1}{2}\right)b({\boldsymbol{e}}_{i},{\boldsymbol{u}},{\boldsymbol{m}})-\left(\alpha_{P}+\frac{1}{2}\right)b({\boldsymbol{e}}_{i},{\boldsymbol{u}},{\boldsymbol{m}})-\frac{1}{2}\left({\boldsymbol{e}}_{i}{\cdot}\nabla\rho,{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\boldsymbol{u}}{\cdot}{\boldsymbol{u}}}\right)\\\
=&\
(1-\alpha_{\boldsymbol{m}})((\nabla{\cdot\,}{\boldsymbol{u}}){\boldsymbol{m}},{\boldsymbol{e}}_{i})-b({\boldsymbol{e}}_{i},{\boldsymbol{u}},{\boldsymbol{m}})-\frac{1}{2}\left({\boldsymbol{e}}_{i}{\cdot}\nabla\rho,{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\boldsymbol{u}}{\cdot}{\boldsymbol{u}}}\right)\\\
=&\
(1-\alpha_{\boldsymbol{m}})((\nabla{\cdot\,}{\boldsymbol{u}}){\boldsymbol{m}},{\boldsymbol{e}}_{i})-b(\rho{\boldsymbol{e}}_{i},{\boldsymbol{u}},{\boldsymbol{u}})-\frac{1}{2}\left({\boldsymbol{e}}_{i}{\cdot}\nabla\rho,{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\boldsymbol{u}}{\cdot}{\boldsymbol{u}}}\right)\\\
=&\
(1-\alpha_{\boldsymbol{m}})((\nabla{\cdot\,}{\boldsymbol{u}}){\boldsymbol{m}},{\boldsymbol{e}}_{i})+\frac{1}{2}((\nabla{\cdot\,}(\rho{\boldsymbol{e}}_{i})){\boldsymbol{u}},{\boldsymbol{u}})-\frac{1}{2}\left({\boldsymbol{e}}_{i}{\cdot}\nabla\rho,{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\boldsymbol{u}}{\cdot}{\boldsymbol{u}}}\right)\\\
=&\
(1-\alpha_{\boldsymbol{m}})((\nabla{\cdot\,}{\boldsymbol{u}}){\boldsymbol{m}},{\boldsymbol{e}}_{i}),\end{split}$
(3.24)

where (2.6) was used in the third equality and (2.4) was used in the second
equality. Thus, momentum is conserved if $\alpha_{\boldsymbol{m}}=1$.
Angular momentum: Conservation of angular momentum is derived by setting
${\boldsymbol{v}}={\boldsymbol{\phi}}_{i}$ inside (3.5) to obtain

$\begin{split}{\partial}_{t}\int_{\Omega}({\boldsymbol{m}}{\times}{\boldsymbol{x}})_{i}\,\mathrm{d}{\boldsymbol{x}}=&\
({\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\partial}_{t}}{\boldsymbol{m}},{\boldsymbol{\phi}}_{i})=-b({\boldsymbol{u}},{\boldsymbol{m}},{\boldsymbol{\phi}}_{i})-\alpha_{\boldsymbol{m}}((\nabla{\cdot\,}{\boldsymbol{u}}){\boldsymbol{m}},{\boldsymbol{\phi}}_{i})-(\nabla
P,{\boldsymbol{\phi}}_{i})\\\
&-\left(\alpha_{P}-\frac{1}{2}\right)b({\boldsymbol{\phi}}_{i},{\boldsymbol{m}},{\boldsymbol{u}})-\left(\alpha_{P}+\frac{1}{2}\right)b({\boldsymbol{\phi}}_{i},{\boldsymbol{u}},{\boldsymbol{m}})+\frac{1}{2}\alpha_{\rho}((\nabla{\cdot\,}{\boldsymbol{\phi}}_{i}){\boldsymbol{m}},{\boldsymbol{u}})\\\
&-\frac{1}{2}\left(\left({\boldsymbol{\phi}}_{i}{\cdot}\nabla\rho,{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\boldsymbol{u}}{\cdot}{\boldsymbol{u}}}\right)+\alpha_{\rho}\left((\nabla{\cdot\,}{\boldsymbol{\phi}}_{i})(\rho-\bar{\rho}),{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\boldsymbol{u}}{\cdot}{\boldsymbol{u}}}\right)\right)\\\
\\\
=&-b({\boldsymbol{u}},{\boldsymbol{m}},{\boldsymbol{\phi}}_{i})-\alpha_{\boldsymbol{m}}((\nabla{\cdot\,}{\boldsymbol{u}}){\boldsymbol{m}},{\boldsymbol{\phi}}_{i})-(\nabla
P,{\boldsymbol{\phi}}_{i})\\\
&-\left(\alpha_{P}-\frac{1}{2}\right)b({\boldsymbol{\phi}}_{i},{\boldsymbol{m}},{\boldsymbol{u}})-\left(\alpha_{P}+\frac{1}{2}\right)b({\boldsymbol{\phi}}_{i},{\boldsymbol{u}},{\boldsymbol{m}})-\frac{1}{2}\left({\boldsymbol{\phi}}_{i}{\cdot}\nabla\rho,{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\boldsymbol{u}}{\cdot}{\boldsymbol{u}}}\right),\end{split}$
(3.25)

since $\nabla{\cdot\,}{\boldsymbol{\phi}}_{i}=0$. Using (2.5) on
$\left(\alpha_{P}-\frac{1}{2}\right)b({\boldsymbol{\phi}}_{i},{\boldsymbol{m}},{\boldsymbol{u}})$
and $b({\boldsymbol{u}},{\boldsymbol{m}},{\boldsymbol{\phi}}_{i})$ and
integration by parts on the pressure term yield

$\begin{split}{\partial}_{t}\int_{\Omega}({\boldsymbol{m}}{\times}{\boldsymbol{x}})_{i}\,\mathrm{d}{\boldsymbol{x}}=&\
b({\boldsymbol{u}},{\boldsymbol{\phi}}_{i},{\boldsymbol{m}})+((\nabla{\cdot\,}{\boldsymbol{u}}){\boldsymbol{m}},{\boldsymbol{\phi}}_{i})-\alpha_{\boldsymbol{m}}((\nabla{\cdot\,}{\boldsymbol{u}}){\boldsymbol{m}},{\boldsymbol{\phi}}_{i})+(P,\nabla{\cdot\,}{\boldsymbol{\phi}}_{i})\\\
&+\left(\alpha_{P}-\frac{1}{2}\right)b({\boldsymbol{\phi}}_{i},{\boldsymbol{u}},{\boldsymbol{m}})-\left(\alpha_{P}+\frac{1}{2}\right)b({\boldsymbol{\phi}}_{i},{\boldsymbol{u}},{\boldsymbol{m}})-\frac{1}{2}\left({\boldsymbol{\phi}}_{i}{\cdot}\nabla\rho,{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\boldsymbol{u}}{\cdot}{\boldsymbol{u}}}\right)\\\
=&\
b({\boldsymbol{u}},{\boldsymbol{\phi}}_{i},{\boldsymbol{m}})+(1-\alpha_{\boldsymbol{m}})((\nabla{\cdot\,}{\boldsymbol{u}}){\boldsymbol{m}},{\boldsymbol{\phi}}_{i})-b({\boldsymbol{\phi}}_{i},{\boldsymbol{u}},{\boldsymbol{m}})-\frac{1}{2}\left({\boldsymbol{\phi}}_{i}{\cdot}\nabla\rho,{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\boldsymbol{u}}{\cdot}{\boldsymbol{u}}}\right)\\\
=&\
b({\boldsymbol{u}},{\boldsymbol{\phi}}_{i},{\boldsymbol{m}})+(1-\alpha_{\boldsymbol{m}})((\nabla{\cdot\,}{\boldsymbol{u}}){\boldsymbol{m}},{\boldsymbol{\phi}}_{i})-b(\rho{\boldsymbol{\phi}}_{i},{\boldsymbol{u}},{\boldsymbol{u}})-\frac{1}{2}\left({\boldsymbol{\phi}}_{i}{\cdot}\nabla\rho,{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\boldsymbol{u}}{\cdot}{\boldsymbol{u}}}\right)\\\
=&\
b({\boldsymbol{u}},{\boldsymbol{\phi}}_{i},{\boldsymbol{m}})+(1-\alpha_{\boldsymbol{m}})((\nabla{\cdot\,}{\boldsymbol{u}}){\boldsymbol{m}},{\boldsymbol{\phi}}_{i})+\frac{1}{2}((\nabla{\cdot\,}(\rho{\boldsymbol{\phi}}_{i})){\boldsymbol{u}},{\boldsymbol{u}})-\frac{1}{2}\left({\boldsymbol{\phi}}_{i}{\cdot}\nabla\rho,{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\boldsymbol{u}}{\cdot}{\boldsymbol{u}}}\right)\\\
=&\
b({\boldsymbol{u}},{\boldsymbol{\phi}}_{i},{\boldsymbol{m}})+(1-\alpha_{\boldsymbol{m}})((\nabla{\cdot\,}{\boldsymbol{u}}){\boldsymbol{m}},{\boldsymbol{\phi}}_{i}),\end{split}$
(3.26)

where (2.6) was used in the third equality and (2.4) was used in the second
equality. Next, by using (2.2) it can be shown that

$\begin{split}b({\boldsymbol{u}},{\boldsymbol{\phi}}_{i},{\boldsymbol{m}})&=b({\boldsymbol{u}},{\boldsymbol{\phi}}_{i},{\boldsymbol{u}}\rho)=\frac{1}{2}b({\boldsymbol{u}},{\boldsymbol{\phi}}_{i},{\boldsymbol{u}}\rho)+\frac{1}{2}b({\boldsymbol{u}},{\boldsymbol{\phi}}_{i},{\boldsymbol{u}}\rho)\\\
&=\frac{1}{2}\left((\nabla{\boldsymbol{\phi}}_{i})^{\top}{\boldsymbol{u}},{\boldsymbol{u}}\rho\right)+\frac{1}{2}\left((\nabla{\boldsymbol{\phi}}_{i})^{\top}{\boldsymbol{u}},{\boldsymbol{u}}\rho\right)=\frac{1}{2}\left((\nabla{\boldsymbol{\phi}}_{i})^{\top}{\boldsymbol{u}},{\boldsymbol{u}}\rho\right)+\frac{1}{2}((\nabla{\boldsymbol{\phi}}_{i}){\boldsymbol{u}},{\boldsymbol{u}}\rho)=0,\end{split}$
(3.27)

since
$\nabla{\boldsymbol{\phi}}_{i}+(\nabla{\boldsymbol{\phi}}_{i})^{\top}=0$.
Since $b({\boldsymbol{u}},{\boldsymbol{\phi}}_{i},{\boldsymbol{m}})=0$,
angular momentum is conserved if $\alpha_{\boldsymbol{m}}=1$.

∎

##  4 Properties of viscous regularizations.

In this section, we investigate how viscous regularizations of the
incompressible Euler equations (2.1) affect conservation properties. We note
that the purpose of this section is to investigate the properties of the model
and not to investigate the properties of the numerical method. To simplify the
analysis we assume that $\nabla{\cdot\,}{\boldsymbol{u}}=0$ holds pointwise
and assume periodic boundary conditions. We consider a viscous regularization
of the following form

$\displaystyle{\partial}_{t}\rho+{\boldsymbol{u}}{\cdot}\nabla\rho$
$\displaystyle=\nabla{\cdot\,}{\boldsymbol{f}}_{\rho},$ (4.1)
$\displaystyle{\partial}_{t}{\boldsymbol{m}}+{\boldsymbol{u}}{\cdot}\nabla{\boldsymbol{m}}+\nabla
p$
$\displaystyle={\boldsymbol{f}}+\nabla{\cdot\,}\left(\mu\left(\nabla{\boldsymbol{u}}+(\nabla{\boldsymbol{u}})^{\top}\right)\right)+{\boldsymbol{f}}_{{\boldsymbol{m}}},\quad$
$\displaystyle({\boldsymbol{x}},t)\in\Omega{\times}(0,T],$
$\displaystyle\nabla{\cdot\,}{\boldsymbol{u}}$ $\displaystyle=0,$
$\displaystyle{\boldsymbol{u}}({\boldsymbol{x}},0)$
$\displaystyle={\boldsymbol{u}}_{0}({\boldsymbol{x}}),$
$\displaystyle\rho({\boldsymbol{x}},0)$
$\displaystyle=\rho_{0}({\boldsymbol{x}}),\quad$
$\displaystyle{\boldsymbol{x}}\in\Omega,$

where $\mu\geq 0$ is the dynamic viscosity coefficient which in general is
space-dependent, ${\boldsymbol{f}}_{\rho}\mathrel{\mathop{\mathchar
58\relax}}=\kappa\nabla\rho$, where $\kappa\geq 0$ is a mass diffusivity
coefficient which also is space-dependent. In this work, we define
${\boldsymbol{f}}_{\boldsymbol{m}}$ as

${\boldsymbol{f}}_{{\boldsymbol{m}}}({\boldsymbol{u}},{\boldsymbol{f}}_{\rho})\mathrel{\mathop{\mathchar
58\relax}}=A_{1}(\nabla{\boldsymbol{f}}_{\rho}){\boldsymbol{u}}+A_{2}(\nabla{\boldsymbol{f}}_{\rho})^{\top}{\boldsymbol{u}}+A_{3}(\nabla{\boldsymbol{u}}){\boldsymbol{f}}_{\rho}+A_{4}(\nabla{\boldsymbol{u}})^{\top}{\boldsymbol{f}}_{\rho}+A_{5}(\nabla{\cdot\,}{\boldsymbol{f}}_{\rho}){\boldsymbol{u}},$
(4.2)

where $A_{i}\in{\mathbb{R}}$. Depending on how the constants $A_{i}$ are
chosen, different properties of the model are obtained and this is summarized
in Theorem 4.1.

######  Theorem 4.1.

If ${\boldsymbol{f}}=0$, the system (4.1):

  1. 1.

Dissipates kinetic energy if $A_{1}+A_{2}-A_{3}=0\text{ and
}A_{5}-\frac{1}{2}A_{4}=\frac{1}{2}$.

  2. 2.

Conserves momentum if $A_{1}-A_{3}=0\text{ and }A_{5}-A_{4}=0$.

  3. 3.

Conserves angular momentum if $A_{1}-A_{3}=0\text{ and }A_{5}-A_{4}=0\text{
and }A_{2}=A_{4}$.

###### Proof.

We divide the proof into several sections. Extending the proofs to other
boundary conditions such as no-slip (${\boldsymbol{u}}=0|_{{\partial}\Omega}$)
can be done by using the proof technique in [5, Sec 3.1.2] or [16, Sec 3.1.2].

Kinetic energy: Repeating similar steps on (4.1) as in the proof of Theorem
3.1 and assuming $\nabla{\cdot\,}{\boldsymbol{u}}=0$ yields

$\begin{split}\frac{1}{2}{\partial}_{t}\int_{\Omega}\rho{\boldsymbol{u}}{\cdot}{\boldsymbol{u}}\,\mathrm{d}{\boldsymbol{x}}=&\
({\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\partial}_{t}}{\boldsymbol{m}},{\boldsymbol{u}})-\frac{1}{2}\left({\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\partial}_{t}}\rho,{\boldsymbol{u}}{\cdot}{\boldsymbol{u}}\right)=\left({\boldsymbol{f}}_{{\boldsymbol{m}}},{\boldsymbol{u}}\right)-\frac{1}{2}\left(\left(\nabla{\cdot\,}{\boldsymbol{f}}_{\rho}\right){\boldsymbol{u}},{\boldsymbol{u}}\right)+\left(\nabla{\cdot\,}\left(\mu\left(\nabla{\boldsymbol{u}}+(\nabla{\boldsymbol{u}})^{\top}\right)\right),{\boldsymbol{u}}\right),\\\
\end{split}$ (4.3)

where only the viscous terms determine if kinetic energy is dissipated. Next,
we obtain

$\begin{split}\frac{1}{2}{\partial}_{t}\int_{\Omega}\rho{\boldsymbol{u}}{\cdot}{\boldsymbol{u}}\,\mathrm{d}{\boldsymbol{x}}=&\
A_{1}b({\boldsymbol{u}},{\boldsymbol{f}}_{\rho},{\boldsymbol{u}})+A_{2}b({\boldsymbol{u}},{\boldsymbol{f}}_{\rho},{\boldsymbol{u}})+A_{3}b({\boldsymbol{u}},{\boldsymbol{u}},{\boldsymbol{f}}_{\rho})+A_{4}b({\boldsymbol{f}}_{\rho},{\boldsymbol{u}},{\boldsymbol{u}})\\\
&+A_{5}((\nabla{\cdot\,}{\boldsymbol{f}}_{\rho}){\boldsymbol{u}},{\boldsymbol{u}})-\frac{1}{2}((\nabla{\cdot\,}{\boldsymbol{f}}_{\rho}){\boldsymbol{u}},{\boldsymbol{u}})+\left(\nabla{\cdot\,}\left(\mu\left(\nabla{\boldsymbol{u}}+(\nabla{\boldsymbol{u}})^{\top}\right)\right),{\boldsymbol{u}}\right),\\\
=&\
A_{1}b({\boldsymbol{u}},{\boldsymbol{f}}_{\rho},{\boldsymbol{u}})+A_{2}b({\boldsymbol{u}},{\boldsymbol{f}}_{\rho},{\boldsymbol{u}})-A_{3}b({\boldsymbol{u}},{\boldsymbol{f}}_{\rho},{\boldsymbol{u}})-\frac{1}{2}A_{4}((\nabla{\cdot\,}{\boldsymbol{f}}_{\rho}){\boldsymbol{u}},{\boldsymbol{u}})\\\
&+A_{5}((\nabla{\cdot\,}{\boldsymbol{f}}_{\rho}){\boldsymbol{u}},{\boldsymbol{u}})-\frac{1}{2}((\nabla{\cdot\,}{\boldsymbol{f}}_{\rho}){\boldsymbol{u}},{\boldsymbol{u}})-\left(\mu\left(\nabla{\boldsymbol{u}}+(\nabla{\boldsymbol{u}})^{\top}\right),\nabla{\boldsymbol{u}}\right),\end{split}$
(4.4)

by using integration by parts on the $\mu$ term and using (2.5) on the $A_{3}$
term and (2.6) on the $A_{4}$ term. Using that the contraction between a
symmetric and anti-symmetric matrix is zero [18, Ch 11.2.1], i.e.,
$(A+A^{\top})\mathrel{\mathop{\mathchar 58\relax}}(A-A^{\top})=0$ where $A$ is
a square matrix, one can show that

$\begin{split}\left(\mu\left(\nabla{\boldsymbol{u}}+\left(\nabla{\boldsymbol{u}}\right)^{\top}\right),\nabla{\boldsymbol{u}}\right)=&\frac{1}{2}\left(\mu\left(\nabla{\boldsymbol{u}}+\left(\nabla{\boldsymbol{u}}\right)^{\top}\right),\nabla{\boldsymbol{u}}+(\nabla{\boldsymbol{u}})^{\top}\right)\\\
&+\frac{1}{2}\left(\mu\left(\nabla{\boldsymbol{u}}+\left(\nabla{\boldsymbol{u}}\right)^{\top}\right),\nabla{\boldsymbol{u}}-(\nabla{\boldsymbol{u}})^{\top}\right)\\\
=&\frac{1}{2}\left(\mu\left(\nabla{\boldsymbol{u}}+\left(\nabla{\boldsymbol{u}}\right)^{\top}\right),\nabla{\boldsymbol{u}}+(\nabla{\boldsymbol{u}})^{\top}\right).\end{split}$
(4.5)

Inserting (4.5) into (4.4) gives

$\begin{split}\frac{1}{2}{\partial}_{t}\int_{\Omega}\rho{\boldsymbol{u}}{\cdot}{\boldsymbol{u}}\,\mathrm{d}{\boldsymbol{x}}=&\
A_{1}b({\boldsymbol{u}},{\boldsymbol{f}}_{\rho},{\boldsymbol{u}})+A_{2}b({\boldsymbol{u}},{\boldsymbol{f}}_{\rho},{\boldsymbol{u}})-A_{3}b({\boldsymbol{u}},{\boldsymbol{f}}_{\rho},{\boldsymbol{u}})-\frac{1}{2}A_{4}((\nabla{\cdot\,}{\boldsymbol{f}}_{\rho}){\boldsymbol{u}},{\boldsymbol{u}})\\\
&+A_{5}((\nabla{\cdot\,}{\boldsymbol{f}}_{\rho}){\boldsymbol{u}},{\boldsymbol{u}})-\frac{1}{2}((\nabla{\cdot\,}{\boldsymbol{f}}_{\rho}){\boldsymbol{u}},{\boldsymbol{u}})-\frac{1}{2}\|\sqrt{\mu}\left(\nabla{\boldsymbol{u}}+(\nabla{\boldsymbol{u}})^{\top}\right)\|^{2}.\end{split}$
(4.6)

From (4.6), the conditions necessary for kinetic energy to be dissipated are

$A_{1}+A_{2}-A_{3}=0\text{ and }A_{5}-\frac{1}{2}A_{4}=\frac{1}{2}.$ (4.7)

Momentum: Repeating similar steps on (4.1) as in the proof of Theorem 3.1 and
assuming $\nabla{\cdot\,}{\boldsymbol{u}}=0$ yields

${\partial}_{t}\int_{\Omega}{\boldsymbol{m}}_{i}\,\mathrm{d}{\boldsymbol{x}}=({\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\partial}_{t}}{\boldsymbol{m}},{\boldsymbol{e}}_{i})=\left(\nabla{\cdot\,}\left(\mu\left(\nabla{\boldsymbol{u}}+(\nabla{\boldsymbol{u}})^{\top}\right)\right),{\boldsymbol{e}}_{i}\right)+({\boldsymbol{f}}_{{\boldsymbol{m}}},{\boldsymbol{e}}_{i}),$
(4.8)

where only the viscous terms determine if momentum is conserved. Performing
integration by parts on the $\mu$ term yields

${\partial}_{t}\int_{\Omega}{\boldsymbol{m}}_{i}\,\mathrm{d}{\boldsymbol{x}}=A_{1}b({\boldsymbol{e}}_{i},{\boldsymbol{f}}_{\rho},{\boldsymbol{u}})+A_{2}b({\boldsymbol{u}},{\boldsymbol{f}}_{\rho},{\boldsymbol{e}}_{i})+A_{3}b({\boldsymbol{e}}_{i},{\boldsymbol{u}},{\boldsymbol{f}}_{\rho})+A_{4}b({\boldsymbol{f}}_{\rho},{\boldsymbol{u}},{\boldsymbol{e}}_{i})+A_{5}((\nabla{\cdot\,}{\boldsymbol{f}}_{\rho}){\boldsymbol{u}},{\boldsymbol{e}}_{i}).$
(4.9)

Using (2.5) on the $A_{3}$, $A_{2}$ and $A_{4}$ terms inside (4.9) gives

${\partial}_{t}\int_{\Omega}{\boldsymbol{m}}_{i}\,\mathrm{d}{\boldsymbol{x}}=A_{1}b({\boldsymbol{e}}_{i},{\boldsymbol{f}}_{\rho},{\boldsymbol{u}})-A_{3}b({\boldsymbol{e}}_{i},{\boldsymbol{f}}_{\rho},{\boldsymbol{u}})-A_{4}((\nabla{\cdot\,}{\boldsymbol{f}}_{\rho}){\boldsymbol{u}},{\boldsymbol{e}}_{i})+A_{5}((\nabla{\cdot\,}{\boldsymbol{f}}_{\rho}){\boldsymbol{u}},{\boldsymbol{e}}_{i}),$
(4.10)

which leads to the conditions necessary for momentum to be conserved

$A_{1}-A_{3}=0\text{ and }A_{5}-A_{4}=0.$ (4.11)

Angular momentum: Repeating similar steps on (4.1) as in the proof of Theorem
3.1 and assuming $\nabla{\cdot\,}{\boldsymbol{u}}=0$ yields

${\partial}_{t}\int_{\Omega}({\boldsymbol{m}}{\times}{\boldsymbol{x}})_{i}\,\mathrm{d}{\boldsymbol{x}}=\left(\nabla{\cdot\,}\left(\mu\left(\nabla{\boldsymbol{u}}+(\nabla{\boldsymbol{u}})^{\top}\right)\right),{\boldsymbol{\phi}}_{i}\right)+({\boldsymbol{f}}_{{\boldsymbol{m}}},{\boldsymbol{\phi}}_{i}).$
(4.12)

Performing integration by parts on the $\mu$ term yields

${\partial}_{t}\int_{\Omega}({\boldsymbol{m}}{\times}{\boldsymbol{x}})_{i}\,\mathrm{d}{\boldsymbol{x}}=-(\mu\left(\nabla{\boldsymbol{u}}+(\nabla{\boldsymbol{u}})^{\top}\right),\nabla{\boldsymbol{\phi}}_{i})+({\boldsymbol{f}}_{{\boldsymbol{m}}},{\boldsymbol{\phi}}_{i}).$
(4.13)

Next, using that the contraction between a symmetric and anti-symmetric matrix
is zero [18, Ch 11.2.1], one can show that

$\begin{split}\left(\mu\left(\nabla{\boldsymbol{u}}+(\nabla{\boldsymbol{u}})^{\top}\right),\nabla{\boldsymbol{\phi}}_{i}\right)=&\left(\frac{1}{2}\mu\left(\nabla{\boldsymbol{u}}+\left(\nabla{\boldsymbol{u}}\right)^{\top}\right),\nabla{\boldsymbol{\phi}}_{i}+(\nabla{\boldsymbol{\phi}}_{i})^{\top}\right)\\\
&+\left(\frac{1}{2}\mu\left(\nabla{\boldsymbol{u}}+\left(\nabla{\boldsymbol{u}}\right)^{\top}\right),\nabla{\boldsymbol{\phi}}_{i}-(\nabla{\boldsymbol{\phi}}_{i})^{\top}\right)\\\
=&\left(\frac{1}{2}\mu\left(\nabla{\boldsymbol{u}}+\left(\nabla{\boldsymbol{u}}\right)^{\top}\right),\nabla{\boldsymbol{\phi}}_{i}+(\nabla{\boldsymbol{\phi}}_{i})^{\top}\right).\end{split}$
(4.14)

Since $\nabla{\boldsymbol{\phi}}_{i}=-\nabla{\boldsymbol{\phi}}_{i}^{\top}$ we
are left with

$\begin{split}{\partial}_{t}\int_{\Omega}({\boldsymbol{m}}{\times}{\boldsymbol{x}})_{i}\,\mathrm{d}{\boldsymbol{x}}=&\
({\boldsymbol{f}}_{{\boldsymbol{m}}},{\boldsymbol{\phi}}_{i})=A_{1}b({\boldsymbol{\phi}}_{i},{\boldsymbol{f}}_{\rho},{\boldsymbol{u}})+A_{2}b({\boldsymbol{u}},{\boldsymbol{f}}_{\rho},{\boldsymbol{\phi}}_{i})\\\
&+A_{3}b({\boldsymbol{\phi}}_{i},{\boldsymbol{u}},{\boldsymbol{f}}_{\rho})+A_{4}b({\boldsymbol{f}}_{\rho},{\boldsymbol{u}},{\boldsymbol{\phi}}_{i})+A_{5}((\nabla{\cdot\,}{\boldsymbol{f}}_{\rho}){\boldsymbol{u}},{\boldsymbol{\phi}}_{i}).\end{split}$
(4.15)

Using (2.5) on the $A_{3}$, $A_{2}$ and $A_{4}$ terms and using that
$\nabla{\cdot\,}{\boldsymbol{u}}=\nabla{\cdot\,}{\boldsymbol{\phi}}_{i}=0$
gives

$\begin{split}{\partial}_{t}\int_{\Omega}({\boldsymbol{m}}{\times}{\boldsymbol{x}})_{i}\,\mathrm{d}{\boldsymbol{x}}=&\
A_{1}b({\boldsymbol{\phi}}_{i},{\boldsymbol{f}}_{\rho},{\boldsymbol{u}})-A_{2}b({\boldsymbol{u}},{\boldsymbol{\phi}}_{i},{\boldsymbol{f}}_{\rho})-A_{3}b({\boldsymbol{\phi}}_{i},{\boldsymbol{f}}_{\rho},{\boldsymbol{u}})\\\
&-A_{4}b({\boldsymbol{f}}_{\rho},{\boldsymbol{\phi}}_{i},{\boldsymbol{u}})-A_{4}((\nabla{\cdot\,}{\boldsymbol{f}}_{\rho}){\boldsymbol{u}},{\boldsymbol{\phi}}_{i})+A_{5}((\nabla{\cdot\,}{\boldsymbol{f}}_{\rho}){\boldsymbol{u}},{\boldsymbol{\phi}}_{i}).\end{split}$
(4.16)

Next, using the definition of the trilinear form (2.2) yields

$\begin{split}{\partial}_{t}\int_{\Omega}({\boldsymbol{m}}{\times}{\boldsymbol{x}})_{i}\,\mathrm{d}{\boldsymbol{x}}=&\
A_{1}b({\boldsymbol{\phi}}_{i},{\boldsymbol{f}}_{\rho},{\boldsymbol{u}})-A_{2}\left((\nabla{\boldsymbol{\phi}}_{i})^{\top}{\boldsymbol{u}},{\boldsymbol{f}}_{\rho}\right)-A_{3}b({\boldsymbol{\phi}}_{i},{\boldsymbol{f}}_{\rho},{\boldsymbol{u}})\\\
&-A_{4}((\nabla{\boldsymbol{\phi}}_{i}){\boldsymbol{u}},{\boldsymbol{f}}_{\rho})-A_{4}((\nabla{\cdot\,}{\boldsymbol{f}}_{\rho}){\boldsymbol{u}},{\boldsymbol{\phi}}_{i})+A_{5}((\nabla{\cdot\,}{\boldsymbol{f}}_{\rho}){\boldsymbol{u}},{\boldsymbol{\phi}}_{i}),\end{split}$
(4.17)

which leads to the conditions necessary for angular momentum to be conserved

$A_{1}-A_{3}=0\text{ and }A_{5}-A_{4}=0\text{ and }A_{2}=A_{4},$ (4.18)

since
$\nabla{\boldsymbol{\phi}}_{i}+(\nabla{\boldsymbol{\phi}}_{i})^{\top}=0$.

∎

###  4.1 Discussion.

One implication of Theorem 4.1 is that there is no viscous regularization
$(\kappa>0)$ that simultaneously dissipates kinetic energy and conserves
angular momentum. Another implication is that a unique viscous regularization
dissipates kinetic energy and conserves momentum. This can be seen by
combining the kinetic energy and the momentum conditions in Theorem 4.1 to
arrive at $A_{5}=A_{4}=1$, $A_{2}=0$ and $A_{1}=A_{3}$. Setting these
conditions inside (4.2) and performing some simplifications lead to the
following viscous flux

${\boldsymbol{f}}_{GP}\mathrel{\mathop{\mathchar
58\relax}}=(\nabla{\boldsymbol{u}})^{\top}{\boldsymbol{f}}_{\rho}+(\nabla{\cdot\,}{\boldsymbol{f}}_{\rho}){\boldsymbol{u}}+A_{1}(\nabla{\boldsymbol{f}}_{\rho}){\boldsymbol{u}}+A_{1}(\nabla{\boldsymbol{u}}){\boldsymbol{f}}_{\rho}=\nabla{\cdot\,}({\boldsymbol{f}}_{\rho}\otimes{\boldsymbol{u}})+A_{1}\nabla({\boldsymbol{f}}_{\rho}{\cdot}{\boldsymbol{u}}),$
(4.19)

where the second term does not change the model since it only results that a
modified pressure is solved for, i.e.,
$P=p-A_{1}{\boldsymbol{f}}_{\rho}{\cdot}{\boldsymbol{u}}$, and does not affect
the behavior of density and velocity. Since that term is not important, we are
left with the so-called Guermond-Popov viscous flux [12]. Compared to (4.19),
the Guermond-Popov flux [12] is simplified since it now is in the
incompressible setting instead of the compressible setting. This means that,
if we only consider the density and momentum contribution of [12], it is
exactly the same as (4.19) and (4.1).

Combining the last two conditions of Theorem 4.1 leads to a family of viscous
regularizations that conserve momentum and angular momentum. Below two
examples of these are given

$\displaystyle{\boldsymbol{f}}_{0}\mathrel{\mathop{\mathchar 58\relax}}=0,$
(4.20)
$\displaystyle{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\boldsymbol{f}}_{KS}\mathrel{\mathop{\mathchar
58\relax}}=\nabla{\cdot\,}({\boldsymbol{f}}_{\rho}\otimes{\boldsymbol{u}}+{\boldsymbol{u}}\otimes{\boldsymbol{f}}_{\rho}),}$
(4.21)

where ${\boldsymbol{f}}_{0}$ has previously been used in the variable density
flow literature [3]. From Ficks law, Kažihov and Smagulov [17] derived and
proposed (4.21). They showed kinetic energy dissipation assuming that $\kappa$
is constant and

$\kappa<\frac{2\mu}{\max_{\Omega}\rho-\min_{\Omega}\rho}.$ (4.22)

It is not clear if (4.22) is true for all fluid pairs at all temperatures. For
example, consider hydrogen gas $\text{H}_{2}$ and sulfur hexafluoride gas
$\text{SF}_{6}$ at room temperature. Values of dynamic viscosity, density and
pairwise mass diffusivity are given in Table 3. Inserting these values into
(4.22) shows an order of magnitude violation:

$\begin{split}\frac{2\mu}{\max_{\Omega}\rho-\min_{\Omega}\rho}&=\frac{2\cdot
15.3\cdot 10^{-6}}{6.07-0.0838}\text{ [m${}^{2}$/sec]}\\\ &=5.1\cdot
10^{-6}\text{ [m${}^{2}$/sec],}\\\ \\\ \kappa=4.12\cdot 10^{-5}\text{
[m${}^{2}$/sec]}&\nless 5.1\cdot 10^{-6}\text{ [m${}^{2}$/sec].}\end{split}$
(4.23)

Table 3: Experimental values from [15] of dynamic viscosity at 300 K, 1 kPa and pairwise mass diffusion at 20 ℃, 1 atm. Density from ideal gas law. |  $\mu$ [Pa $\cdot$ sec] |  $\rho$ [kg/m3] |  $\kappa$ [m2/sec]
---|---|---|---
$\text{H}_{2}$ | $8.9\cdot 10^{-6}$ | 0.0838 | $4.12\cdot 10^{-5}$
$\text{SF}_{6}$ | $15.3\cdot 10^{-6}$ | 6.07 | $4.12\cdot 10^{-5}$

##  5 Fully discrete conservative method.

In this section, we describe the fully discrete method approximating the
incompressible Navier-Stokes equations (4.1). Let
$(\rho^{n},{\boldsymbol{u}}^{n},P^{n})\in{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\mathcal{M}}{\times}\mathbf{{\mathbfcal
V}}{\times}{\mathcal{Q}}}$ be solutions at time $t_{n}$, where the time-levels
are defined as $0=t_{0}<t_{1}<\dots<t_{N_{t}}=T$ and ${N_{t}}$ denotes the
total number of time-levels. Let $\Delta t_{n+j}\mathrel{\mathop{\mathchar
58\relax}}=t_{n+j+1}-t_{n+j}$ denote the time step. We also define
${\boldsymbol{u}}^{n+1/2}\mathrel{\mathop{\mathchar
58\relax}}=\frac{{\boldsymbol{u}}^{n}+{\boldsymbol{u}}^{n+1}}{2}$.

We propose using the second-order accurate time-stepping method based on the
modified Crank-Nicolson method by Gawlik and Gay-Balmaz [8] which has
favorable fully discrete properties. Applying this method on (3.5) yields:
Given $\left(\rho^{n},{\boldsymbol{u}}^{n}\right)$, find
$\left(\rho^{n+1},{\boldsymbol{u}}^{n+1},P^{n+1/2}\right)\in{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\mathcal{M}}{\times}\mathbf{{\mathbfcal
V}}{\times}{\mathcal{Q}}}$ such that

$\displaystyle\left(\frac{\rho^{n+1}-\rho^{n}}{\Delta
t_{n}}+{\boldsymbol{u}}^{n+1/2}{\cdot}\nabla\rho^{n+1/2}+\alpha_{\rho}\left(\nabla{\cdot\,}{\boldsymbol{u}}^{n+1/2}\right)\left(\rho^{n+1/2}-\bar{\rho}^{n+1/2}\right),w\right)$
(5.1) $\displaystyle+\left(\kappa\nabla\rho^{n+1/2},\nabla w\right)$
$\displaystyle=0,\quad$ $\displaystyle\forall w\in{\mathcal{M}},$
$\displaystyle\left(\frac{{\boldsymbol{m}}^{n+1}-{\boldsymbol{m}}^{n}}{\Delta
t_{n}},{\boldsymbol{v}}\right)+b\left({\boldsymbol{u}}^{n+1/2},{\boldsymbol{m}}^{n+1/2},{\boldsymbol{v}}\right)$
$\displaystyle+\alpha_{\boldsymbol{m}}\left(\left(\nabla{\cdot\,}{\boldsymbol{u}}^{n+1/2}\right){\boldsymbol{m}}^{n+1/2},{\boldsymbol{v}}\right)+\left(\nabla
P^{n+1/2},{\boldsymbol{v}}\right)$
$\displaystyle+\left(\alpha_{P}-\frac{1}{2}\right)b\left({\boldsymbol{v}},{\boldsymbol{m}}^{n+1/2},{\boldsymbol{u}}^{n+1/2}\right)$
$\displaystyle+\left(\alpha_{P}+\frac{1}{2}\right)b\left({\boldsymbol{v}},{\boldsymbol{u}}^{n+1/2},{\boldsymbol{m}}^{n+1/2}\right)+\frac{\alpha_{\rho}}{2}\left(\nabla\left({\boldsymbol{m}}^{n+1/2}{\cdot}{\boldsymbol{u}}^{n+1/2}\right),{\boldsymbol{v}}\right)$
$\displaystyle+\frac{1}{2}\left(\left({\boldsymbol{v}}{\cdot}\nabla\rho^{n+1/2},{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\boldsymbol{u}}^{n}{\cdot}{\boldsymbol{u}}^{n+1}}\right)-\alpha_{\rho}\left(\nabla\left(\left(\rho^{n+1/2}-\bar{\rho}^{n+1/2}\right){\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\boldsymbol{u}}^{n}{\cdot}{\boldsymbol{u}}^{n+1}}\right),{\boldsymbol{v}}\right)\right)$
$\displaystyle+\left(\mu\left(\nabla{\boldsymbol{u}}^{n+1/2}+\left(\nabla{\boldsymbol{u}}^{n+1/2}\right)^{\top}\right),\nabla{\boldsymbol{v}}\right)$
$\displaystyle-\left({\boldsymbol{f}}^{n+1/2},{\boldsymbol{v}}\right)-\left({\boldsymbol{f}}_{{\boldsymbol{m}}}\left({\boldsymbol{u}}^{n+1/2},{\boldsymbol{f}}_{\rho}^{n+1/2}\right),{\boldsymbol{v}}\right)$
$\displaystyle=0,\quad$ $\displaystyle\forall{\boldsymbol{v}}\in{\mathbfcal
V},$ $\displaystyle\left(\nabla{\cdot\,}{\boldsymbol{u}}^{n+1/2},q\right)$
$\displaystyle=0,\quad$ $\displaystyle\forall q\in{\mathcal{Q}}.$

Note that the term
${\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\boldsymbol{u}}{\cdot}{\boldsymbol{u}}}$
has been discretized as
${\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\boldsymbol{u}}^{n}{\cdot}{\boldsymbol{u}}^{n+1}}$
which is different from the usual Crank-Nicolson method.

####  5.0.1 Properties of the fully discrete method.

This section is a fully discrete analog of Section 3.2 and aims to verify that
the fully discrete method satisfies the properties of Table 1. Provided that
the conditions of Table 1 are fulfilled, all of the properties are retained at
the fully discrete level except conservation of momentum and angular momentum.
To conserve momentum and angular momentum, the term
${\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\boldsymbol{u}}^{n}{\cdot}{\boldsymbol{u}}^{n+1}}$
should be discretized as
${\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\boldsymbol{u}}^{n+1/2}{\cdot}{\boldsymbol{u}}^{n+1/2}}$
instead but this will lead to a loss of conservation of kinetic energy.
Additional techniques are required to conserve both kinetic energy and
momentum and angular momentum. For example, recently Zhang et al. [30] used
the so-called scalar auxiliary variable (SAV) technique to enhance their time-
stepping scheme to be able to conserve both kinetic energy and momentum.

######  Proposition 5.1.

The properties of Table 1, except momentum and angular momentum, hold at the
fully discrete level for (5.1) provided that the conditions inside the table
are fulfilled.

Mass: Repeating the steps of the proof of Theorem 3.1 on (5.1) yields

$\int_{\Omega}\rho^{n+1}\,\mathrm{d}{\boldsymbol{x}}=\int_{\Omega}\rho^{n}\,\mathrm{d}{\boldsymbol{x}}+\Delta
t_{n}\left(1-\alpha_{\rho}\right)\left(\nabla{\cdot\,}{\boldsymbol{u}}^{n+1/2},\rho^{n+1/2}-\bar{\rho}^{n+1/2}\right),$
(5.2)

and similar steps can be performed to show that
$\left(\rho^{n+1/2}-\bar{\rho}^{n+1/2}\right)\in{\mathcal{Q}}$ provided that
$k_{\rho}\leq k_{P}$. Mass is conserved if $k_{\rho}\leq k_{P}$ or if
$\alpha_{\rho}=1$.

Squared density: Setting $w=\rho^{n+1/2}-\bar{\rho}^{n+1/2}$ inside (5.1) and
repeating the steps of the proof of Theorem 3.1 yields

$\left(\rho^{n+1}-\rho^{n},\rho^{n+1/2}-\bar{\rho}^{n+1/2}\right)=\Delta
t_{n}\left(\frac{1}{2}-\alpha_{\rho}\right)\left(\left(\nabla{\cdot\,}{\boldsymbol{u}}^{n+1/2}\right)\left(\rho^{n+1/2}-\bar{\rho}^{n+1/2}\right),\rho^{n+1/2}-\bar{\rho}^{n+1/2}\right).$
(5.3)

Simplifying the left-hand side finally yields

$\begin{split}&\frac{1}{2}\|\rho^{n+1}\|^{2}-\left(\rho^{n+1},\bar{\rho}^{n+1/2}\right)=\frac{1}{2}\|\rho^{n}\|^{2}-\left(\rho^{n},\bar{\rho}^{n+1/2}\right)\\\
&+\Delta
t_{n}\left(\frac{1}{2}-\alpha_{\rho}\right)\left(\left(\nabla{\cdot\,}{\boldsymbol{u}}^{n+1/2}\right)\left(\rho^{n+1/2}-\bar{\rho}^{n+1/2}\right),\rho^{n+1/2}-\bar{\rho}^{n+1/2}\right),\end{split}$
(5.4)

which is a fully discrete equivalent of (3.12). If
$\bar{\rho}^{n+1/2}=\frac{1}{|\Omega|}\int_{\Omega}\rho^{n+1/2}\,\mathrm{d}{\boldsymbol{x}}$
then

$\frac{1}{2}\|\rho^{n+1}\|^{2}-\left(\rho^{n+1},\bar{\rho}^{n+1/2}\right)-\left(\frac{1}{2}\|\rho^{n}\|^{2}-\left(\rho^{n},\bar{\rho}^{n+1/2}\right)\right),$
(5.5)

simplifies to

$\frac{1}{2}\left(\|\rho^{n+1}\|^{2}-\left(\int_{\Omega}\rho^{n+1}\,\mathrm{d}{\boldsymbol{x}}\right)^{2}\right)-\frac{1}{2}\left(\|\rho^{n}\|^{2}-\left(\int_{\Omega}\rho^{n}\,\mathrm{d}{\boldsymbol{x}}\right)^{2}\right),$
(5.6)

which is a fully discrete equivalent of (3.13).

Lastly, similar steps can be performed to show that for $A\in{\mathbb{R}}$

$\begin{split}&\frac{1}{2}\|\rho^{n+1}\|^{2}-\left(\rho^{n+1},A\right)=\frac{1}{2}\|\rho^{n}\|^{2}-\left(\rho^{n},A\right)\\\
&+\Delta
t_{n}\left(\nabla{\cdot\,}{\boldsymbol{u}}^{n+1/2},\frac{1}{2}\left(\rho^{n+1/2}-A\right)^{2}-\alpha_{\rho}\left(\rho^{n+1/2}-\bar{\rho}^{n+1/2}\right)\left(\rho^{n+1/2}-A\right)\right).\end{split}$
(5.7)

There exists an $A$ such that
$\left(\frac{1}{2}\left(\rho^{n+1/2}-A\right)^{2}-\alpha_{\rho}\left(\rho^{n+1/2}-\bar{\rho}^{n+1/2}\right)\left(\rho^{n+1/2}-A\right)\right)\in{\mathcal{Q}}$
provided that $2k_{\rho}\leq k_{P}$. This means that squared density is
conserved if $2k_{\rho}\leq k_{P}$ by using the weak divergence-free condition
(5.1) and that mass is conserved.

Shift invariance: Substituting $\rho^{n+j}$ with $\rho^{n+j}+c$ and repeating
the steps of the proof of Theorem 3.1 shows that the density update (5.1) is
invariant if $\alpha_{\rho}=0$ or if
$\bar{\rho}^{n+1/2}=\frac{1}{|\Omega|}\int_{\Omega}\rho^{n+1/2}\,\mathrm{d}{\boldsymbol{x}}$.
Kinetic energy: We begin by using the identity from [8, Sec 3] and the
definition of the $L^{2}$-projection (3.18) to obtain

$\begin{split}\frac{1}{2}\|\sqrt{\rho^{n+1}}{\boldsymbol{u}}^{n+1}\|^{2}-\frac{1}{2}\|\sqrt{\rho^{n}}{\boldsymbol{u}}^{n}\|^{2}&=\left(\rho^{n+1}{\boldsymbol{u}}^{n+1}-\rho^{n}{\boldsymbol{u}}^{n},{\boldsymbol{u}}^{n+1/2}\right)-\frac{1}{2}\left(\rho^{n+1}-\rho^{n},{\boldsymbol{u}}^{n}{\cdot}{\boldsymbol{u}}^{n+1}\right)\\\
&=\left(\rho^{n+1}{\boldsymbol{u}}^{n+1}-\rho^{n}{\boldsymbol{u}}^{n},{\boldsymbol{u}}^{n+1/2}\right)-\frac{1}{2}\left(\rho^{n+1}-\rho^{n},\overline{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\boldsymbol{u}}^{n}{\cdot}{\boldsymbol{u}}^{n+1}}}\right).\end{split}$
(5.8)

Next, repeating the steps of the proof of Theorem 3.1 on (5.1) leads to

$\begin{split}\frac{1}{2}\left(\|\sqrt{\rho^{n+1}}{\boldsymbol{u}}^{n+1}\|^{2}-\|\sqrt{\rho^{n}}{\boldsymbol{u}}^{n}\|^{2}\right)=\Delta
t_{n}\left({\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\alpha_{P}+\frac{1}{2}-\alpha_{\boldsymbol{m}}+\frac{1}{2}\alpha_{\rho}}\right)\left(\left(\nabla{\cdot\,}{\boldsymbol{u}}^{n+1/2}\right){\boldsymbol{m}}^{n+1/2},{\boldsymbol{u}}^{n+1/2}\right),\end{split}$
(5.9)

which shows that kinetic energy is conserved if
$\alpha_{\boldsymbol{m}}-\alpha_{P}-\alpha_{\rho}/2=1/2$.
Momentum: Setting ${\boldsymbol{v}}={\boldsymbol{e}}_{i}$ inside (5.1) and
repeating the steps of the proof of Theorem 3.1 yields

$\begin{split}\int_{\Omega}{\boldsymbol{m}}_{i}^{n+1}\,\mathrm{d}{\boldsymbol{x}}=\int_{\Omega}{\boldsymbol{m}}_{i}^{n}\,\mathrm{d}{\boldsymbol{x}}+{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\Delta
t_{n}}\left(1-\alpha_{\boldsymbol{m}}\right)\left(\left(\nabla{\cdot\,}{\boldsymbol{u}}^{n+1/2}\right){\boldsymbol{m}}^{n+1/2},{\boldsymbol{e}}_{i}\right)\\\
+\frac{1}{2}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\Delta
t_{n}}\left({\boldsymbol{e}}_{i}{\cdot}\nabla\rho^{n+1/2},{\boldsymbol{u}}^{n+1/2}{\cdot}{\boldsymbol{u}}^{n+1/2}-{\boldsymbol{u}}^{n}{\cdot}{\boldsymbol{u}}^{n+1}\right).\end{split}$
(5.10)

This means that momentum is not conserved due to the last term in (5.10).
Angular momentum: Setting ${\boldsymbol{v}}={\boldsymbol{\phi}}_{i}$ inside
(5.1) and repeating the steps of the proof of Theorem 3.1 yields

$\begin{split}\int_{\Omega}\left({\boldsymbol{m}}{\times}{\boldsymbol{x}}\right)_{i}^{n+1}\,\mathrm{d}{\boldsymbol{x}}=\int_{\Omega}\left({\boldsymbol{m}}{\times}{\boldsymbol{x}}\right)_{i}^{n}\,\mathrm{d}{\boldsymbol{x}}+{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\Delta
t_{n}}\left(1-\alpha_{\boldsymbol{m}}\right)\left(\left(\nabla{\cdot\,}{\boldsymbol{u}}^{n+1/2}\right){\boldsymbol{m}}^{n+1/2},{\boldsymbol{\phi}}_{i}\right)\\\
+\frac{1}{2}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\Delta
t_{n}}\left({\boldsymbol{\phi}}_{i}{\cdot}\nabla\rho^{n+1/2},{\boldsymbol{u}}^{n+1/2}{\cdot}{\boldsymbol{u}}^{n+1/2}-{\boldsymbol{u}}^{n}{\cdot}{\boldsymbol{u}}^{n+1}\right).\end{split}$
(5.11)

This means that angular momentum is not conserved due to the last term in
(5.11).

###  5.1 Adaptive time-stepping

We use a simple adaptive time-stepping algorithm based on the CFL condition.
Given user-defined parameters CFL, $s_{\min}$ and $s_{\max}$, the next time
step is computed using the following adaptive algorithm:

/* Compute time step increment based on CFL condition */

$s_{\text{cfl}}=\textrm{CFL}\min_{{\boldsymbol{x}}\in\Omega}h({\boldsymbol{x}})/\left(\|{\boldsymbol{u}}^{n}\|_{l^{2}}\Delta
t_{n-1}\right)$

/* Make sure the next time step $\Delta t_{n}$ is bounded by
$s_{\text{max}}\Delta t_{n-1}$ */

$s=\min{\left(s_{\text{cfl}},s_{\text{max}}\right)}$

$\Delta t_{n}=s\Delta t_{n-1}$

if  _$s <s_{\text{min}}$_ then

Repeat previous time step with $\Delta t_{n}$ instead

end if

return _Time step $\Delta t_{n}$ and flag whether to repeat the time step or
not_

Algorithm 1 A simplified version of the time step control algorithm from [11,
Sec 5.4] used to compute the next time step $\Delta t_{n}$ given user-defined
parameters CFL, $s_{\text{max}},s_{\text{min}}$.

##  6 Numerical validation.

The goal of this section is to numerically compare the formulations from Table
2. The table includes newly developed formulations and formulations previously
known from the literature. Another goal is to numerically verify the viscous
regularization from Section 4. The nonlinear system is solved with the
relative tolerance of $10^{-14}$. Our code is implemented in FEniCS [2], an
open-source finite element library. We perform experiments using high-order
Lagrange elements in space. More specifically, the SI-MEDMAC uses
$(\rho,{\boldsymbol{u}},P)\in{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\mathbb{P}}_{2}{\times}{\mathbb{P}}_{3}{\times}{\mathbb{P}}_{2}}$
and the other formulations uses
$(\rho,{\boldsymbol{u}},P)\in{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\mathbb{P}}_{3}{\times}{\mathbb{P}}_{3}{\times}{\mathbb{P}}_{2}}$.

###  6.1 Accuracy test.

In this section, we verify the accuracy of the proposed formulations using a
manufactured solution on a unit disk. The time step is computed using CFL =
0.025 and the initial time step is used throughout the whole simulation. A
small time step is chosen so that the error in time is negligible. Since the
error will be dominated by spatial errors, the expected convergence rate is
then 4 for velocity and 3 for pressure. The expected rate for density is 4 for
all formulations except SI-MEDMAC where we expect a rate of 3. We follow the
setup from [14] where the forcing function ${\boldsymbol{f}}$ is chosen to
obtain the following exact solution

$\begin{split}&\rho({\boldsymbol{x}},t)=2+x\cos(\sin(t))+y\sin(\sin(t)),\\\
&{\boldsymbol{u}}({\boldsymbol{x}},t)=\begin{bmatrix}-y\cos(t)\\\
x\cos(t)\end{bmatrix},\\\ &p({\boldsymbol{x}},t)=\sin(x)\sin(y)\sin(t),\\\
&P({\boldsymbol{x}},t)=p-\alpha_{\rho}\rho{\boldsymbol{u}}{\cdot}{\boldsymbol{u}}-\frac{1}{2}\alpha_{\rho}\bar{\rho}{\boldsymbol{u}}{\cdot}{\boldsymbol{u}}-\frac{1}{|\Omega|}\int_{\Omega}\left(p-\alpha_{\rho}\rho{\boldsymbol{u}}{\cdot}{\boldsymbol{u}}-\frac{1}{2}\alpha_{\rho}\bar{\rho}{\boldsymbol{u}}{\cdot}{\boldsymbol{u}}\right)\,\mathrm{d}{\boldsymbol{x}}.\end{split}$
(6.1)

Note that we have included the modified pressure $P$ above which is the one
that is solved for in (5.1). Dirichlet boundary conditions are imposed for the
velocity. The boundary conditions are imposed strongly through the linear
system using the routine implemented in FEniCS [2]. We perform a convergence
study using a series of unstructured meshes. The dynamic viscosity is set to
$\mu=0.01$ and we set $\kappa=0$. The termination time is set to $T=1$.

The convergence results are presented in Figure 1 where the $L^{2}$-error is
plotted against the number of degrees of freedom (DOF). All the errors are
computed using a high-order quadrature and are relative, i.e., they are
normalized with their corresponding norm. Overall, the velocity errors
converge with the expected convergence rate 4 and all formulations have
similar errors. Similarly, the pressure errors converge with the expected rate
of 3 with similar errors for all formulations. However, for density, the
convergence rate drops to around 3.5 for all formulations except the LSI-EMAC
formulation. The reason for this is that the added term in the density update
$\alpha_{\rho}(\nabla{\cdot\,}{\boldsymbol{u}})(\rho-\bar{\rho})$ pollutes the
accuracy of $\rho$. If ${\boldsymbol{u}}$ has an error $\mathcal{O}(h^{4})$,
then $\nabla{\cdot\,}{\boldsymbol{u}}=\mathcal{O}(h^{3})$. This means that a
$\mathcal{O}(h^{3})$ perturbation has been added to the density update and
this explains why density converges with a suboptimal rate of $\approx 3.5$
instead of 4. Even though the SI-MEDMAC uses $\rho\in{\mathbb{P}}_{2}$, we
observe that it converges with the superoptimal rate $\approx 3.5$. SI-MEDMAC
is the most accurate formulation for density even though it uses a lower
polynomial degree.

Another important observation is that the shift-invariant formulations are
more accurate than their counterpart, i.e., SI-EDMAC is more accurate than
EDMAC, and SI-MEMAC is more accurate than MEMAC. More specifically, the error
for density is roughly 1 order of magnitude lower for SI compared to the non-
SI formulations.

(a) Density

(b) Velocity

(c) Modified pressure

Figure 1: Manufactured solution. Relative $L^{2}$-error versus degrees of
freedom. Convergence study for $T=1$, $\mu=0.01$, CFL = 0.025. All
formulations used
$(\rho,{\boldsymbol{u}},P)\in{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\mathbb{P}}_{3}{\times}{\mathbb{P}}_{3}{\times}{\mathbb{P}}_{2}}$,
except SI-MEDMAC which used
$(\rho,{\boldsymbol{u}},P)\in{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\mathbb{P}}_{2}{\times}{\mathbb{P}}_{3}{\times}{\mathbb{P}}_{2}}$.

###  6.2 Gresho problem with variable density.

In this section, we numerically verify the conserved properties of the
formulations under investigation (see Table 2) for the inviscid case. To this
end, we consider the so-called Gresho problem on the domain,
$\Omega=\\{(x,y)\in(-0.5,-0.5){\times}(0.5,0.5)\\}$. The problem was
originally conceived for constant density incompressible flow [10] and we
modify the problem so that density is variable and obtain the following exact
solution

$\begin{split}&\rho=5+0.1\left(1-\tanh\left(\frac{x^{2}+y^{2}}{r_{0}^{2}}-1\right)\right),\\\
&r\leq 0.2\mathrel{\mathop{\mathchar
58\relax}}\begin{dcases}{\boldsymbol{u}}=\begin{pmatrix}-5y\\\
5x\end{pmatrix}\\\ p=12.5r^{2}+C_{1}\end{dcases},\\\
&r>0.4\mathrel{\mathop{\mathchar
58\relax}}\begin{dcases}{\boldsymbol{u}}=\begin{pmatrix}0\\\ 0\end{pmatrix}\\\
p=0\end{dcases},\\\ 0.2\leq&r\leq 0.4\mathrel{\mathop{\mathchar
58\relax}}\begin{dcases}{\boldsymbol{u}}=\begin{pmatrix}\frac{-2y}{r}+5y\\\
\frac{2x}{r}-5x\end{pmatrix}\\\
p=12.5r^{2}-20r+4\log(r)+C_{2}\end{dcases},\end{split}$ (6.2)

where $r=\sqrt{x^{2}+y^{2}}$ and we set $r_{0}=1/8$ and

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

The main difficulty with this benchmark is that velocity is non-smooth at
$r=0.2$ and $r=0.4$, which will make $\nabla{\cdot\,}{\boldsymbol{u}}\neq 0$
even when $h\rightarrow 0$. This means that the problem will test if the
numerical scheme can conserve the properties considered in this work even when
the divergence-free condition is severely violated. To be able to run the
simulations longer without code crashes due to positivity issues, the density
profile was chosen to be smooth. It is important to emphasize that a
stabilized scheme for the incompressible variable density Euler equations
would be able to solve this problem without numerical instabilities. Since the
purpose of this test is to verify Table 1 numerically, no stabilization is
used. We use a coarse mesh (4705 ${\mathbb{P}}_{3}$ and 2113
${\mathbb{P}}_{2}$ nodes). The time step is computed using CFL = 0.5 and the
initial time step is used throughout the whole simulation. Periodic boundary
conditions are enforced. The periodic boundary mapping is done via built-in
functions in FEniCS [2].

The results are presented in Figure 2 and show that the proved properties
agree with the numerical results. The minimum of $\rho$ is computed node-wise
and the $L^{2}$-errors are not normalized. More specifically, we observe that
kinetic energy, momentum and angular momentum are conserved for all EMAC
formulations. There are some small violations of this due to cancellation
errors, momentum errors from the modified Crank-Nicolson method and when
density is close to the vacuum stage. Overall the errors are similar to what
has been obtained by other researchers for the constant density Gresho
problem, e.g., [5, Sec 4.1] and [16, Sec 6.2]. In Figure 2, by ‘Squared
density*’ we refer to the modified energy of the advection equation given by
$\int_{\Omega}\rho^{2}\,\mathrm{d}{\boldsymbol{x}}-\frac{1}{|\Omega|}\left(\int_{\Omega}\rho\,\mathrm{d}{\boldsymbol{x}}\right)^{2}$.

One important observation is that the shift-invariant property plays a big
role, where the shift-invariant formulations are the most robust. For
instance, compare SI-MEMAC and SI-EDMAC to MEMAC and EDMAC. The SI-EDMAC and
SI-MEDMAC formulations are the most robust and the only formulations that
avoided code crash. The SI-MEDMAC formulation is the most accurate for density
even though it used fewer degrees of freedom, i.e., 2113 ${\mathbb{P}}_{2}$
nodes compared to 4705 ${\mathbb{P}}_{3}$ nodes.

Figure 2: Gresho problem with variable density and periodic boundary
conditions using $\text{CFL}=0.5$ and 4705 ${\mathbb{P}}_{3}$ nodes.
$\kappa=\nu=\mu=0$. No stabilization is used. All formulations used
$(\rho,{\boldsymbol{u}},P)\in{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\mathbb{P}}_{3}{\times}{\mathbb{P}}_{3}{\times}{\mathbb{P}}_{2}}$,
except SI-MEDMAC which used
$(\rho,{\boldsymbol{u}},P)\in{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\mathbb{P}}_{2}{\times}{\mathbb{P}}_{3}{\times}{\mathbb{P}}_{2}}$.

###  6.3 Viscous test.

In this test, we numerically verify the properties of the viscous
regularization. Since we are only interested in the properties of the model
and not the discretization, we modify the Gresho problem from Section 6.2 so
that ${\boldsymbol{u}}$ is smooth, meaning that
$\nabla{\cdot\,}{\boldsymbol{u}}\rightarrow 0$ as $h\rightarrow 0$. This means
that for a fine mesh, the divergence-free condition will not affect the
conserved properties much. The only factor that will meaningfully affect the
computed properties is viscous regularization. We consider the following
initial condition

$\begin{split}&\rho=1+0.5\left(1-\tanh\left(\frac{x^{2}+y^{2}}{r_{0}^{2}}-1\right)\right),\\\
&{\boldsymbol{u}}=0.1\left(1-\tanh\left(\frac{x^{2}+y^{2}}{r_{0}^{2}}-1\right)\right)\begin{pmatrix}-5y\\\
5x\end{pmatrix},\end{split}$ (6.4)

on the domain $\Omega=\\{(x,y)\in(-1.25,-1.25){\times}(1.25,1.25)\\}$ using
periodic boundary conditions. The initial condition (6.4) is an exact solution
iff $\kappa=\mu=0$. For $\mu\neq 0$ or $\kappa\neq 0$, the analytic solution
is not known. We are interested in the behavior of the solution when
$\kappa>0$ and $\mu=0$ since we want to verify Theorem 4.1. If $\mu=0$, the
first condition of Theorem 4.1 is equivalent to kinetic energy conservation.

We set $r_{0}=0.125$, $\kappa=0.01$, $\mu=0$. The time step is computed using
CFL = 0.2 and the initial time step is used throughout the whole simulation.
We use a fine mesh (73728 ${\mathbb{P}}_{3}$ nodes) so that the discretization
error is small since we are only interested in the viscous contribution. Since
the discretization error is small, the formulations in Table 2 will produce
very similar results. We use the SI-MEDMAC formulation.

The results are presented in Figure 3 using three different viscous
regularizations (4.19)-(4.21). The results agree with the theory, i.e.,

  * •

${\boldsymbol{f}}_{GP}$ conserves kinetic energy and momentum but not angular
momentum.

  * •

${\boldsymbol{f}}_{KS}$ and ${\boldsymbol{f}}_{0}$ conserve momentum and
angular momentum but not kinetic energy.

Figure 3: Comparison of viscous regularizations (4.19)-(4.21) using
$\mu=\nu=0$, $\kappa=0.01$ for a smooth variant of the Gresho problem with
variable density and periodic boundary conditions. The solution was obtained
using $\text{CFL}=0.2$, 73728 ${\mathbb{P}}_{3}$ nodes and no stabilization
was used.

###  6.4 Lock-exchange 2D.

In this section, we compare some different viscous regularizations for a more
realistic problem. The problem we consider is the so-called lock-exchange
problem which was numerically studied extensively by Bartholomew and Laizet
[3] and Birman et al. [4]. Notably, this configuration has also been the
subject of physical experimental exploration by Lowe et al. [19], yielding
observations that align with the findings put forth by [4]. The reason why
this problem is of interest to us is that it is one of the few case studies in
the variable density flow literature which includes mass diffusivity in the
model, i.e., $\kappa>0$. The classical setup is a heavy fluid with density
$\rho_{1}$ and a lighter fluid with density $\rho_{2}$ separated by a vertical
barrier. The heavy fluid is located to the left of the barrier and the lighter
fluid is to the right. The barrier is removed at $t=0$ and then the heavy gas
moves to the right and the lighter fluid moves to the left. Following the
setup from [3, 4], the computational domain is set to
$\Omega=\\{(x,y)\in(-L/2,L/2){\times}(0,32L)\\}$ with a characteristic
velocity scale set to $u=\sqrt{g(\rho_{1}-\rho_{2})/\rho_{1}L}$ and the
Reynolds number is defined as
$Re=\rho_{1}L^{3/2}\sqrt{g(\rho_{1}-\rho_{2})/\rho_{1}}/\mu$. Following the
reference we set $L=1$, $\rho_{2}=0.7$ and $\rho_{1}=1.0$. The Reynolds number
is set to $4000$ yielding a constant kinematic viscosity coefficient $\nu$. In
summary, we set $\mu=\rho\nu$. The forcing function is set to
${\boldsymbol{f}}=(0,-\rho g)$ with $g=1$ to yield a downward gravitational
force. The initial condition is initially regularized using the error function
and is given by

$\rho^{0}({\boldsymbol{x}})=\frac{1}{2}\left(\frac{\rho_{2}}{\rho_{1}}+1\right)-\frac{1}{2}\left(1-\frac{\rho_{2}}{\rho_{1}}\right)\text{erf}\left(x_{0}\sqrt{Re}\right),$
(6.5)

where $x_{0}=14$ is the location of the barrier at $t=0$. Following the
reference, the mass diffusivity coefficient is set to $\kappa=\nu$ which means
that the so-called Schmidt number which governs the relation of molecular
diffusivity and kinematic viscosity is set to 1. Also following the reference,
we set slip boundary conditions on all boundaries which are imposed strongly
through the linear system. A Neumann condition
${\boldsymbol{n}}{\cdot}\nabla\rho|_{\partial\Omega}=0,$ is automatically
imposed weakly on all boundaries through the weak form (5.1).

The computed densities are presented in Figures 4 \- 6 at times $T=1,3,5,7,10$
(at time-scale $\\\
L^{1/2}\left(g(\rho_{1}-\rho_{2})/\rho_{1}\right)^{-1/2}$). The contour lines
at the final time step are plotted in Figure 7. The results were computed
using the SI-MEDMAC formulation using 1851121 ${\mathbb{P}}_{3}$ nodes on a
uniform mesh. The parameters for the time-stepping are set to
$\text{CFL}=0.1$, $s_{max}=1.02$, $s_{min}=0.98$. The results using the
momentum and angular momentum conserving viscous flux ${\boldsymbol{f}}_{0}$
match with the reference [3, Fig 4] since we use the same viscous flux on a
fine mesh. In B we extend the benchmark to 3D.

|
---|---

Figure 4: Momentum and angular momentum conserving viscous flux:
${\boldsymbol{f}}_{0}$. Lock-exchange problem at $Re=4000$, Schmidt number 1
and density ratio $1/0.7$. Density $\rho$ at times $t=1,3,5,7,10$.

|
---|---

Figure 5: Momentum and angular momentum conserving viscous flux:
${\boldsymbol{f}}_{KS}$. Lock-exchange problem at $Re=4000$, Schmidt number 1
and density ratio $1/0.7$. Density $\rho$ at times $t=1,3,5,7,10$.

|
---|---

Figure 6: Guermond-Popov viscous flux: ${\boldsymbol{f}}_{GP}$. Lock-exchange
problem at $Re=4000$, Schmidt number 1 and density ratio $1/0.7$. Density
$\rho$ at times $t=1,3,5,7,10$.

|  |
---|---|---

(a) Momentum and angular momentum conserving viscous flux:
${\boldsymbol{f}}_{0}$.

(b) Kažihov and Smagulov [17] viscous model: ${\boldsymbol{f}}_{KS}$.

(c) Guermond-Popov viscous flux: ${\boldsymbol{f}}_{GP}$.

Figure 7: Lock-exchange problem at $Re=4000$, Schmidt number 1 and density
ratio $1/0.7$. Contour lines of density at $T=10$ for three different viscous
regularizations.

##  7 Conclusion.

By modifying the governing equations in a consistent way, we introduce a new
formulation for the incompressible variable density Navier-Stokes equations
that, when discretized by a Galerkin method, is shift-invariant and mass,
squared density, kinetic energy, momentum and angular momentum conserving (SI-
MEDMAC). The formulation we propose for the inviscid case (2.1) is

$\displaystyle{\partial}_{t}\rho+{\boldsymbol{u}}{\cdot}\nabla\rho+\frac{1}{2}(\nabla{\cdot\,}{\boldsymbol{u}})(\rho-\bar{\rho})$
$\displaystyle=0,$ (7.1)
$\displaystyle{\partial}_{t}{\boldsymbol{m}}+{\boldsymbol{u}}{\cdot}\nabla{\boldsymbol{m}}+(\nabla{\cdot\,}{\boldsymbol{u}}){\boldsymbol{m}}+\nabla
P+(\nabla{\boldsymbol{u}}){\boldsymbol{m}}$
$\displaystyle+\frac{1}{2}\left(\nabla\rho({\boldsymbol{u}}{\cdot}{\boldsymbol{u}})-\frac{1}{2}\nabla((\rho-\bar{\rho}){\boldsymbol{u}}{\cdot}{\boldsymbol{u}})\right)$
$\displaystyle={\boldsymbol{f}},$

where
$P=p-\frac{1}{4}\rho{\boldsymbol{u}}{\cdot}{\boldsymbol{u}}-\frac{1}{4}\bar{\rho}{\boldsymbol{u}}{\cdot}{\boldsymbol{u}}$
is a modified pressure and
$\bar{\rho}=\frac{1}{|\Omega|}\int_{\Omega}\rho\,\mathrm{d}{\boldsymbol{x}}$.
If $\rho=1$, the formulation simplifies to the EMAC formulation by Charnyi et
al. [5]. One important finding is that by ensuring that the polynomial degree
of the mixed FEM formulation satisfies $k_{\rho}\leq k_{P}$, mass conservation
is automatically achieved. Similarly, by setting $2k_{\rho}\leq k_{P}$,
squared density conservation can be automatically achieved. We numerically
compare the new formulation to some previously considered in the literature
(see Table 2). The new formulation performs the best on the benchmarks
considered in this manuscript, both in terms of accuracy for smooth problems
and in terms of robustness.

Lastly, we consider the effect of viscous regularizations and our theoretical
findings are summarized in Theorem 4.1. One implication of Theorem 4.1 is
that, if mass diffusion is included in the model, there is no viscous
regularization that simultaneously dissipates kinetic energy and conserves
angular momentum. Another implication is that a unique viscous regularization
dissipates kinetic energy and conserves momentum.

As additional work, one could extend the formulation to more realistic
multiphase flow models where surface tension is included such as Cahn-Hilliard
or similar. One interesting research direction would be to consider
positivity-preserving discretizations that are also shift-invariant. It would
also be interesting to investigate which viscous regularization in Section 4
gives the best match with experimental data from real-world experiments. For
underresolved flows, i.e., $\kappa\ll|{\boldsymbol{u}}|h$ or
$\mu\ll\rho|{\boldsymbol{u}}|h$, we recommend that a stabilization technique
[29, 21, 23, 24] is used.

## Acknowledgments.

The computations were enabled by resources in project SNIC 2022/22-428
provided by the Swedish National Infrastructure for Computing (SNIC) at
UPPMAX, partially funded by the Swedish Research Council through grant
agreement no. 2018-05973. The first author was partly supported by the Center
for Interdisciplinary Mathematics, Uppsala University. We extend our sincere
gratitude to the anonymous reviewers whose constructive comments greatly
increased the quality of this manuscript.

##  Appendix A Different formulations of the model problem

Since it is common to write the governing equations in primitive form, we
reiterate the results of this manuscript for the governing equations below

$\displaystyle{\partial}_{t}\rho+{\boldsymbol{u}}{\cdot}\nabla\rho$
$\displaystyle=\nabla{\cdot\,}{\boldsymbol{f}}_{\rho},$ (A.1)
$\displaystyle\rho{\partial}_{t}{\boldsymbol{u}}+{\boldsymbol{m}}{\cdot}\nabla{\boldsymbol{u}}+\nabla
p$
$\displaystyle={\boldsymbol{f}}+\nabla{\cdot\,}\left(\mu\left(\nabla{\boldsymbol{u}}+(\nabla{\boldsymbol{u}})^{\top}\right)\right)+{\boldsymbol{f}}_{{\boldsymbol{m}}}-(\nabla{\cdot\,}{\boldsymbol{f}}_{\rho}){\boldsymbol{u}},\quad$
$\displaystyle({\boldsymbol{x}},t)\in\Omega{\times}(0,T],$
$\displaystyle\nabla{\cdot\,}{\boldsymbol{u}}$ $\displaystyle=0,$
$\displaystyle{\boldsymbol{u}}({\boldsymbol{x}},0)$
$\displaystyle={\boldsymbol{u}}_{0}({\boldsymbol{x}}),$
$\displaystyle\rho({\boldsymbol{x}},0)$
$\displaystyle=\rho_{0}({\boldsymbol{x}}),\quad$
$\displaystyle{\boldsymbol{x}}\in\Omega.$

Similarly, it is also common to write the time-derivative as
$\sqrt{\rho}{\partial}_{t}\left(\sqrt{\rho}{\boldsymbol{u}}\right)$ [13]
instead of ${\partial}_{t}{\boldsymbol{m}}$ which leads to

$\displaystyle{\partial}_{t}\rho+{\boldsymbol{u}}{\cdot}\nabla\rho=$
$\displaystyle\nabla{\cdot\,}{\boldsymbol{f}}_{\rho},$ (A.2)
$\displaystyle\sqrt{\rho}{\partial}_{t}\left(\sqrt{\rho}{\boldsymbol{u}}\right)+\frac{1}{2}{\boldsymbol{u}}{\cdot}\nabla{\boldsymbol{m}}+\frac{1}{2}{\boldsymbol{m}}{\cdot}\nabla{\boldsymbol{u}}+\nabla
p=$
$\displaystyle{\boldsymbol{f}}+\nabla{\cdot\,}\left(\mu\left(\nabla{\boldsymbol{u}}+(\nabla{\boldsymbol{u}})^{\top}\right)\right)$
$\displaystyle+{\boldsymbol{f}}_{{\boldsymbol{m}}}-\frac{1}{2}(\nabla{\cdot\,}{\boldsymbol{f}}_{\rho}){\boldsymbol{u}},\
$ $\displaystyle({\boldsymbol{x}},t)\in\Omega{\times}(0,T],$
$\displaystyle\nabla{\cdot\,}{\boldsymbol{u}}=$ $\displaystyle 0,$
$\displaystyle{\boldsymbol{u}}({\boldsymbol{x}},0)=$
$\displaystyle{\boldsymbol{u}}_{0}({\boldsymbol{x}}),$
$\displaystyle\rho({\boldsymbol{x}},0)=$
$\displaystyle\rho_{0}({\boldsymbol{x}}),$
$\displaystyle{\boldsymbol{x}}\in\Omega.$

Many Navier-Stokes solvers use these formulations as their starting point and
all the theoretical results of this manuscript directly translate to these
forms provided that the formulations are modified accordingly. For instance,
Theorem 4.1 holds for the above formulations since the term
$(\nabla{\cdot\,}{\boldsymbol{f}}_{\rho}){\boldsymbol{u}}$ is included.
Additionally, all the conservativeness given when
$\nabla{\cdot\,}{\boldsymbol{u}}\neq 0$ holds provided the formulations are
slightly modified, see the forthcoming sections.

###  A.1 Velocity formulation.

If solving the problem in primitive form (A.1) is preferred we propose the
following formulation: Find $\\\
(\rho,{\boldsymbol{m}},P)\in{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\mathcal{M}}{\times}\mathbf{{\mathbfcal
V}}{\times}{\mathcal{Q}}}$ such that

$\displaystyle({\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\partial}_{t}}\rho,w)+({\boldsymbol{u}}{\cdot}\nabla\rho,w)+\alpha_{\rho}((\nabla{\cdot\,}{\boldsymbol{u}})(\rho-\bar{\rho}),w)+(\kappa\nabla\rho,\nabla
w)$ $\displaystyle=0,\quad$ $\displaystyle\forall w\in{\mathcal{M}},$ (A.3)
$\displaystyle(\rho{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\partial}_{t}}{\boldsymbol{u}},{\boldsymbol{v}})+b({\boldsymbol{m}},{\boldsymbol{u}},{\boldsymbol{v}})+\alpha_{\boldsymbol{m}}((\nabla{\cdot\,}{\boldsymbol{u}}){\boldsymbol{m}},{\boldsymbol{u}})+(\nabla
P,{\boldsymbol{v}})$
$\displaystyle+\left(\alpha_{P}+\frac{1}{2}\right)b({\boldsymbol{v}},{\boldsymbol{m}},{\boldsymbol{u}})+\left(\alpha_{P}-\frac{1}{2}\right)b({\boldsymbol{v}},{\boldsymbol{u}},{\boldsymbol{m}})-\frac{1}{2}\alpha_{\rho}(\nabla({\boldsymbol{m}}{\cdot}{\boldsymbol{u}}),{\boldsymbol{v}})$
$\displaystyle-\frac{1}{2}\left(({\boldsymbol{v}}{\cdot}\nabla\rho,{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\boldsymbol{u}}{\cdot}{\boldsymbol{u}}})-\alpha_{\rho}\left(\nabla\left(\left(\rho-\bar{\rho}\right){\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\boldsymbol{u}}{\cdot}{\boldsymbol{u}}}\right),{\boldsymbol{v}}\right)\right)+\alpha_{\rho}((\nabla{\cdot\,}{\boldsymbol{u}})\bar{\rho}{\boldsymbol{u}},{\boldsymbol{v}})+\alpha_{\rho}(\nabla\left(\bar{\rho}{\boldsymbol{u}}{\cdot}{\boldsymbol{u}}\right),{\boldsymbol{v}})$
$\displaystyle+\left(\mu\left(\nabla{\boldsymbol{u}}+(\nabla{\boldsymbol{u}})^{\top}\right),\nabla{\boldsymbol{v}}\right)$
$\displaystyle-({\boldsymbol{f}},{\boldsymbol{v}})-({\boldsymbol{f}}_{m},{\boldsymbol{v}})-b(\kappa\nabla\rho,{\boldsymbol{u}},{\boldsymbol{v}})-b(\kappa\nabla\rho,{\boldsymbol{v}},{\boldsymbol{u}})$
$\displaystyle=0,\quad$
$\displaystyle\forall{\boldsymbol{v}}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\in}{\mathbfcal
V},$ $\displaystyle(\nabla{\cdot\,}{\boldsymbol{u}},q)$
$\displaystyle=0,\quad$ $\displaystyle\forall q\in{\mathcal{Q}},$

where
$P=p-\alpha_{P}\rho{\boldsymbol{u}}{\cdot}{\boldsymbol{u}}-\frac{1}{2}\alpha_{\rho}\bar{\rho}{\boldsymbol{u}}{\cdot}{\boldsymbol{u}}$.
Note that integration by parts was performed on
$(\nabla{\cdot\,}{\boldsymbol{f}}_{\rho}){\boldsymbol{u}}$ inside (A.1). The
properties of (A.3) are summarized in Table 4 and Table 1.

###  A.2 Fully discrete approximation.

Similar to Section 5 the standard Crank-Nicolson method can be applied on
(A.3).
$\rho{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\partial}_{t}}{\boldsymbol{u}}$
is discretized as
$\frac{\rho^{n+1/2}\left({\boldsymbol{u}}^{n+1}-{\boldsymbol{u}}^{n}\right)}{\Delta
t_{n}}$ and
${\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\boldsymbol{u}}{\cdot}{\boldsymbol{u}}}$
is discretized as
${\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\boldsymbol{u}}^{n+1/2}{\cdot}{\boldsymbol{u}}^{n+1/2}}$.
This leads to all properties at the fully discrete level except kinetic
energy. By instead discretizing
${\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\boldsymbol{u}}{\cdot}{\boldsymbol{u}}}$
and
${\frac{1}{2}\left(|{\boldsymbol{u}}^{n}|^{2}+|{\boldsymbol{u}}^{n+1}|^{2}\right)}$
all properties except momentum and angular momentum are fulfilled.

###  A.3 Mixed formulation.

If solving the problem in the mixed form (A.2) is preferred we propose the
following formulation: Find
$(\rho,{\boldsymbol{m}},P)\in{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\mathcal{M}}{\times}\mathbf{{\mathbfcal
V}}{\times}{\mathcal{Q}}}$ such that

$\displaystyle({\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\partial}_{t}}\rho,w)+({\boldsymbol{u}}{\cdot}\nabla\rho,w)+\alpha_{\rho}((\nabla{\cdot\,}{\boldsymbol{u}})(\rho-\bar{\rho}),w)$
$\displaystyle=0,\quad$ $\displaystyle\forall w\in{\mathcal{M}},$ (A.4)
$\displaystyle(\sqrt{\rho}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\partial}_{t}}(\sqrt{\rho}{\boldsymbol{u}}),{\boldsymbol{v}})+\frac{1}{2}(b({\boldsymbol{m}},{\boldsymbol{u}},{\boldsymbol{v}})+b({\boldsymbol{u}},{\boldsymbol{m}},{\boldsymbol{v}}))+\alpha_{\boldsymbol{m}}((\nabla{\cdot\,}{\boldsymbol{u}}){\boldsymbol{m}},{\boldsymbol{u}})+(\nabla
P,{\boldsymbol{v}})$
$\displaystyle+\alpha_{P}b({\boldsymbol{v}},{\boldsymbol{m}},{\boldsymbol{u}})+\alpha_{P}b({\boldsymbol{v}},{\boldsymbol{u}},{\boldsymbol{m}})+\frac{1}{2}\alpha_{\rho}((\nabla{\cdot\,}{\boldsymbol{u}})\bar{\rho}{\boldsymbol{u}},{\boldsymbol{v}})+\frac{1}{2}\alpha_{\rho}(\nabla\left(\bar{\rho}{\boldsymbol{u}}{\cdot}{\boldsymbol{u}}\right),{\boldsymbol{v}})$
$\displaystyle+\left(\mu\left(\nabla{\boldsymbol{u}}+(\nabla{\boldsymbol{u}})^{\top}\right),\nabla{\boldsymbol{v}}\right)$
$\displaystyle=({\boldsymbol{f}},{\boldsymbol{v}})+({\boldsymbol{f}}_{m},{\boldsymbol{v}})+\frac{1}{2}\left(b(\kappa\nabla\rho,{\boldsymbol{u}},{\boldsymbol{v}})+b(\kappa\nabla\rho,{\boldsymbol{v}},{\boldsymbol{u}})\right)$
$\displaystyle=0,\quad$
$\displaystyle\forall{\boldsymbol{v}}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\in}{\mathbfcal
V},$ $\displaystyle(\nabla{\cdot\,}{\boldsymbol{u}},q)$
$\displaystyle=0,\quad$ $\displaystyle\forall q\in{\mathcal{Q}},$

where
$P=p-\alpha_{P}{\boldsymbol{m}}{\cdot}{\boldsymbol{u}}-\frac{1}{2}\alpha_{\rho}\bar{\rho}{\boldsymbol{u}}{\cdot}{\boldsymbol{u}}$.
Note that integration by parts was performed on
$(\nabla{\cdot\,}{\boldsymbol{f}}_{\rho}){\boldsymbol{u}}$ inside (A.1). The
properties of (A.4) are summarized in Table 4 and Table 1.

###  A.4 Fully discrete approximation

Similar to Section 5 the standard Crank-Nicolson method can be applied on
(A.4).
$\sqrt{\rho}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\partial}_{t}}(\sqrt{\rho}{\boldsymbol{u}})$
is discretized as
$\rho^{n+1/2}\frac{{\boldsymbol{u}}^{n+1}-{\boldsymbol{u}}^{n}}{\Delta
t_{n}}+\frac{1}{2}{\boldsymbol{u}}^{n+1/2}\frac{\rho^{n+1}-\rho^{n}}{\Delta
t_{n}}$ and this leads to all properties at the fully discrete level except
kinetic energy.

Table 4: Semi-discrete conservation properties for the inviscid model problem. Note that the conditions on $\alpha_{\rho},\alpha_{\boldsymbol{m}},\alpha_{P}$ are different for the different formulations. | ${\partial}_{t}\int_{\Omega}\rho{\boldsymbol{u}}{\cdot}{\boldsymbol{u}}\,\mathrm{d}{\boldsymbol{x}}=0$ | ${\partial}_{t}\int_{\Omega}{\boldsymbol{m}}\,\mathrm{d}{\boldsymbol{x}}=0$ | ${\partial}_{t}\int_{\Omega}{\boldsymbol{m}}{\times}{\boldsymbol{x}}\,\mathrm{d}{\boldsymbol{x}}=0$
---|---|---|---
${\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\partial}_{t}}{\boldsymbol{m}}\eqref{eq:mom_update_general}$ | $\alpha_{\boldsymbol{m}}-\alpha_{P}-\alpha_{\rho}/2=1/2$ | $\alpha_{\boldsymbol{m}}=1$ | $\alpha_{\boldsymbol{m}}=1$
$\sqrt{\rho}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\partial}_{t}}(\sqrt{\rho}{\boldsymbol{u}})$ (A.4) | $\alpha_{\boldsymbol{m}}-\alpha_{P}=1/2$ | $\alpha_{\rho}/2+\alpha_{\boldsymbol{m}}=1$ | $\alpha_{\rho}/2+\alpha_{\boldsymbol{m}}=1$
$\rho{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\partial}_{t}}{\boldsymbol{u}}$ (A.3) | $\alpha_{\boldsymbol{m}}-\alpha_{P}+\alpha_{\rho}/2=1/2$ | $\alpha_{\rho}+\alpha_{\boldsymbol{m}}=1$ | $\alpha_{\rho}+\alpha_{\boldsymbol{m}}=1$

##  Appendix B Lock-exchange 3D.

In this section, we consider the lock-exchange problem in 3D. Even though the
theory is developed for the Crank-Nicolson scheme, one can use the viscous
regularization and formulations developed in this work with different time-
integration schemes such as high-order (backward differentiation formula) BDF
methods. To this end, we consider the variable time-step BDF4 method from [21]
which we modify to be fully implicit.

We follow the setup of Bartholomew and Laizet [3]. The setup is similar to the
2D case with a few exceptions. No-slip boundary conditions are used on the top
and bottom boundaries and slip is used on the remaining boundaries. These
boundary conditions are imposed strongly through the linear systems. A Neumann
condition ${\boldsymbol{n}}{\cdot}\nabla\rho|_{\partial\Omega}=0,$ is
automatically imposed weakly on all boundaries through the weak form (5.1).
The Reynolds number is computed in the same as in 2D, but now it is set to
$Re=2236$ instead. The domain is set to
$\Omega=\\{(x,y)\in(-1L,17L){\times}(0,2L){\times}(0,2L)\\}$.

A random perturbation is added to the initial velocity

${\boldsymbol{u}}(0,{\boldsymbol{x}})=X\cdot
0.05\cdot\exp{(-25x^{2})}+(0,0,0)^{\top},$ (B.1)

where $X\sim U(-1,1)$, a uniform distribution of random numbers in $(-1,1)$.
The authors of [3] added this perturbation ‘to simulate that a physical
barrier is removed’.

We compute the solution using 15,696,120 ${\mathbb{P}}_{3}$ nodes on a uniform
mesh which corresponds to a similar mesh resolution as Bartholomew and Laizet
[3]. We also set $\text{CFL}=0.1$ and use the LSI-EMAC formulation using
${\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\mathbb{P}}_{3}{\times}{\mathbb{P}}_{3}{\times}{\mathbb{P}}_{2}}$
elements. The nonlinear system is solved using the relative tolerance
$10^{-7}$. The results are presented at $T=20$ in Figure 8. In the figure, we
compare the Guermond-Popov viscous flux (${\boldsymbol{f}}_{GP}$) and the
momentum and angular momentum conserving flux that Bartholomew and Laizet [3]
use (${\boldsymbol{f}}_{0}$). The simulations closely match Bartholomew and
Laizet [3], but there is no pointwise match due to the random perturbations in
the initial condition and differences in the numerical method.

Figure 8: Lock-exchange problem in 3D at $Re=2236$ and Schmidt number 1.
Density $\rho$ at $T=20$. Left: Momentum and angular momentum conserving
viscous flux ${\boldsymbol{f}}_{0}$. Right: Guermond-Popov viscous flux
${\boldsymbol{f}}_{GP}$.

## References

## References

  * Almgren et al. [1998] A. S. Almgren, J. B. Bell, P. Colella, L. H. Howell, and M. L. Welcome.  A conservative adaptive projection method for the variable density incompressible Navier-Stokes equations.  _J. Comput. Phys._ , 142(1):1–46, 1998.  ISSN 0021-9991.  doi: 10.1006/jcph.1998.5890.  URL https://doi.org/10.1006/jcph.1998.5890.
  * Alnaes et al. [2015] M. Alnaes, J. Blechta, J. Hake, A. Johansson, B. Kehlet, A. Logg, C. Richardson, J. Ring, M. Rognes, and G. Wells.  The FEniCS project version 1.5.  _Arch. Numer. Softw._ , 3:9–23, 01 2015.  doi: 10.11588/ans.2015.100.20553.  URL https://doi.org/10.11588/ans.2015.100.20553.
  * Bartholomew and Laizet [2019] P. Bartholomew and S. Laizet.  A new highly scalable, high-order accurate framework for variable-density flows: application to non-Boussinesq gravity currents.  _Comput. Phys. Commun._ , 242:83–94, 2019.  ISSN 0010-4655.  doi: 10.1016/j.cpc.2019.03.019.  URL https://doi.org/10.1016/j.cpc.2019.03.019.
  * Birman et al. [2005] V. K. Birman, J. E. Martin, and E. Meiburg.  The non-Boussinesq lock-exchange problem. II. High-resolution simulations.  _J. Fluid Mech._ , 537:125–144, 2005.  ISSN 0022-1120.  doi: 10.1017/S0022112005005033.  URL https://doi.org/10.1017/S0022112005005033.
  * Charnyi et al. [2017] S. Charnyi, T. Heister, M. A. Olshanskii, and L. G. Rebholz.  On conservation laws of Navier-Stokes Galerkin discretizations.  _J. Comput. Phys._ , 337:289–308, 2017.  ISSN 0021-9991.  doi: 10.1016/j.jcp.2017.02.039.  URL https://doi.org/10.1016/j.jcp.2017.02.039.
  * Charnyi et al. [2019] S. Charnyi, T. Heister, M. A. Olshanskii, and L. G. Rebholz.  Efficient discretizations for the EMAC formulation of the incompressible Navier-Stokes equations.  _Appl. Numer. Math._ , 141:220–233, 2019.  ISSN 0168-9274.  doi: 10.1016/j.apnum.2018.11.013.  URL https://doi.org/10.1016/j.apnum.2018.11.013.
  * Evans and Hughes [2013] J. A. Evans and T. J. R. Hughes.  Isogeometric divergence-conforming B-splines for the unsteady Navier-Stokes equations.  _J. Comput. Phys._ , 241:141–167, 2013.  ISSN 0021-9991.  doi: 10.1016/j.jcp.2013.01.006.  URL https://doi.org/10.1016/j.jcp.2013.01.006.
  * Gawlik and Gay-Balmaz [2020] E. S. Gawlik and F. Gay-Balmaz.  A conservative finite element method for the incompressible Euler equations with variable density.  _J. Comput. Phys._ , 412:109439, 15, 2020.  ISSN 0021-9991.  doi: 10.1016/j.jcp.2020.109439.  URL https://doi.org/10.1016/j.jcp.2020.109439.
  * Girault and Raviart [1986] V. Girault and P.-A. Raviart.  _Finite element methods for Navier-Stokes equations_ , volume 5 of _Springer Series in Computational Mathematics_.  Springer-Verlag, Berlin, 1986.  ISBN 3-540-15796-4.  doi: 10.1007/978-3-642-61623-5.  URL https://doi.org/10.1007/978-3-642-61623-5.  Theory and algorithms.
  * Gresho and Chan [1990] P. M. Gresho and S. T. Chan.  On the theory of semi-implicit projection methods for viscous incompressible flow and its implementation via a finite element method that also introduces a nearly consistent mass matrix. II. Implementation.  volume 11, pages 621–659. 1990.  doi: 10.1002/fld.1650110510.  URL https://doi.org/10.1002/fld.1650110510.  Computational methods in flow analysis (Okayama, 1988).
  * Guermond and Minev [2019] J.-L. Guermond and P. Minev.  High-order adaptive time stepping for the incompressible Navier-Stokes equations.  _SIAM J. Sci. Comput._ , 41(2):A770–A788, 2019\.  ISSN 1064-8275.  doi: 10.1137/18M1209301.  URL https://doi.org/10.1137/18M1209301.
  * Guermond and Popov [2014] J.-L. Guermond and B. Popov.  Viscous regularization of the Euler equations and entropy principles.  _SIAM J. Appl. Math._ , 74(2):284–305, 2014.  ISSN 0036-1399.  doi: 10.1137/120903312.  URL https://doi.org/10.1137/120903312.
  * Guermond and Quartapelle [2000] J.-L. Guermond and L. Quartapelle.  A projection FEM for variable density incompressible flows.  _J. Comput. Phys._ , 165(1):167–188, 2000.  ISSN 0021-9991.  doi: 10.1006/jcph.2000.6609.  URL https://doi.org/10.1006/jcph.2000.6609.
  * Guermond and Salgado [2009] J.-L. Guermond and A. Salgado.  A splitting method for incompressible flows with variable density based on a pressure Poisson equation.  _J. Comput. Phys._ , 228(8):2834–2846, 2009.  ISSN 0021-9991.  doi: 10.1016/j.jcp.2008.12.036.  URL http://dx.doi.org/10.1016/j.jcp.2008.12.036.
  * Haynes [2015] W. M. Haynes.  _CRC Handbook of Chemistry and Physics 2015-2016_.  CRC Handbook of Chemistry & Physics. CRC Press, 96 edition, 2015.  ISBN 1482260964; 9781482260960.  URL libgen.li/file.php?md5=c7c8a202cf41d366ba493814f9966517.
  * Ingimarson et al. [2023] S. Ingimarson, M. Neda, L. G. 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(2):176–198, 2023.  ISSN 2617-8710.  doi: https://doi.org/10.4208/ijnam2023-1008.  URL http://global-sci.org/intro/article_detail/ijnam/21353.html.
  * Kažihov and Smagulov [1977] A. V. Kažihov and v. Smagulov.  The correctness of boundary value problems in a certain diffusion model of an inhomogeneous fluid.  _Dokl. Akad. Nauk SSSR_ , 234(2):330–332, 1977\.  ISSN 0002-3264.
  * Larson and Bengzon [2013] M. G. Larson and F. Bengzon.  _The finite element method: theory, implementation, and applications_ , volume 10 of _Texts in Computational Science and Engineering_.  Springer, Heidelberg, 2013.  ISBN 978-3-642-33286-9; 978-3-642-33287-6.  doi: 10.1007/978-3-642-33287-6.  URL https://doi.org/10.1007/978-3-642-33287-6.
  * Lowe et al. [2005] R. J. Lowe, J. W. Rottman, and P. F. Linden.  The non-Boussinesq lock-exchange problem. Part 1. Theory and experiments.  _J. Fluid Mech._ , 537:101–124, 2005.  doi: 10.1017/S0022112005005069.
  * Lundgren and Nazarov [2022a] L. Lundgren and M. Nazarov.  A high-order artificial compressibility method based on Taylor series time-stepping for variable density flow.  _J. Comput. Appl. Math._ , page 114846, 2022a.  ISSN 0377-0427.  doi: https://doi.org/10.1016/j.cam.2022.114846.  URL https://www.sciencedirect.com/science/article/pii/S0377042722004447.
  * Lundgren and Nazarov [2022b] L. Lundgren and M. Nazarov.  A high-order residual-based viscosity finite element method for incompressible variable density flow, 2022b.  URL https://ssrn.com/abstract=4268034.
  * Manzanero et al. [2020] J. Manzanero, G. Rubio, D. A. Kopriva, E. Ferrer, and E. Valero.  An entropy-stable discontinuous Galerkin approximation for the incompressible Navier-Stokes equations with variable density and artificial compressibility.  _J. Comput. Phys._ , 408:109241, 32, 2020.  ISSN 0021-9991.  doi: 10.1016/j.jcp.2020.109241.  URL https://doi.org/10.1016/j.jcp.2020.109241.
  * Nazarov [2013] M. Nazarov.  Convergence of a residual based artificial viscosity finite element method.  _Comput. Math. Appl._ , 65(4):616–626, 2013.  ISSN 0898-1221.  doi: 10.1016/j.camwa.2012.11.003.  URL http://dx.doi.org/10.1016/j.camwa.2012.11.003.
  * Nazarov and Hoffman [2013] M. Nazarov and J. Hoffman.  Residual-based artificial viscosity for simulation of turbulent compressible flow using adaptive finite element methods.  _Internat. J. Numer. Methods Fluids_ , 71(3):339–357, 2013.  ISSN 0271-2091.  doi: 10.1002/fld.3663.  URL http://dx.doi.org/10.1002/fld.3663.
  * Olshanskii and Rebholz [2020] M. A. Olshanskii and L. G. Rebholz.  Longer time accuracy for incompressible Navier-Stokes simulations with the EMAC formulation.  _Comput. Methods Appl. Mech. Engrg._ , 372:113369, 17, 2020\.  ISSN 0045-7825.  doi: 10.1016/j.cma.2020.113369.  URL https://doi.org/10.1016/j.cma.2020.113369.
  * Pyo and Shen [2007] J.-H. Pyo and J. Shen.  Gauge-Uzawa methods for incompressible flows with variable density.  _J. Comput. Phys._ , 221(1):181–197, 2007.  ISSN 0021-9991.  doi: 10.1016/j.jcp.2006.06.013.  URL https://doi.org/10.1016/j.jcp.2006.06.013.
  * Raviart and Thomas [1977] P.-A. Raviart and J. M. Thomas.  A mixed finite element method for 2nd order elliptic problems.  In _Mathematical aspects of finite element methods (Proc. Conf., Consiglio Naz. delle Ricerche (C.N.R.), Rome, 1975)_ , Lecture Notes in Math., Vol. 606, pages 292–315. Springer, Berlin, 1977.
  * Scott and Vogelius [1985] L. R. Scott and M. Vogelius.  Conforming finite element methods for incompressible and nearly incompressible continua.  In _Large-scale computations in fluid mechanics, Part 2 (La Jolla, Calif., 1983)_ , volume 22 of _Lectures in Appl. Math._ , pages 221–244. Amer. Math. Soc., Providence, RI, 1985.  doi: 10.1051/m2an/1985190101111.  URL https://doi.org/10.1051/m2an/1985190101111.
  * Stiernström et al. [2021] V. Stiernström, L. Lundgren, M. Nazarov, and K. Mattsson.  A residual-based artificial viscosity finite difference method for scalar conservation laws.  _J. Comput. Phys._ , 430:Paper No. 110100, 29, 2021.  ISSN 0021-9991.  doi: 10.1016/j.jcp.2020.110100.  URL https://doi.org/10.1016/j.jcp.2020.110100.
  * Zhang et al. [2023] Y. Zhang, H. Dong, and K. Wang.  Mass, momentum and energy identical-relation-preserving scheme for the Navier-Stokes equations with variable density.  _Comput. Math. Appl._ , 137:73—-92, 2023.  ISSN 0898-1221.  doi: 10.1016/j.camwa.2023.02.004.  URL https://doi.org/10.1016/j.camwa.2023.02.004.