A hybridizable discontinuous Galerkin method for the Navier--Stokes equations with pointwise divergence-free velocity field

TL;DR

This paper presents a hybridizable discontinuous Galerkin method for the incompressible Navier--Stokes equations that ensures the velocity field is divergence-free at every point, improving stability and conservation properties.

Contribution

It introduces a modified function space approach enabling a simple, divergence-free velocity approximation that is momentum conserving, energy stable, and pressure-robust.

Findings

Method achieves pointwise divergence-free velocity fields.

Numerical examples confirm stability and accuracy in 2D and 3D.

Applicable for various polynomial approximation orders.

Abstract

We introduce a hybridizable discontinuous Galerkin method for the incompressible Navier--Stokes equations for which the approximate velocity field is pointwise divergence-free. The method builds on the method presented by Labeur and Wells [SIAM J. Sci. Comput., vol. 34 (2012), pp. A889--A913]. We show that with modifications of the function spaces in the method of Labeur and Wells it is possible to formulate a simple method with pointwise divergence-free velocity fields which is momentum conserving, energy stable, and pressure-robust. Theoretical results are supported by two- and three-dimensional numerical examples and for different orders of polynomial approximation.

Tables

Table 1. Table 1: Summary of the properties of the method of  (Labeur and Wells, 2012 ) and the proposed method of this paper. The skew-symmetric and divergence forms refer to different formulations of the momentum equation. In  (Labeur and Wells, 2012 ) both an equal- and mixed-order velocity-pressure approximation are introduced.

Formulation	mass	momentum	energy	pressure
	conserving	conserving	stable	robust
Equal order (Labeur and Wells, 2012)	$\times$	$✓$	only in	$\times$
			skew-symmetric
			form
Mixed order (Labeur and Wells, 2012)	$✓$	only in	only in	$\times$
		divergence	skew-symmetric
		form	form
Proposed method	$✓$	$✓$	$✓$	$✓$

Table 2. Table 2: Computed velocity, pressure and divergence errors in the L 2 superscript 𝐿 2 L^{2} norm for the HDG method applied to the Kovasznay problem.

	$k=2$
Cells	$\parallel u_h-u\parallel$	rate	$\parallel p_h-p\parallel$	rate	$\parallel \nabla \cdot u_h\parallel$
64	1.8e-2	-	1.6e-2	-	3.8e-14
256	2.2e-3	3.0	4.0e-3	2.0	6.7e-14
1024	2.8e-4	3.0	9.8e-4	2.0	1.3e-13
4096	3.5e-5	3.0	2.4e-4	2.0	2.5e-13
	$k=3$
Cells	$\parallel u_h-u\parallel$	rate	$\parallel p_h-p\parallel$	rate	$\parallel \nabla \cdot u_h\parallel$
64	1.4e-3	-	2.0e-3	-	1.9e-13
256	9.4e-5	3.9	2.0e-4	3.3	6.1e-13
1024	5.8e-6	4.0	2.3e-5	3.1	7.8e-13
4096	3.6e-7	4.0	2.8e-6	3.1	1.6e-12

Table 3. Table 3: Computed errors in the L 2 superscript 𝐿 2 L^{2} norm for the HDG method with the position-dependent Coriolis forcing term with different viscosity values. Note that the pressure error does not depend on the viscosity.

	$\nu =0.001$, $k=2$				$\nu =1$, $k=2$
Cells	$\parallel u_h-u\parallel$	$\parallel \nabla \cdot u_h\parallel$	$\parallel p_h-p\parallel$	rate	$\parallel u_h-u\parallel$	$\parallel \nabla \cdot u_h\parallel$	$\parallel p_h-p\parallel$	rate
64	3.7e-15	4.6e-14	9.0e-4	-	1.3e-14	4.6e-14	9.0e-4	-
256	5.1e-15	9.3e-14	2.3e-4	2.0	9.9e-15	9.2e-14	2.3e-4	2.0
1024	6.4e-15	1.8e-13	5.6e-5	2.0	7.8e-14	1.9e-13	5.6e-5	2.0
4096	2.2e-14	6.3e-13	1.4e-5	2.0	2.3e-13	3.7e-13	1.4e-5	2.0

Table 4. Table 4: Computed errors in the L 2 superscript 𝐿 2 L^{2} norm for the HDG method with different viscosity values. A comparison of the proposed method with the method of  (Labeur and Wells, 2012 ) . Note that the errors in the velocity for the proposed method do not depend on ν 𝜈 \nu , in contrast with the method of  (Labeur and Wells, 2012 ) .

	Proposed Method: $\nu =0.001$
Cells	$\parallel u_h-u\parallel$	rate	$\parallel \nabla \left(u_h-u\right)\parallel$	rate	$\parallel p_h-p\parallel$	rate	$\parallel \nabla \cdot u_h\parallel$
128	4.2e-6		2.8e-4		4.5e-4		1.5e-14
512	2.4e-7	4.1	3.4e-5	3.0	5.7e-5	3.0	4.3e-14
2048	1.4e-8	4.1	4.2e-6	3.0	7.2e-6	3.0	2.3e-14
8192	8.5e-10	4.1	5.2e-7	3.0	9.0e-7	3.0	2.4e-14
	Proposed Method: $\nu =1$
Cells	$\parallel u_h-u\parallel$	rate	$\parallel \nabla \left(u_h-u\right)\parallel$	rate	$\parallel p_h-p\parallel$	rate	$\parallel \nabla \cdot u_h\parallel$
128	4.2e-6	-	2.8e-4	-	6.5e-4	-	2.2e-15
512	2.4e-7	4.1	3.4e-5	3.0	7.6e-5	3.1	4.5e-15
2048	1.4e-8	4.1	4.2e-6	3.0	9.1e-6	3.1	8.8e-15
8192	8.5e-10	4.1	5.2e-7	3.0	1.1e-6	3.0	1.8e-14
	Method of (Labeur and Wells, 2012): $\nu =0.001$
Cells	$\parallel u_h-u\parallel$	rate	$\parallel \nabla \left(u_h-u\right)\parallel$	rate	$\parallel p_h-p\parallel$	rate	$\parallel \nabla \cdot u_h\parallel$
128	6.4e-4	-	5.3e-2	-	4.6e-4	-	2.9e-14
512	4.1e-5	4.0	6.7e-3	3.0	5.8e-5	3.0	4.8e-14
2048	2.6e-6	4.0	8.4e-4	3.0	7.2e-6	3.0	1.7e-14
8192	1.6e-7	4.0	1.1e-4	3.0	9.0e-7	3.0	2.3e-14
	Method of (Labeur and Wells, 2012): $\nu =1$
Cells	$\parallel u_h-u\parallel$	rate	$\parallel \nabla \left(u_h-u\right)\parallel$	rate	$\parallel p_h-p\parallel$	rate	$\parallel \nabla \cdot u_h\parallel$
128	3.9e-6	-	2.8e-4	-	6.2e-4	-	2.2e-15
512	2.3e-7	4.1	3.4e-5	3.0	7.3e-5	3.1	4.4e-15
2048	1.4e-8	4.1	4.2e-6	3.0	8.8e-6	3.1	8.9e-15
8192	8.3e-10	4.0	5.2e-7	3.0	1.1e-6	3.0	1.8e-14

Equations

$${\partial }_{t}u+\nabla \cdot \sigma$$

$$\nabla \cdot u$$

$$\sigma \text{mathchar}58={\sigma }_a+{\sigma }_d\text{mbox}with{\sigma }_a\text{mathchar}58=u\otimes u\text{mbox}and{\sigma }_d\text{mathchar}58=p\mathbb{I}-\nu \nabla u,$$

$$u$$

$$\sigma \cdot n-\max \left(u\cdot n,0\right)u$$

$$u(x,0)$$

$${\partial }_{t}u+(1-\chi )u\cdot \nabla u+\chi \nabla \cdot {\sigma }_a+\nabla \cdot {\sigma }_d=f,$$

$$V_h$$

$${\bar{V}}_h$$

$$Q_h$$

$${\bar{Q}}_h$$

$$0=$$

$$\begin{gathered} {\int }_{\Omega }f\cdot v_{h}dx={\int }_{\Omega }{\partial }_t{u}_h\cdot v_{h}dx-\sum _K{\int }_K{\sigma }_h\text{mathchar}58\nabla v_{h}dx+\sum _K{\int }_{\partial K}{\hat{\sigma }}_h\text{mathchar}58\left(v_h\otimes n\right)ds \\ +\sum _K{\int }_{\partial K}\nu \left(\left({\bar{u}}_h-u_h\right)\otimes n\right)\text{mathchar}58\nabla v_{h}ds\forall v_h\in V_h, \end{gathered}$$

$${\int }_{{\Gamma }_N}h\cdot {\bar{v}}_{h}ds=\sum _K{\int }_{\partial K}{\hat{\sigma }}_h\text{mathchar}58\left({\bar{v}}_h\otimes n\right)ds-{\int }_{{\Gamma }_N}\left(1-\lambda \right)\left({\bar{u}}_h\cdot n\right){\bar{u}}_h\cdot {\bar{v}}_{h}ds\forall {\bar{v}}_h\in {\bar{V}}_h,$$

$${\hat{\sigma }}_{a,h}\text{mathchar}58={\sigma }_{a,h}+\left({\bar{u}}_h-u_h\right)\otimes \lambda u_h,$$

$${\hat{\sigma }}_{d,h}\text{mathchar}58={\bar{p}}_h\mathbb{I}-\nu \nabla u_h-\frac{\nu \alpha }{h_K}\left({\bar{u}}_h-u_h\right)\otimes n,$$

$$\nabla \cdot u_h=0\forall x\in K,\forall K\in \mathcal{T},$$

$$[[u_h]]$$

$$u_h\cdot n$$

$$0={\int }_K{q}_h\nabla \cdot u_{h}dx\forall q_h\in P_{k-1}(K),\forall K\in \mathcal{T}.$$

$$0=\sum _{F\in {\mathcal{F}}_I}{\int }_F[[u_h]]{\bar{q}}_{h}ds+\sum _{F\in {\mathcal{F}}_B}{\int }_F\left(u_h-{\bar{u}}_h\right)\cdot n{\bar{q}}_{h}ds\forall {\bar{q}}_h\in {\bar{Q}}_h.$$

$$\frac{d}{dt}{\int }_K{u}_{h}dx={\int }_{K}fdx-{\int }_{\partial K}{\hat{\sigma }}_{h}nds\forall K\in \mathcal{T}.$$

$$\frac{d}{dt}{\int }_{\Omega }{u}_{h}dx={\int }_{\Omega }fdx-{\int }_{\partial \Omega }(1-\lambda )({\bar{u}}_h\cdot n){\bar{u}}_{h}ds-{\int }_{\partial \Omega }hds.$$

$$\frac{d}{dt}{\int }_K{u}_h\cdot e_{j}dx+{\int }_{\partial K}\left({\hat{\sigma }}_h\cdot n\right)\cdot e_{j}ds={\int }_{K}f\cdot e_{j}dx,$$

$$\frac{d}{dt}\sum _K{\int }_K{\left|u_h\right|}^{2}dx\le 0.$$

$$\begin{gathered} \sum _K\frac{1}{2}{\int }_K{\partial }_t{\left|u_h\right|}^{2}dx+\sum _K\frac{1}{2}{\int }_{\partial K}\left(u_h\cdot n\right){\left|u_h\right|}^{2}ds \\ -\sum _K\frac{1}{2}{\int }_{\partial K}\left(u_h\cdot n\right){\left|{\bar{u}}_h\right|}^{2}ds+\sum _K\frac{1}{2}{\int }_{\partial K}\left|u_h\cdot n\right|{\left|u_h-{\bar{u}}_h\right|}^{2}ds \\ +\sum _K{\int }_K\nu {\left|\nabla u_h\right|}^{2}dx+\sum _K{\int }_{\partial K}\frac{\nu \alpha }{h_K}{\left|{\bar{u}}_h-u_h\right|}^{2}ds \\ +2\sum _K{\int }_{\partial K}\nu \left(\nabla u_h\cdot n\right)\cdot \left({\bar{u}}_h-u_h\right)ds+{\int }_{{\Gamma }_N}(1-\lambda )({\bar{u}}_h\cdot n){\left|{\bar{u}}_h\right|}^{2}ds \\ -\sum _K{\int }_K\left(u_h\otimes u_h\right)\text{mathchar}58\nabla u_{h}dx=0, \end{gathered}$$

$$-\sum _K\frac{1}{2}{\int }_{\partial K}\left(u_h\cdot n\right){\left|{\bar{u}}_h\right|}^{2}ds=-\frac{1}{2}{\int }_{{\Gamma }_N}\left({\bar{u}}_h\cdot n\right){\left|{\bar{u}}_h\right|}^{2}ds.$$

$$-\sum _K{\int }_K\left(u_h\otimes u_h\right)\text{mathchar}58\nabla u_{h}dx=-\frac{1}{2}\sum _K{\int }_{\partial K}\left(u_h\cdot n\right){\left|u_h\right|}^{2}ds.$$

$$\begin{gathered} \frac{1}{2}\sum _K{\int }_K{\partial }_t{\left|u_h\right|}^{2}dx=-\frac{1}{2}\sum _K{\int }_{\partial K}\left|u_h\cdot n\right|{\left|u_h-{\bar{u}}_h\right|}^{2}ds \\ -\sum _K{\int }_K\nu {\left|\nabla u_h\right|}^{2}dx-\sum _K{\int }_{\partial K}\frac{\nu \alpha }{h_K}{\left|{\bar{u}}_h-u_h\right|}^{2}ds \\ -2\sum _K{\int }_{\partial K}\nu \left(\nabla u_{h}n\right)\cdot \left({\bar{u}}_h-u_h\right)ds-\frac{1}{2}{\int }_{{\Gamma }_N}\left|{\bar{u}}_h\cdot n\right|{\left|{\bar{u}}_h\right|}^{2}ds, \end{gathered}$$

$${\int }_{{\Gamma }_N}(1-\lambda )({\bar{u}}_h\cdot n){\left|{\bar{u}}_h\right|}^{2}ds-\frac{1}{2}{\int }_{{\Gamma }_N}\left({\bar{u}}_h\cdot n\right){\left|{\bar{u}}_h\right|}^{2}ds=\frac{1}{2}{\int }_{{\Gamma }_N}\left|{\bar{u}}_h\cdot n\right|{\left|{\bar{u}}_h\right|}^{2}ds.$$

$$\begin{gathered} \sum _K{\int }_K\nu {\left|\nabla u_h\right|}^{2}dx+\sum _K{\int }_{\partial K}\frac{\nu \alpha }{h_K}{\left|{\bar{u}}_h-u_h\right|}^{2}ds \\ \ge 2\left|\sum _K{\int }_{\partial K}\nu \left(\nabla u_h\cdot n\right)\cdot \left({\bar{u}}_h-u_h\right)ds\right|, \end{gathered}$$

$${\sigma }_h^{n+\theta }={\sigma }_{a,h}^{n+\theta }+{\sigma }_{d,h}^{n+\theta }\text{mbox}where{\sigma }_{a,h}^{n+\theta }=u_h^{n+\theta }\otimes u_h^n,$$

$${\hat{\sigma }}_h^{n+\theta }={\hat{\sigma }}_{a,h}^{n+\theta }+{\hat{\sigma }}_{d,h}^{n+\theta }\text{mbox}where{\hat{\sigma }}_{a,h}^{n+\theta }={\sigma }_{a,h}^{n+\theta }+({\bar{u}}_h^{n+\theta }-u_h^{n+\theta })\otimes \lambda u_h^n.$$

$$0=$$

Full text

††journal:

1 Introduction

Numerous finite element methods for the incompressible Navier–Stokes equations result in approximate velocity fields that are not pointwise divergence-free. This lack of pointwise satisfaction of the continuity equation typically leads to violation of conservation laws beyond just mass conservation, such as conservation of energy. A key issue is that, in the absence of a pointwise solenoidal velocity field, the conservative and advective format of the Navier–Stokes equations are not equivalent. The review paper by John et al (2017) presents cases for the Stokes limit where the lack of pointwise enforcement of the continuity equation can lead to large solution errors. Elements that are stable (in sense of the inf-sup condition), but do not enforce the continuity equation pointwise, such as the Taylor–Hood, Crouzeix–Raviart, and MINI elements, can suffer from large errors in the pressure, which in turn can pollute the velocity approximation. The concept of ‘pressure-robustness’ to explain the aforementioned issues is discussed by John et al (2017). A second issue is when a computed velocity field that is not pointwise divergence-free is used as the advective velocity in a transport solver. The lack of pointwise incompressibility can lead to spurious results and can compromise stability of the transport equation.

Discontinuous Galerkin (DG) finite element methods provide a natural framework for handling the advective term in the Navier–Stokes equations, and have been studied extensively in this context, e.g. (Bassi et al, 2006; Cockburn et al, 2004; Di Pietro and Ern, 2012; Ferrer and Willden, 2011; Rhebergen et al, 2013; Shahbazi et al, 2007). A difficulty in the construction of DG methods for the Navier–Stokes equations is that it is not possible to have both an energy-stable and locally momentum conserving method unless the approximate velocity is exactly divergence-free (Cockburn et al, 2004, p. 1068). To overcome this problem, a post-processing operator was introduced by Cockburn et al (2004). The operator, which is a slight modification of the Brezzi–Douglas–Marini interpolation operator (see e.g. (Boffi et al, 2013)), applied to the DG approximate velocity field generates a post-processed velocity that is pointwise divergence-free. Key to the operator is that it can be applied element-wise and is therefore inexpensive to apply. A second issue with DG methods, and a common criticism, is that the number of degrees-of-freedom on a given mesh is considerably larger than for a conforming method. This is especially the case in three spatial dimensions.

An approach to representing pointwise divergence-free velocity fields is to use a $H({\rm div})$-conforming velocity field, in which the normal component of the velocity is continuous across facets, together with a discontinuous pressure field from an appropriate space. Such a velocity space can be constructed by using a $H({\rm div})$-conforming finite element space, or by enforcing the desired continuity via hybridization (Boffi et al, 2013). However, construction of $H({\rm div})$-conforming methods for the Navier–Stokes (and Stokes) equations is not straightforward as the tangential components of the viscous stress on cell facets must be appropriately handled. Moreover, for advection dominated flows it is not immediately clear how the advective terms can be appropriately stabilized. Examples of hybridization for the Stokes equations can be found in (Carrero et al, 2005; Cockburn and Gopalakrishnan, 2005a, b), and for the Navier–Stokes equations in (Lehrenfeld and Schöberl, 2016).

A synthesis of discontinuous Galerkin and hybridized methods has lead to the development of hybridizable Discontinuous Galerkin (HDG) finite element methods (Cockburn et al, 2009; Labeur and Wells, 2007). These methods were introduced with the purpose of reducing the computational cost of DG methods on a given mesh, while retaining the attractive conservation and stability properties of DG methods. This is achieved as follows. The governing equations are posed cell-wise in terms of the approximate fields on a cell and numerical fluxes, in which the latter depends on traces of the approximate fields and fields that are defined only on facets. Fields defined on a cell are not coupled directly to fields on neighboring cells, but ‘communicate’ only via the fields that are defined on facets. By coupling degrees of freedom on a cell only to degrees of freedom of the facet functions, cell degrees of freedom can be eliminated in favor of facet degrees of freedom only. The result is that the HDG global system of algebraic equations is significantly smaller than those obtained using DG.

It has been shown that, after post-processing, solutions obtained by HDG methods may show super-convergence results for elliptic problems (for polynomial approximations of order $k$, the order of accuracy is order $k+2$ in the $L^{2}$-norm). This property has been exploited also in the context of the Navier–Stokes equations by, e.g., (Cesmelioglu et al, 2016; Nguyen et al, 2011). Although the velocity field is not automatically pointwise divergence-free, a post-processing is applied that results in an approximate velocity field that is exactly divergence-free and $H({\rm div})$-conforming and super-converges for low Reynolds number flows. Super-convergence is, however, lost when the flow is convection dominated.

We use the HDG approach to construct a simple discretization of the Navier–Stokes equations in which the computed velocity field is $H({\rm div})$-conforming and pointwise divergence-free. To achieve this, we first note that unlike many other HDG methods for incompressible flows (Cesmelioglu et al, 2016; Cockburn and Gopalakrishnan, 2009; Cockburn et al, 2011; Lehrenfeld and Schöberl, 2016; Nguyen et al, 2010, 2011; Qiu and Shi, 2016), the HDG methods of Labeur and Wells (2012) and Rhebergen and Cockburn (2012) involve facet unknowns for the pressure. The pressure field on a cell plays the role of cell-wise Lagrange multiplier to enforce the continuity equation, whereas the facet pressure unknowns play the role of Lagrange multipliers enforcing continuity of the normal component of the velocity across cell boundaries (Rhebergen and Wells, 2017). It was shown already in (Labeur and Wells, 2012) that if the polynomial approximation of the element pressure on simplices is one order lower than the polynomial approximation of the velocity that the approximate velocity field is exactly divergence-free on cells. However, the method in (Labeur and Wells, 2012) could not simultaneously satisfy mass conservation, momentum conservation and energy stability. This shortcoming is due to the computed velocity field for the method in (Labeur and Wells, 2012) not being $H({\rm div})$-conforming. We note that fast solvers for the Stokes part of the problem are developed and analysed in (Rhebergen and Wells, 2018).

In this paper we show that if the facet pressure space is chosen appropriately, we obtain approximate velocity fields that are $H({\rm div})$-conforming and pointwise divergence-free. We are guided in this by the stability analysis in (Rhebergen and Wells, 2017) for the Stokes problem, which provides guidance on the permissible function spaces. The consequences of this modification of the method of (Labeur and Wells, 2012) are profound: the method proposed in this work results in a scheme that is both mass and momentum conserving (locally and globally), energy stable and pressure-robust. We summarize properties of the proposed method and those of (Labeur and Wells, 2012) in table 1.

The remainder of this paper is organized as follows. Section 2 briefly introduces the Navier–Stokes problem, which is followed by the main result of this paper in section 3; a momentum conserving and energy stable HDG method for the Navier–Stokes equations with pointwise solenoidal and $H({\rm div})$-conforming velocity field. Numerical results are presented in section 4 and conclusions are drawn in section 5.

2 Incompressible Navier–Stokes problem

Let $\Omega\subset\mathbb{R}^{d}$ be a polygonal ($d=2$) or polyhedral ($d=3$) domain with boundary outward unit normal $n$, and let the time interval of interest be given by $I=(0,t_{N}]$. Given the kinematic viscosity $\nu\in\mathbb{R}^{+}$ and forcing term $f\mathrel{\mathop{\mathchar 58\relax}}\Omega\times I\to\mathbb{R}^{d}$, the Navier–Stokes equations for the velocity field $u\mathrel{\mathop{\mathchar 58\relax}}\Omega\times I\to\mathbb{R}^{d}$ and kinematic pressure field $p\mathrel{\mathop{\mathchar 58\relax}}\Omega\times I\to\mathbb{R}$ are given by

[TABLE]

where $\sigma$ is the momentum flux:

[TABLE]

and $\mathbb{I}$ is the identity tensor and $(a\otimes b)_{ij}=a_{i}b_{j}$.

We partition the boundary of $\Omega$ such that $\partial\Omega=\Gamma_{D}\cup\Gamma_{N}$ and $\Gamma_{D}\cap\Gamma_{N}=\emptyset$. Given $h\mathrel{\mathop{\mathchar 58\relax}}\Gamma_{N}\times I\to\mathbb{R}^{d}$ and a solenoidal initial velocity field $u_{0}\mathrel{\mathop{\mathchar 58\relax}}\Omega\to\mathbb{R}^{d}$, we prescribe the following boundary and initial conditions:

[TABLE]

On inflow parts of $\Gamma_{N}$ ($u\cdot n<0$) we impose the total momentum flux, i.e., $\sigma\cdot n=h$. On outflow parts of $\Gamma_{N}$ ($u\cdot n\geq 0$), only the diffusive part of the momentum flux is prescribed, i.e., $\sigma_{d}\cdot n=h$.

Equation 1a is the conservative form of the Navier–Stokes equation. With satisfaction of the incompressibility constraint, eq. 1b, the momentum equation (1a) can be equivalently expressed as:

[TABLE]

where $\chi\in[0,1]$. For numerous finite element methods, the approximate velocity field is not pointwise or locally (in a weak sense) solenoidal. In such cases, it can be shown that momentum is conserved if $\chi=1$, while energy stability can be proven if $\chi=1/2$. For stabilized finite element methods in which the continuity equation is not satisfied locally, manipulations of the advective term can be applied to achieve momentum conservation (Hughes and Wells, 2005).

The mass conserving (mixed-order) hybridizable discontinuous Galerkin method of Labeur and Wells (2012) is based on a weak formulation of eq. 4. It was proven to be locally momentum conserving for $\chi=1$ and energy stable for $\chi=1/2$, but in their analysis both properties could not be satisfied simultaneously. We will prove how the method can be formulated such that mass and momentum conservation, and energy stability can be satisfied simultaneously, and the method be made invariant with respect to $\chi$.

3 A hybridizable discontinuous Galerkin method

We present a hybridizable discontinuous Galerkin method for the Navier–Stokes problem for which the approximate velocity field is pointwise divergence-free.

3.1 Preliminaries

Let $\mathcal{T}\mathrel{\mathop{\mathchar 58\relax}}=\mathinner{\left\{K\right\}}$ be a triangulation of the domain $\Omega$ into non-overlapping simplex cells $K$. The boundary of a cell is denoted by $\partial K$ and the outward unit normal vector on $\partial K$ by $n$. Two adjacent cells $K^{+}$ and $K^{-}$ share an interior facet $F\mathrel{\mathop{\mathchar 58\relax}}=\partial K^{+}\cap\partial K^{-}$. A facet of $\partial K$ that lies on the boundary of the domain $\partial\Omega$ is called a boundary facet. The sets of interior and boundary facets are denoted by $\mathcal{F}_{I}$ and $\mathcal{F}_{B}$, respectively. The set of all facets is denoted by $\mathcal{F}\mathrel{\mathop{\mathchar 58\relax}}=\mathcal{F}_{I}\cup\mathcal{F}_{B}$.

3.2 Semi-discrete formulation

Consider the following finite element spaces:

[TABLE]

where $P_{l}(D)$ denotes the space of polynomials of degree $l>0$ on a domain $D$. Note that the spaces $V_{h}$ and $Q_{h}$ are defined on the whole domain $\mathcal{T}$, whereas the spaces $\bar{V}_{h}$ and $\bar{Q}_{h}$ are defined only on facets of the triangulation.

The spaces $V_{h}$ and $Q_{h}$ are discontinuous across cell boundaries, hence the trace of a function $a\in V_{h}$ may be double-valued on cell boundaries. At an interior facet, $F$, we denote the traces of $a\in V_{h}$ by $a^{+}$ and $a^{-}$. We introduce the jump operator $\llbracket a\rrbracket\mathrel{\mathop{\mathchar 58\relax}}=a^{+}\cdot n^{+}+a^{-}\cdot n^{-}$, where $n^{\pm}$ the outward unit normal on $\partial K^{\pm}$.

We now state the weak formulation of the proposed method: given a forcing term $f\in\mathinner{\left[L^{2}(\Omega)\right]}^{d}$, boundary condition $h\in\mathinner{\left[L^{2}(\Gamma_{N})\right]}^{d}$ and viscosity $\nu$, find $u_{h},\bar{u}_{h},p_{h},\bar{p}_{h}\in V_{h}\times\bar{V}_{h}\times Q_{h}\times\bar{Q}_{h}$ such that:

[TABLE]

where $\hat{\sigma}_{h}\mathrel{\mathop{\mathchar 58\relax}}=\hat{\sigma}_{a,h}+\hat{\sigma}_{d,h}$ is the ‘numerical flux’ on cell facets. The advective part of the numerical flux is given by:

[TABLE]

where $\lambda$ is an indicator function that takes on a value of unity on inflow cell boundaries (where $u_{h}\cdot n<0$) and a value of zero on outflow cell facets (where $u_{h}\cdot n\geq 0$). This definition of the numerical flux provides upwinding of the advective component of the flux. The diffusive part of the numerical flux is defined as

[TABLE]

where $\alpha>0$ is a penalty parameter as is typical of Nitsche and interior penalty methods. It is proven in (Wells, 2011; Rhebergen and Wells, 2017) that $\alpha$ needs to be sufficiently large to ensure stability.

A key feature of this formulation, and what distinguishes it from standard discontinuous Galerkin methods, is that functions on cells (functions in $V_{h}$ and $Q_{h}$) are not coupled across facets directly via the numerical flux. Rather, fields on neighboring cells are coupled via the facet functions $\bar{u}_{h}$ and $\bar{p}_{h}$. The fields $u_{h}$ and $p_{h}$ can therefore be eliminated locally via static condensation, resulting in a global system of equations in terms of the facet functions only. This substantially reduces the size of the global systems compared to a standard discontinuous Galerkin method on the same mesh, yet still permits the natural incorporation of upwinding and cell-wise balances.

The weak formulation presented here is the weak formulation of Labeur and Wells (2012) with conservative form of the advection term ($\chi=1$ in eq. 4). The key difference is that we have been more prescriptive on the relationships between the finite element spaces in eq. 5, and we will prove that this leads to some appealing properties. In particular, the spaces in eq. 5 are such that: for $u_{h}\in\mathinner{\left[P_{k}(K)\right]}^{d}$, $\nabla\cdot u_{h}\in P_{k-1}(K)$ and $u_{h}\cdot n\in P_{k}(F)$; and for $\bar{u}_{h}\in\mathinner{\left[P_{k}(F)\right]}^{d}$, $\bar{u}_{h}\cdot n\in P_{k}(F)$. Furthermore, the function spaces have been chosen such that the resulting method is inf-sup stable, see (Rhebergen and Wells, 2017). The resulting weak formulation can be shown to be equivalent to a weak formulation in which the approximate velocity field lies in the Brezzi–Douglas–Marini (BDM) finite element space (Rhebergen and Wells, 2017, Section 3.4). Hybridization of other $H({\rm div})$ conforming finite element spaces, see e.g. (Boffi et al, 2013), are also possible.

Proposition 1 (mass conservation)

If $u_{h}\in V_{h}$ and $\bar{u}_{h}\in\bar{V}_{h}$ satisfy eq. 6, with $V_{h}$ and $\bar{V}_{h}$ defined in eq. 5, then

[TABLE]

and

[TABLE]

Proof

Applying integration-by-parts to eq. 6a:

[TABLE]

Since $q_{h}$, $\nabla\cdot u_{h}\in P_{k-1}(K)$, pointwise satisfaction of the continuity equation, eq. 9, follows.

It follows from eq. 6b that:

[TABLE]

Since $\bar{q}_{h}$, $u_{h}\cdot n$, $\bar{u}_{h}\cdot n\in P_{k}(F)$, eq. 10 follows. ∎

Proposition 1 is a stronger statement of mass conservation than in Labeur and Wells (2012, Proposition 4.2), in which mass conservation for the mixed-order case was proved locally (cell-wise) in an integral sense only. Under certain conditions, implementations in (Labeur and Wells, 2012) satisfy eq. 9, but not eq. 10. We will show that this difference is critical for the formulation in this work as it allows simultaneous satisfaction of momentum conservation and energy stability.

We next show momentum conservation for the semi-discrete weak formulation in terms of the numerical flux.

Proposition 2 (momentum conservation)

Let $u_{h},\bar{u}_{h},p_{h},\bar{p}_{h}\in V_{h}\times\bar{V}_{h}\times Q_{h}\times\bar{Q}_{h}$ satisfy eq. 6. Then,

[TABLE]

Furthermore, if $\Gamma_{D}=\emptyset$,

[TABLE]

Proof

In eq. 6c, set $v_{h}=e_{j}$ on $K$, where $e_{j}$ is a canonical unit basis vector, and set $v_{h}=0$ on $\mathcal{T}\backslash K$ in eq. 6c:

[TABLE]

which proves eq. 13. Equation 14 follows immediately by setting $v_{h}=e_{j}$ in eq. 6c, $\bar{v}_{h}=-e_{j}$ in eq. 6d and summing the two results. ∎

We next prove that the method is also globally energy stable.

Proposition 3 (global energy stability)

If $u_{h},\bar{u}_{h},p_{h},\bar{p}_{h}\in V_{h}\times\bar{V}_{h}\times Q_{h}\times\bar{Q}_{h}$ satisfy eq. 6, for homogeneous boundary conditions, $f=0$ and for a suitably large $\alpha$:

[TABLE]

Proof

Setting $q_{h}=-p_{h}$, $\bar{q}_{h}=-\bar{p}_{h}$, $v_{h}=u_{h}$ and $\bar{v}_{h}=-\bar{u}_{h}$ in eqs. 6a, 6b, 6c and 6d and inserting the expressions for the numerical fluxes (eqs. 2, 7 and 8), and summing:

[TABLE]

where we have used that $\lambda u_{h}\cdot n=\mathinner{\left(u_{h}\cdot n-|u_{h}\cdot n|\right)}/2$, and applied integration-by-parts to the pressure gradient terms. Since $\bar{u}_{h}$ is single-valued on facets, the normal component of $u_{h}$ is continuous across facets and $\bar{u}_{h}\cdot n=u_{h}\cdot n$ on the domain boundary (see proposition 1), the third integral on the left-hand side of eq. 17 can be simplified:

[TABLE]

We consider now the last term on the left-hand side of eq. 17. On each cell $K$ it holds that $-u_{h}\otimes u_{h}\mathrel{\mathop{\mathchar 58\relax}}\nabla u_{h}=(\nabla\cdot u_{h})(u_{h}\cdot u_{h})/2-\nabla\cdot((u_{h}\otimes u_{h})\cdot u_{h})/2=-\nabla\cdot((u_{h}\otimes u_{h})\cdot u_{h})/2$, since $\nabla\cdot u_{h}=0$ (by proposition 1). It follows that

[TABLE]

Combining eqs. 17, 18 and 19,

[TABLE]

where we have used that

[TABLE]

It can be proven that there exists an $\alpha>0$, independent of $h_{K}$, such that

[TABLE]

(see (Wells, 2011, Lemma 5.2) and (Rhebergen and Wells, 2017, Lemma 4.2)). Therefore, the right-hand side of eq. 20 is non-positive, proving eq. 16. ∎

The key results that enable us to prove global energy stability for this conservative form of the Navier–Stokes equations are: (a) the pointwise solenoidal velocity field; and (b) continuity of the normal component of the velocity field across facets. The latter point is not fulfilled by the method in (Labeur and Wells, 2012).

3.3 A fully-discrete weak formulation

We now consider a fully-discrete formulation. We partition the time interval $I$ into an ordered series of time levels $0=t^{0}<t^{1}<\cdots<t^{N}$. The difference between each time level is denoted by $\Delta t^{n}=t^{n+1}-t^{n}$. To discretize in time, we consider the $\theta$-method and denote midpoint values of a function $y$ by $y^{n+\theta}\mathrel{\mathop{\mathchar 58\relax}}=(1-\theta)y^{n}+\theta y^{n+1}$. Following Labeur and Wells (2012), the convective velocity will be evaluated at the current time $t^{n}$, thereby linearizing the problem, i.e.:

[TABLE]

and

[TABLE]

The time-discrete counterpart of eq. 6 is: given $u_{h}^{n},\bar{u}_{h}^{n},p_{h}^{n},\bar{p}_{h}^{n}\in V_{h}\times\bar{V}_{h}\times Q_{h}\times\bar{Q}_{h}$ at time $t^{n}$, the forcing term $f^{n+\theta}\in\mathinner{\left[L^{2}(\Omega)\right]}^{d}$, the boundary condition $h^{n+\theta}\in\mathinner{\left[L^{2}(\Gamma_{N})\right]}^{d}$, and the viscosity $\nu$, find $u_{h}^{n+1},\bar{u}_{h}^{n+1},p_{h}^{n+1},\bar{p}_{h}^{n+1}\in V_{h}\times\bar{V}_{h}\times Q_{h}\times\bar{Q}_{h}$ such that mass conservation,

[TABLE]

are satisfied for all $v_{h},\bar{v}_{h},q_{h},\bar{q}_{h}\in V_{h}\times\bar{V}_{h}\times Q_{h}\times\bar{Q}_{h}$. Here $\lambda$ is evaluated using the known velocity field at time $t^{n}$.

In section 3.2 we proved that the semi-discrete formulation eq. 6 is momentum conserving, energy stable and exactly mass conserving when using the function spaces given by eq. 5. We show next that the fully-discrete formulation given by eq. 25 inherits these properties.

Proposition 4 (fully-discrete mass conservation)

If $u_{h}^{n+1}\in V_{h}$ and $\bar{u}_{h}^{n+1}\in\bar{V}_{h}$ satisfy eq. 25, then

[TABLE]

and

[TABLE]

Proof

The proof is similar to that of proposition 1 and therefore omitted. ∎

Proposition 5 (fully-discrete momentum conservation)

If $u_{h}^{n},\bar{u}_{h}^{n},p_{h}^{n},\bar{p}_{h}^{n}\in V_{h}\times\bar{V}_{h}\times Q_{h}\times\bar{Q}_{h}$ and $u_{h}^{n+1},\bar{u}_{h}^{n+1},p_{h}^{n+1},\bar{p}_{h}^{n+1}\in V_{h}\times\bar{V}_{h}\times Q_{h}\times\bar{Q}_{h}$ satisfy eq. 25, then

[TABLE]

Furthermore, if $\Gamma_{D}=\emptyset$,

[TABLE]

Proof

The proof is similar to that of proposition 2 and therefore omitted. ∎

Proposition 6 (fully-discrete energy stability)

If $u_{h}^{n},\bar{u}_{h}^{n},p_{h}^{n},\bar{p}_{h}^{n}\in V_{h}\times\bar{V}_{h}\times Q_{h}\times\bar{Q}_{h}$ and $u_{h}^{n+1},\bar{u}_{h}^{n+1},p_{h}^{n+1},\bar{p}_{h}^{n+1}\in V_{h}\times\bar{V}_{h}\times Q_{h}\times\bar{Q}_{h}$ satisfy eq. 25, then with homogeneous boundary conditions, no forcing terms, for suitably large $\alpha$, and $\theta\geq 1/2$,

[TABLE]

Proof

Setting $q_{h}=-\theta p_{h}^{n+\theta}$, $\bar{q}_{h}=-\theta\bar{p}_{h}^{n+\theta}$, $v_{h}=u_{h}^{n+\theta}$ and $\bar{v}_{h}=-\bar{u}_{h}^{n+\theta}$, in eqs. 25a, 25b, 25c and 25d, adding the results, using the expressions for the diffusive fluxes, given by eqs. 2 and 8, partial integration of the pressure gradient terms and using that $\nabla\cdot u_{h}^{n}=0$ by proposition 4, we obtain, using the same steps as in the proof of proposition 3,

[TABLE]

The first term on the left-hand side of eq. 31 can be reformulated as

[TABLE]

Inserting this expression into eq. 31:

[TABLE]

As in proposition 3, there exists an $\alpha>0$, independent of $h_{K}$, such that the right hand side of eq. 33 is non-positive. The result follows. ∎

4 Numerical examples

We now demonstrate the performance of the method for a selection of numerical examples, paying close attention to mass and momentum conservation, and energy stability.

For all stationary examples considered, exact solutions are known. For the stationary examples we use a fixed-point iteration with stopping criterion $|e_{p}^{i+1}-e_{p}^{i}|/(e_{p}^{i+1}+e_{p}^{i})\leq{\rm TOL}$, where $e_{p}^{i}$ is the pressure error in the $L^{2}$ norm at the $i$th iterate, and ${\rm TOL}$ is a given tolerance that we set to $10^{-4}$. All unsteady examples use $\theta=1$. In all examples we set the penalty parameter to be $\alpha=6k^{2}$.

In the implementation we apply cell-wise static condensation such that only the degrees-of-freedom associated with the facet spaces appear in the global system. Compared to standard discontinuous Galerkin methods, this significantly reduces the size of the global system. We could eliminate the facet pressure field and use a BDM element, see (Lehrenfeld and Schöberl, 2016), and the BDM normal velocity in place of $\bar{u}_{h}\cdot n$. However, we feel that handling all fields in a hybridized framework offers some simplicity.

Examples have been implemented using the NGSolve finite element library (Schöberl, 2014). All examples use unstructured simplicial meshes.

4.1 Kovasznay flow

We consider the steady, two-dimensional analytical solution of the Navier–Stokes equations from Kovasznay (Kovasznay, 1948) on a domain $\Omega=\mathinner{\left(-0.5,1\right)}\times\mathinner{\left(-0.5,1.5\right)}$. For a Reynolds number $\mathit{Re}$, let the viscosity be given by $\nu=1/\mathit{Re}$. The solution to the Kovasznay problem is:

[TABLE]

where $C$ is an arbitrary constant, and where

[TABLE]

We choose $C$ such that the mean pressure on $\Omega$ is zero. The Kovasznay flow solution in eq. 34 is used to set Dirichlet boundary conditions for the velocity on $\partial\Omega$.

The $L^{2}$-error and rates of convergence are presented in table 2 for $\mathit{Re}=40$ using a series of refined meshes. Optimal rates of convergence are observed for both the velocity field (order $k+1$) and pressure field (order $k$). The divergence of the approximate velocity field is of machine precision in all cases.

4.2 Position-dependent Coriolis force

We now consider the test case from (Linke and Merdon, 2016a, Section 3.2). In particular, we consider on the unit square $(0,1)\times(0,1)$ the steady Navier–Stokes equations augmented with a position-dependent Coriolis force: $\nabla\cdot\sigma+2C\times u=0$ and $\nabla\cdot u=0$, where we set $2C\times u=-2x_{2}(-u_{2},u_{1})$. On boundaries we set $u=(1,0)$. The exact solution to this problem is given by $p=x_{2}^{2}-1/3$ and $u=(1,0)$.

It was shown in (Linke and Merdon, 2016a) that the Scott–Vogelius finite element, in which the velocity is approximated in divergence-free function spaces, is able to produce the exact velocity field while the velocity computed using a Taylor–Hood finite element method is polluted by the pressure error, in part due to the approximate velocity field not being exactly divergence-free. Furthermore, it is shown in (Linke and Merdon, 2016a) that as $\nu\to 0$, the velocity error increases for the Taylor–Hood finite element method.

In table 3 we show the results obtained using the HDG method presented in section 3 for $k=2$. It shows the computed error in the $L^{2}$ norm for the velocity, pressure and divergence errors. Errors in the velocity and velocity divergence are of machine precision, regardless of $\nu$. The HDG method therefore obtains the same quality of solution as produced using the Scott–Vogelius finite element in (Linke and Merdon, 2016a). We do not consider the $k=3$ case because for this discretization the pressure is approximated by quadratic polynomials and so the pressure error is also of machine precision.

4.3 Pressure-robustness

We next demonstrate that our method is pressure-robust and compare the results with those obtained using the method of (Labeur and Wells, 2012). For this we use a test case proposed in (Lederer et al, 2017, Section 6.1). On the unit square $(0,1)\times(0,1)$ we consider the steady Navier–Stokes equations where the boundary conditions and source terms are such that the exact solution is given by $u=\mathrm{curl}\zeta$, with $\zeta=x_{1}^{2}(x_{1}-1)^{2}x_{2}^{2}(x_{2}-1)^{2}$ and $p=x_{1}^{7}+x_{2}^{7}-1/4$. We choose $k=3$ and vary the viscosity $\nu$.

For a mixed velocity-pressure approximation, it can be proven that the method of (Labeur and Wells, 2012) results in an approximate velocity field that is pointwise divergence free, but not $H({\rm div})$-conforming. As such, the method of (Labeur and Wells, 2012) cannot be shown to be pressure-robust. This is confirmed by the results presented in table 4. For the method of (Labeur and Wells, 2012) it is observed that the error in the velocity field depends on $\nu^{-1}\mathinner{\!\left\lVert p_{h}-p\right\rVert}$, while the proposed method is pressure robust. The velocity error does not change with viscosity in table 4.

4.4 Transient higher-order potential flow

In this test, taken from (Linke and Merdon, 2016b, Section 6.6), we solve the time dependent Navier–Stokes equations eq. 1 on the domain $\Omega=[-1,1]^{2}$. This test case studies the time-dependent exact velocity $u(t)=\min(t,1)\nabla\chi$ where $\chi$ is a smooth harmonic potential given by $\chi=x_{1}^{3}x_{2}-x_{2}^{3}x_{1}$. The pressure gradient then satisfies $\nabla p=-\nabla\mathinner{\!\left\lvert u\right\rvert}^{2}/2-\partial_{t}\mathinner{\left(\min(t,1)\nabla\chi\right)}$. We impose the exact velocity solution as Dirichlet boundary condition on all of $\partial\Omega$.

For the simulations we used a grid with 2048 cells, set the time step equal to $\Delta t=0.01$ and compute the solution on the time interval $[0,2]$. Figure 1 shows the velocity and pressure errors as a function of time. We used both $k=2$ and $k=3$, and consider $\nu=1/500$ and $\nu=1/2000$. We observe that the error in pressure and velocity is more or less the same regardless of $\nu$.

Over the computational time interval, using $k=2$ or $k=3$ on a mesh with 2048 cells, for either $\nu=1/500$ and $\nu=1/2000$, the $L^{2}$-norm of the divergence reaches $1.4\times 10^{-10}$ in one point but is otherwise always of the order $10^{-11}$. The momentum balance, in absolute value, never exceeds $3.4\times 10^{-12}$.

4.5 Two-dimensional flow past a circular obstacle

In this test case we consider flow past a circular obstacle (see e.g. (Lehrenfeld and Schöberl, 2016; Schäfer et al, 1996)). The domain is a rectangular channel, $[0,2.2]\times[0,0.41]$, with a circular obstacle of radius $r=0.05$ centered at $(0.2,0.2)$. On the inflow boundary ($x_{1}=0$) we prescribe the $x_{1}$-component of the velocity to be $u_{1}=6x_{2}(0.41-x_{2})/0.41^{2}$. The $x_{2}$-component of the velocity is prescribed as $u_{2}=0$. Homogeneous Dirichlet boundary conditions are applied on the walls ($x_{2}=0$ and $x_{2}=0.41$), and on the obstacle. On the outflow boundary ($x_{1}=2.2$) we prescribe $\sigma_{d}\cdot n=0$. The viscosity is set as $\nu=10^{-3}$. We choose $k=3$ and set $\Delta t=5\times 10^{-5}$ so that the spatial discretization error dominates the temporal discretization error. For the initial condition, we impose the steady Stokes solution of this problem. The mesh of the domain has 6784 cells and we consider the time interval $[0,5]$.

At each time step we compute the drag and lift coefficients, which are defined as

[TABLE]

where $e_{1}$ and $e_{2}$ are unit vectors in the $x_{1}$ and $x_{2}$ directions, respectively, and $\Gamma_{C}$ is the surface of the circular object. We compute a maximum drag coefficient of $C_{D}=3.23232$ and minimum drag coefficient of $C_{D}=3.16583$. The maximum and minimum lift coefficients we compute are, respectively, $C_{L}=0.98251$ and $C_{L}=-1.02246$. These are comparable to those found in literature (Lehrenfeld and Schöberl, 2016; Schäfer et al, 1996). The velocity magnitude at $t=5$ is shown fig. 2.

4.6 Three-dimensional flow past a cylinder

In this test case we consider three-dimensional flow past a cylinder (see e.g. (Lehrenfeld and Schöberl, 2016; Schäfer et al, 1996)) with a time dependent inflow velocity. The domain is a cuboid shaped channel $[0,2.5]\times[0,0.41]\times[0,0.41]$ with a cylinder of radius $r_{\rm cyl}=0.05$ around the $x_{3}$-axis centered at $(x_{1},x_{2})=(0.5,0.2)$. On the inflow boundary ($x_{1}=0$) we prescribe the $x_{1}$-component of the velocity to be $u_{1}=36\sin(\pi t/8)x_{2}x_{3}(0.41-x_{2})(0.41-x_{3})/0.41^{4}$. The $x_{2}$- and $x_{3}$-components of the velocity are prescribed as $u_{2}=0$ and $u_{3}=0$. We impose homogeneous Dirichlet boundary conditions on the walls ($x_{2}=0$, $x_{2}=0.41$, $x_{3}=0$ and $x_{3}=0.41$) and on the cylinder. On the outflow boundary ($x_{1}=2.5$) we prescribe $\sigma_{d}\cdot n=0$. The viscosity is set as $\nu=10^{-3}$.

We choose $k=3$ and set $\Delta t=5\times 10^{-4}$ so that the spatial discretization error dominates the temporal discretization error. The initial condition is the Stokes solution to this problem. The mesh has 4091 cells and we compute on the time interval $[0,8]$. At each time step we compute the drag and lift coefficients, defined by eq. 36, where $r=0.41r_{\rm cyl}$ and $\Gamma_{C}$ is the surface of the cylinder. We compute maximum drag and lift coefficients of $C_{D}=2.98815$ and $C_{L}=0.00348$, respectively. Compared to Schäfer et al (1996), in which the maximum drag and lift coefficients lie in the intervals $C_{D}\in[3.2000,3.3000]$ and $C_{L}\in[0.0020,0.0040]$, we slightly under-predict the drag coefficient, but the lift coefficient lies within the same interval. Figure 3 shows the velocity magnitude at $t=4$.

5 Conclusions

We have introduced a formulation of a hybridizable discontinuous Galerkin method for the incompressible Navier–Stokes equations that computes velocity fields that are pointwise divergence-free. The construction of solenoidal velocity fields does not require post-processing or the use of finite dimensional spaces of divergence-free functions. The pointwise satisfaction of the continuity equation and the continuity of the normal component of the velocity field across cell facets allows us to prove that the method conserves momentum locally (cell-wise) and is energy stable. This is in contrast with the closely related method in Labeur and Wells (2012) which when satisfying the continuity equation pointwise can satisfy local momentum conservation or global energy stability, but not both simultaneously. The analysis that we present is supported by a range of numerical examples in two and three dimensions.