## 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 -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 -conforming finite element space, or by enforcing the desired continuity via hybridization (Boffi et al, 2013). However, construction of -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 , the order of accuracy is order in the -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 -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 -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 -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 -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.

Formulation | mass | momentum | energy | pressure |
---|---|---|---|---|

conserving | conserving | stable | robust | |

Equal order (Labeur and Wells, 2012) | only in | |||

skew-symmetric | ||||

form | ||||

Mixed order (Labeur and Wells, 2012) | only in | only in | ||

divergence | skew-symmetric | |||

form | form | |||

Proposed method |

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 -conforming velocity field. Numerical results are presented in section 4 and conclusions are drawn in section 5.

## 2 Incompressible Navier–Stokes problem

Let be a polygonal () or polyhedral () domain with boundary outward unit normal , and let the time interval of interest be given by . Given the kinematic viscosity and forcing term , the Navier–Stokes equations for the velocity field and kinematic pressure field are given by equationparentequation

(1a) | |||||

(1b) |

where is the momentum flux:

(2) |

and

is the identity tensor and

.We partition the boundary of such that and . Given and a solenoidal initial velocity field , we prescribe the following boundary and initial conditions: equationparentequation

(3a) | |||||

(3b) | |||||

(3c) |

On inflow parts of () we impose the total momentum flux, i.e., . On outflow parts of (), only the diffusive part of the momentum flux is prescribed, i.e., .

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:

(4) |

where . 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 , while energy stability can be proven if . 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 and energy stable for , 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 .

## 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 be a triangulation of the domain into non-overlapping simplex cells . The boundary of a cell is denoted by

and the outward unit normal vector on

by . Two adjacent cells and share an interior facet . A facet of that lies on the boundary of the domain is called a boundary facet. The sets of interior and boundary facets are denoted by and , respectively. The set of all facets is denoted by .### 3.2 Semi-discrete formulation

Consider the following finite element spaces: equationparentequation

(5a) | ||||

(5b) | ||||

(5c) | ||||

(5d) |

where denotes the space of polynomials of degree on a domain . Note that the spaces and are defined on the whole domain , whereas the spaces and are defined only on facets of the triangulation.

The spaces and are discontinuous across cell boundaries, hence the trace of a function may be double-valued on cell boundaries. At an interior facet, , we denote the traces of by and . We introduce the jump operator , where the outward unit normal on .

We now state the weak formulation of the proposed method: given a forcing term , boundary condition and viscosity , find such that: equationparentequation

(6a) | ||||

(6b) | ||||

and | ||||

(6c) | ||||

(6d) |

where is the ‘numerical flux’ on cell facets. The advective part of the numerical flux is given by:

(7) |

where is an indicator function that takes on a value of unity on inflow cell boundaries (where ) and a value of zero on outflow cell facets (where ). 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

(8) |

where is a penalty parameter as is typical of Nitsche and interior penalty methods. It is proven in (Wells, 2011; Rhebergen and Wells, 2017) that 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 and ) are not coupled across facets directly via the numerical flux. Rather, fields on neighboring cells are coupled via the facet functions and . The fields and 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 ( 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 , and ; and for , . 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 conforming finite element spaces, see e.g. (Boffi et al, 2013), are also possible.

###### Proposition 1 (mass conservation)

###### Proof

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)

###### Proof

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

###### Proposition 3 (global energy stability)

If satisfy eq. 6, for homogeneous boundary conditions, and for a suitably large :

(16) |

###### Proof

Setting , , and in eqs. 6d, 6c, 6b and 6a and inserting the expressions for the numerical fluxes (eqs. 8, 7 and 2), and summing:

(17) |

where we have used that , and applied integration-by-parts to the pressure gradient terms. Since is single-valued on facets, the normal component of is continuous across facets and on the domain boundary (see proposition 1), the third integral on the left-hand side of eq. 17 can be simplified:

(18) |

We consider now the last term on the left-hand side of eq. 17. On each cell it holds that , since (by proposition 1). It follows that

(19) |

(20) |

where we have used that

(21) |

It can be proven that there exists an , independent of , such that

(22) |

(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 into an ordered series of time levels . The difference between each time level is denoted by . To discretize in time, we consider the -method and denote midpoint values of a function by . Following Labeur and Wells (2012), the convective velocity will be evaluated at the current time , thereby linearizing the problem, i.e.:

(23) |

and

(24) |

The time-discrete counterpart of eq. 6 is: given at time , the forcing term , the boundary condition , and the viscosity , find such that mass conservation, equationparentequation

(25a) | ||||

(25b) | ||||

and momentum conservation, | ||||

(25c) | ||||

(25d) |

are satisfied for all . Here is evaluated using the known velocity field at time .

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)

###### Proof

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

###### Proposition 5 (fully-discrete momentum conservation)

###### Proof

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

###### Proposition 6 (fully-discrete energy stability)

If and satisfy eq. 25, then with homogeneous boundary conditions, no forcing terms, for suitably large , and ,

(30) |

###### Proof

Setting , , and , in eqs. 25d, 25c, 25b and 25a, adding the results, using the expressions for the diffusive fluxes, given by eqs. 8 and 2, partial integration of the pressure gradient terms and using that by proposition 4, we obtain, using the same steps as in the proof of proposition 3,

(31) |

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

(32) |

Inserting this expression into eq. 31:

(33) |

As in proposition 3, there exists an , independent of , 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 , where is the pressure error in the norm at the th iterate, and is a given tolerance that we set to . All unsteady examples use . In all examples we set the penalty parameter to be .

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 . 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 . For a Reynolds number , let the viscosity be given by . The solution to the Kovasznay problem is: equationparentequation