A Divergence-Conforming Hybridized Discontinuous Galerkin Method for the Incompressible Reynolds Averaged Navier-Stokes Equations

We introduce a hybridized discontinuous Galerkin method for the incompressible Reynolds Averaged Navier-Stokes equations coupled with the Spalart-Allmaras one equation turbulence model. With a special choice of velocity and pressure spaces for both element and trace degrees of freedom, we arrive at a method which returns point-wise divergence-free mean velocity fields and properly balances momentum and energy. We further examine the use of different polynomial degrees and meshes to see how the order of the scalar eddy viscosity affects the convergence of the mean velocity and pressure fields, specifically for the method of manufactured solutions. As is standard with hybridized discontinuous Galerkin methods, static condensation can be employed to remove the element degrees of freedom and thus dramatically reduce the global number of degrees of freedom. Numerical results illustrate the effectiveness of the proposed methodology.

Key Words: Hybridized discontinuous Galerkin methods; Divergence conforming discretizations; RANS: Reynolds averaged Navier-Stokes; Incompressible flow; Spalart-Allmaras turbulence model; Static condensation

The discontinuous Galerkin (DG) method was originally introduced by Reed and Hill in [1] to solve the neutron transport equation, however with little analysis of its properties. A more rigorous exploration of the method and its properties were discovered and discussed in [2, 3]. Since then, the method has gained popularity and has since been extended to a multitude of other problems governed by partial differential equations (PDEs) detailed in [4] and [5]. DG methods combine the advantages of classical finite volume methods (FVM) with finite element methods (FEM). Much like FVM, DG has a natural framework for dealing with advection problems, or more generally hyperbolic problems, but in addition it has a natural extension to higher order methods, in a similar fashion to FEM. Thus DG is capable of resolving solutions with large gradients, including shocks, as well as dealing with complex geometries. Another important advantage of having a method which is higher order is that more work can be done on processor for element interiors, limiting the overall percentage of work needed for communication in a parallel framework. In fact, it has been regularly used to solve large-scale forward [6, 7] and inverse problems [8]. However, there is a major drawback associated with this class of method. In general, DG methods introduce an increased amount of degrees of freedom (DOFs) which cannot be alleviated by increasing the polynomial order. Consequently, implicit time integration schemes are often intractable and one must resort to explicit schemes hindering the range of acceptable Courant-Friedrichs-Lewy (CFL) condition numbers. This can make the method prohibitively expensive in comparison to other existing numerical methods, such as continuous Galerkin (CG) or spectral element methods.

Within the last decade, Cockburn and his collaborators have devised a method which alleviates the high cost of DG methods using hybridization. The resulting approach, referred to as the hybridized discontinuous Galerkin (HDG) method, has been applied successfully to several types of PDEs including, but not limited to, Poisson’s equation [9], advection-diffusion equations [10, 11, 12, 13], the Stokes equations [14, 15, 16, 17], and even the Euler and Navier-Stokes equations [18, 19]. With an HDG method, there are both interior DOFs that reside on the interiors of elements and trace (or skeleton) DOFs that reside on the edges of polygons in 2-D and faces of polyhedra in 3-D. The trace DOFs act as communicators across element boundaries. In fact, static condensation can be employed to write the interior DOFs in terms of the trace DOFs, allowing one to remove the interior DOFs entirely from the discrete system of equations. In this manner, the HDG method allows one to significantly reduce the number of DOFs in the global system as compared with a corresponding DG method, rendering the HDG method competitive with higher-order CG and spectral element methods. Once the trace DOFs are determined, the interior DOFs can be realized locally in an element-by-element manner. Furthermore, with the supposition of convergence of a DG method, the equivalent HDG method is guaranteed to converge [20], with this condition being sufficient but not necessary.

Recently, there has been a concerted effort in further developing HDG methods for the incompressible Navier-Stokes equations and related PDE’s, including the development of stable mixed elements and structure-preserving discretizations [21] as well as super convergent post-processing procedures [22]. There has additionally been efforts to alleviate the associated complexity of developing new HDG discretizations. Bui-Thanh has contributed to this avenue by constructing HDG methods from a Godunov approach [23] as well as from the Rankine-Hugoniot condition [24]. These types of methodologies render the trace unknowns as upwind states, while providing a parameter-free framework.

In this work, we construct a new HDG method for the incompressible Reynolds-Averaged Navier-Stokes (RANS) equations. To construct our method, we begin with the divergence-conforming HDG method recently introduced for the incompressible Navier-Stokes equations by Rhebergen and Wells in [21], and we modify the method to handle a RANS closure model. More specifically, we design an HDG method for the Spalart-Allmaras one-equation turbulence model [25], where we treat the transport equation for the eddy viscosity as an advection-diffusion equation. It is important to note that we introduce trace DOFs for the primal variables, that is, the flow field and turbulence variables, for the sake of implementational ease. This is in contrast with many alternative HDG methods where trace DOFs represent inter-element fluxes.

A defining feature of our HDG method is that it produces a point-wise divergence-free mean velocity field. It has been shown in recent works that discretizations of the incompressible Navier-Stokes equations which exactly preserve the divergence-free constraint are typically more robust than methods which only satisfy the divergence-free constraint in a weak manner [26]. For example, divergence-conforming discretizations simultaneously conserve momentum and energy [27], while discretizations which only weakly satisfy the divergence-free constraint typically conserve either momentum or energy but not both. The velocity error is also decoupled from the pressure error for divergence-conforming discretizations, a property commonly referred to as pressure robustness [28]. Finally, mass conservation is considered to be of prime importance for coupled-flow transport [29], so divergence-free discretizations are particularly attractive for such applications. This last point inspires our current work, as the Spalart-Allmaras model may be viewed as a coupled-flow transport problem wherein the eddy viscosity is a scalar transported by the fluid flow.

An outline of this paper is as follows. In the next section, we recall the strong form of the incompressible RANS equations. In Section 3, we present a semi-discrete divergence-conforming HDG method for the incompressible RANS equations, assuming an underlying linear eddy viscosity model, and we show that the method conserves mass in a point-wise manner, conserves momentum globally, and is energy stable. Motivated by the fact that most turbulence model equations are transport equations, we recall the strong form of the scalar advection-diffusion equation in Section 4, and we present a semi-discrete HDG method for the advection-diffusion equation, assuming an underlying divergence-free velocity field, in Section 5. In Section 6, we present a semi-discrete HDG method for the incompressible RANS equations with the Spalart-Allmaras one-equation turbulence model by combining our previous semi-discrete HDG methods for the incompressible RANS equations and the advection-diffusion equation. In Section 7, we present a staggered time-integration scheme for our HDG method for the incompressible RANS equations with the Spalart-Allmaras model, and we illustrate how static condensation can be employed to reduce computational cost in Section 8. We present illustrative numerical results in Section 9, and finally, in Section 10, we draw conclusions and discuss future research directions.

We begin this paper by recalling the strong form of the incompressible RANS equations. The derivation of the incompressible RANS equations relies on a Reynolds decomposition of the velocity and pressure fields into mean and fluctuating components, viz.:

	$\displaystyle\bf{u}$	$\displaystyle=\overline{\bf{u}}+\bf{u}^{\prime}$		(1)
	$\displaystyle p$	$\displaystyle=\overline{p}+p^{\prime}$

where $\overline{*}$ and $*^{\prime}$ denote the mean and fluctuation respectively. We then obtain the incompressible RANS equations by taking the mean of the incompressible Navier-Stokes equations and exploiting the above Reynolds decompositions. The resulting system is displayed below:

		$\displaystyle\frac{\partial\overline{\bf{u}}}{\partial t}+\nabla\cdot(\overline{\bf{u}}\otimes\overline{\bf{u}})+\frac{1}{\rho}\nabla\overline{p}-\nabla\cdot(2\nu\nabla^{S}\bf{u})+\nabla\cdot(\overline{\bf{u}^{\prime}\otimes\bf{u}^{\prime}})=\bf{f}$		(2)
		$\displaystyle\nabla\cdot\overline{\bf{u}}=0$

Above, $\nabla^{S}\bf{u}=\frac{1}{2}(\nabla\overline{\bf{u}}+(\nabla\overline{\bf{u}})^{T})$ denotes the mean strain rate tensor and $\overline{\bf{u}^{\prime}\otimes\bf{u}^{\prime}}$ denotes the Reynolds stress tensor.

In this paper, we will focus our attention on linear eddy viscosity models. These models assume that the mean strain rate tensor and the anisotropic part of the Reynolds stress tensors are aligned, and hence there exists an eddy viscosity $\nu_{T}$ such that:

$\displaystyle\overline{\bf{u}^{\prime}\otimes\bf{u}^{\prime}}=-2\nu_{T}\nabla^{S}{\overline{\bf{u}}}+\frac{2}{3}k\bf{I}$

(3)

where $k=\frac{1}{2}\textup{tr}\left(\overline{\bf{u}^{\prime}\otimes\bf{u}^{\prime}}\right)$ is the turbulent kinetic energy and $\bf{I}$ is the identity tensor. It is generally assumed that the eddy viscosity is non-negative. Otherwise, the resulting system may be unstable. In what follows, we combine the contributions of the mean pressure field and the isotropic part of the Reynolds stress tensor using a modified pressure $\overline{p}+\frac{2}{3}\rho k$. With an abuse of notation, we use $\overline{p}$ to denote this modified pressure.

We are now ready to state the strong form of the incompressible RANS equations subject to a linear eddy viscosity model. Let $\Omega$ denote a Lipschitz and bounded domain in $\mathbb{R}^{d}$, where $d=2$ for two-dimensional domains and $d=3$ for three-dimensional domains. Let $\partial\Omega$ denote the boundary of $\Omega$ with outward unit normal ${\bf{n}}$. Let $\partial\Omega$ be partitioned into a Dirichlet boundary $\Gamma_{D}$ and a Neumann boundary $\Gamma_{N}$ such that $\partial\Omega=\overline{\Gamma_{D}\cup\Gamma_{N}}$ and $\Gamma_{D}\cap\Gamma_{N}=\emptyset$. We assume a constant density $\rho\in\mathbb{R}^{+}$, a variable kinematic viscosity $\nu:\Omega\times(0,\infty)\rightarrow\mathbb{R}^{+}$, a variable body force ${\bf{f}}:\Omega\times(0,\infty)\rightarrow\mathbb{R}^{d}$, a variable velocity specification on the Dirichlet boundary ${\bf{g}}:\Gamma_{D}\times(0,\infty)\rightarrow\mathbb{R}^{d}$, a variable traction specification on the Neumann boundary ${\bf{h}}:\Gamma_{N}\times(0,\infty)\rightarrow\mathbb{R}^{d}$, and a variable initial velocity $\overline{\bf{u}}_{0}:\Omega\rightarrow\mathbb{R}^{d}$. We also freeze the turbulent viscosity in the following presentation. That is, we assume that $\nu_{T}:\Omega\times(0,\infty)\rightarrow\mathbb{R}^{+}$ is a known function. This will allow us to examine the properties of our discretization for the incompressible RANS equations independent of the chosen eddy viscosity model. We will later make the turbulent viscosity a function of the unknown mean flow field and unknown turbulence variables. With the above assumptions, the strong form is as follows: Strong Form for the Incompressible RANS Equations

Find $\overline{\bf{u}}:\bar{\Omega}\times[0,\infty)\rightarrow\mathbb{R}^{d}$ and $\overline{p}:\Omega\times(0,\infty)\rightarrow\mathbb{R}$ such that: $\displaystyle\frac{\partial\overline{\bf{u}}}{\partial t}+\nabla\cdot(\overline{\bf{u}}\otimes\overline{\bf{u}})+\frac{1}{\rho}\nabla\overline{p}-\nabla\cdot(2(\nu+\nu_{T})\nabla^{S}\overline{\bf{u}})$ $\displaystyle=\bf{f}$ $\displaystyle\text{ in }\Omega\times(0,\infty)$ (4) $\displaystyle\nabla\cdot\overline{\bf{u}}$ $\displaystyle=0$ $\displaystyle\text{ in }\Omega\times(0,\infty)$ $\displaystyle\overline{\bf{u}}$ $\displaystyle={\bf{g}}$ $\displaystyle\text{ on }\Gamma_{D}\times(0,\infty)$ $\displaystyle\left[-\frac{1}{\rho}\overline{p}\textbf{I}+2(\nu+\nu_{T})\nabla^{S}\overline{\bf{u}}\right]\cdot{\bf{n}}-\text{min}(\overline{\bf{u}}\cdot{\bf{n}},0)\overline{\bf{u}}$ $\displaystyle={\bf{h}}$ $\displaystyle\text{ on }\Gamma_{N}\times(0,\infty)$ $\displaystyle\overline{\bf{u}}(\cdot,0)$ $\displaystyle=\overline{\bf{u}}_{0}$ $\displaystyle\text{ in }\Omega$

Note that we have utilized somewhat unorthodox Neumann boundary conditions above. On the outflow parts $\Gamma_{N}$ (i.e., where $\overline{\bf{u}}\cdot{\bf{n}}\geq 0$), we set the traction, i.e., $\left[-\frac{1}{\rho}\overline{p}\textbf{I}+2(\nu+\nu_{T})\nabla^{S}\overline{\bf{u}}\right]\cdot{\bf{n}}={\bf{h}}$, as is standard. However, on the inflow parts of $\Gamma_{N}$, we instead set the sum of momentum, pressure, and diffusive fluxes, i.e., $\left[-(\overline{\bf{u}}\otimes\overline{\bf{u}})-\frac{1}{\rho}\overline{p}\textbf{I}+2(\nu+\nu_{T})\nabla^{S}\overline{\bf{u}}\right]\cdot{\bf{n}}={\bf{h}}$. This yields a well-posed formulation in the presence of backflow [30].

Let $P_{k}(D)$ denote the space of polynomials of degree $k\geq 0$ on a domain $D$. For polynomial degree $k\geq 1$, we consider the following finite element spaces defined on the mesh:

	$\displaystyle V^{h}$	$\displaystyle=\{{\bf{v}}^{h}\in[L^{2}(\Omega)]^{d}:\left.{\bf{v}}^{h}\right|_{\Omega_{e}}\in[P_{k}(\Omega_{e})]^{d}\hskip 5.69054pt\forall\Omega_{e}\in\mathcal{T}\}$		(5)
	$\displaystyle\hat{V}^{h}$	$\displaystyle=\{\hat{\bf{v}}^{h}\in[L^{2}(\tilde{\Gamma})]^{d}:\left.\hat{\bf{v}}^{h}\right|_{F}\in[P_{k}(F)]^{d}\hskip 5.69054pt\forall F\in\mathcal{F}\}$
	$\displaystyle Q^{h}$	$\displaystyle=\{q^{h}\in L^{2}(\Omega):\left.q^{h}\right|_{\Omega_{e}}\in P_{k-1}(\Omega_{e})\hskip 5.69054pt\forall\Omega_{e}\in\mathcal{T}\}$
	$\displaystyle\hat{Q}^{h}$	$\displaystyle=\{\hat{q}^{h}\in L^{2}(\tilde{\Gamma}):\left.\hat{q}^{h}\right|_{F}\in P_{k}(F)\hskip 5.69054pt\forall F\in\mathcal{F}\}$

We will approximate the velocity and pressure fields over element interiors using the spaces $V^{h}$ and $Q^{h}$, respectively. We will approximate the velocity and pressure fields over the mesh skeleton using the spaces $\hat{V}^{h}$ and $\hat{Q}^{h}$, respectively. We will henceforth refer to $V^{h}$ and $Q^{h}$ as the (discrete) velocity and pressure spaces and we will refer to $\hat{V}^{h}$ and $\hat{Q}^{h}$ as the (discrete) velocity and pressure trace spaces. It is important to note that functions in the velocity and pressure spaces are defined on the whole domain $\Omega$, whereas functions in the trace spaces are defined only on the mesh skeleton $\tilde{\Gamma}$. We also introduce the spaces:

	$\displaystyle\hat{V}_{g}^{h}$	$\displaystyle=\{\hat{\bf{v}}^{h}\in\hat{V}^{h}:\hat{\bf{v}}^{h}=\textbf{g}\text{ on }\Gamma_{D}\}$		(6)
	$\displaystyle\hat{V}_{0}^{h}$	$\displaystyle=\{\hat{\bf{v}}^{h}\in\hat{V}^{h}:\hat{\bf{v}}^{h}=\textbf{0}\text{ on }\Gamma_{D}\}$

which will correspond to our velocity trace trial and test spaces.

Note that functions in the velocity and pressure spaces $V^{h}$ and $Q^{h}$ may be discontinuous across cell boundaries. Let $F$ be an interior facet of the mesh shared by two elements $\Omega_{e^{+}}$ and $\Omega_{e^{-}}$. We denote the outward facing normals on $F$ to $\Omega_{e^{+}}$ and $\Omega_{e^{-}}$ as $\bf{n}^{+}$ and $\bf{n}^{-}$, respectively. For ${\bf{y}}\in V^{h}$, we denote the trace of ${\bf{y}}|_{\Omega_{e^{+}}}$ along $F$ as ${\bf{y}}^{+}$, and we denote the trace of ${\bf{y}}|_{\Omega_{e^{-}}}$ along $F$ as ${\bf{y}}^{-}$. The jump operator across $F$ can then be defined as $\llbracket{\bf{y}}\rrbracket={\bf{y}}^{+}\cdot{\bf{n}}^{+}+{\bf{y}}^{-}\cdot{\bf{n}}^{-}$.

With all of the aforementioned terminology in place, we are now ready to state our semi-discrete HDG formulation for the incompressible RANS equations. Our method is the natural extension of the divergence-conforming HDG method of Rhebergen and Wells [21] to the RANS setting, wherein advective fluxes are treated using upwinding and diffusive fluxes are treated using the symmetric interior penalty method.
Semi-Discrete HDG Formulation for the Incompressible RANS Equations

Find $\left(\overline{{\bf{u}}}^{h},\hat{\overline{{\bf{u}}}}^{h},\overline{p}^{h},\hat{\overline{p}}^{h}\right)\in V^{h}\times\hat{V}_{g}^{h}\times Q^{h}\times\hat{Q}^{h}$ such that:

Continuity Equation $\displaystyle\sum_{e}\int_{\Omega_{e}}\frac{1}{\rho}\nabla\cdot\overline{{\bf{u}}}^{h}q^{h}d\Omega=0\hskip 227.62204pt\forall q^{h}\in Q^{h}$ (7) Continuity Conservativity Condition $\displaystyle\sum_{e}\sum_{i}\int_{\Gamma_{e_{i}}}\frac{1}{\rho}\overline{{\bf{u}}}^{h}\cdot{\bf{n}}\hat{q}^{h}d\Gamma-\int_{\partial\Omega}\frac{1}{\rho}\hat{\overline{{\bf{u}}}}^{h}\cdot{\bf{n}}\hat{q}^{h}d\Gamma=0\hskip 126.61476pt\forall\hat{q}^{h}\in\hat{Q}^{h}$ (8) Momentum Equation $\displaystyle\int_{\Omega}\frac{\partial\overline{{\bf{u}}}^{h}}{\partial t}\cdot{\bf{v}}^{h}d\Omega$ (9) $\displaystyle-\sum_{e}$ $\displaystyle\int_{\Omega_{e}}(\overline{{\bf{u}}}^{h}\otimes\overline{{\bf{u}}}^{h}):\nabla{\bf{v}}^{h}d\Omega-\sum_{e}\int_{\Omega_{e}}\frac{1}{\rho}\overline{p}^{h}\textbf{I}:\nabla{\bf{v}}^{h}d\Omega+\sum_{e}\int_{\Omega_{e}}2(\nu+\nu_{T})\nabla^{S}\overline{{\bf{u}}}^{h}:\nabla^{S}{{\bf{v}}}^{h}d\Omega$ $\displaystyle+\sum_{e}\sum_{i}$ $\displaystyle\int_{\Gamma_{e_{i}}}(\overline{{\bf{u}}}^{h}\otimes\overline{{\bf{u}}}^{h}):({\bf{v}}^{h}\otimes{\bf{n}})d\Gamma+\sum_{e}\sum_{i}\int_{\Gamma_{e_{i}}}(\hat{\overline{{\bf{u}}}}^{h}-\overline{{\bf{u}}}^{h})\otimes\lambda\overline{{\bf{u}}}^{h}:({\bf{v}}^{h}\otimes{\bf{n}})d\Gamma$ $\displaystyle+\sum_{e}\sum_{i}$ $\displaystyle\int_{\Gamma_{e_{i}}}\frac{1}{\rho}\hat{\overline{p}}^{h}\textbf{I}:({\bf{v}}^{h}\otimes{\bf{n}})d\Gamma-\sum_{e}\sum_{i}\int_{\Gamma_{e_{i}}}2(\nu+\nu_{T})\nabla^{S}\overline{{\bf{u}}}^{h}:({\bf{v}}^{h}\otimes{\bf{n}})d\Gamma$ $\displaystyle-\sum_{e}\sum_{i}$ $\displaystyle\int_{\Gamma_{e_{i}}}{\frac{2C_{pen}}{h_{e}}(\nu+\nu_{T})}(\hat{\overline{{\bf{u}}}}^{h}-\overline{{\bf{u}}}^{h})\otimes{\bf{n}}:({\bf{v}}^{h}\otimes{\bf{n}})d\Gamma$ $\displaystyle+\sum_{e}\sum_{i}$ $\displaystyle\int_{\Gamma_{e_{i}}}2(\nu+\nu_{T})[(\hat{\overline{{\bf{u}}}}^{h}-\overline{{\bf{u}}}^{h})\otimes{\bf{n}}]:\nabla^{S}{{\bf{v}}}^{h}d\Gamma-\int_{\Omega}{\bf{f}}\cdot{\bf{v}}^{h}d\Omega=0\hskip 56.9055pt\forall{\bf{v}}^{h}\in{V}^{h}$ Momentum Conservativity Condtion $\displaystyle-\sum_{e}\sum_{i}$ $\displaystyle\int_{\Gamma_{e_{i}}}(\overline{{\bf{u}}}^{h}\otimes\overline{{\bf{u}}}^{h}):(\hat{{{\bf{v}}}}^{h}\otimes{\bf{n}})d\Gamma-\sum_{e}\sum_{i}\int_{\Gamma_{e_{i}}}(\hat{\overline{{\bf{u}}}}^{h}-\overline{{\bf{u}}}^{h})\otimes\lambda\overline{{\bf{u}}}^{h}:(\hat{{{\bf{v}}}}^{h}\otimes{\bf{n}})d\Gamma$ (10) $\displaystyle-\sum_{e}\sum_{i}$ $\displaystyle\int_{\Gamma_{e_{i}}}\frac{1}{\rho}\hat{\overline{p}}^{h}\textbf{I}:(\hat{{{\bf{v}}}}^{h}\otimes{\bf{n}})d\Gamma+\sum_{e}\sum_{i}\int_{\Gamma_{e_{i}}}2(\nu+\nu_{T})\nabla^{S}\overline{{\bf{u}}}^{h}:(\hat{{{\bf{v}}}}^{h}\otimes{\bf{n}})d\Gamma$ $\displaystyle+\sum_{e}\sum_{i}$ $\displaystyle\int_{\Gamma_{e_{i}}}{\frac{2C_{pen}}{h_{e}}(\nu+\nu_{T})}(\hat{\overline{{\bf{u}}}}^{h}-\overline{{\bf{u}}}^{h})\otimes{\bf{n}}:(\hat{{{\bf{v}}}}^{h}\otimes{\bf{n}})d\Gamma$ $\displaystyle+$ $\displaystyle\int_{\Gamma_{N}}(1-\lambda)(\hat{\overline{{\bf{u}}}}^{h}\cdot{\bf{n}})(\hat{\overline{{\bf{u}}}}^{h}\cdot\hat{{{\bf{v}}}}^{h})d\Gamma+\int_{\Gamma_{N}}{\bf{h}}\cdot\hat{{{\bf{v}}}}^{h}d\Gamma=0\hskip 81.09035pt\forall\hat{{{\bf{v}}}}^{h}\in\hat{{V}}_{0}^{h}$

Above, on the $i^{\text{th}}$ facet $\Gamma_{e_{i}}$ of the $e^{\text{th}}$ element $\Omega_{e}$, $\lambda$ is an indicator function that takes on a value of unity if the facet is an inflow facet and a value of zero if it is an outflow facet, i.e.:

$\displaystyle\lambda=$

(11)

These definitions result from the fact that we have discretized the incompressible RANS equations using an upwind treatment of the advective component of the flux, as discussed in [21].

The penalty constant $C_{pen}$ arises from the fact that we have discretized the diffusive component of the flux using the symmetric interior penalty method [31]. As is typical with interior penalty methods, the constant $C_{pen}$ must be chosen sufficiently large as to ensure the resulting semi-discrete formulation is energy stable [32]. To arrive at an intelligent choice for $C_{pen}$, one typically turns to trace inequalities. Following Warburton and Hesthaven [33], it can be shown that the following inequality holds in the two-dimensional setting:

$$||\nabla^{S}{\bf{v}^{h}}||_{\partial\Omega_{e}}^{2}\leq(k+1)(k+2)\frac{\text{Perimeter}(\Omega_{e})}{\text{Area}(\Omega_{e})}||\nabla^{S}{\bf{v}^{h}}||_{\Omega_{e}}^{2}$$

(12)

for all $\Omega_{e}\in\mathcal{T}$ and ${\bf{v}}^{h}\in V^{h}$, and a similar result holds in the three-dimensional setting. If we define the mesh size $h_{e}$ to be equal to $\text{Area}(\Omega_{e})/\text{Perimeter}(\Omega_{e})$, it suffices to select $C_{pen}\geq(k+1)(k+2)$. We later demonstrate this choice yields an energy stable method. In all of our later computations, we have selected $C_{pen}=(k+1)(k+2)$. Note that the definition $h_{e}=\text{Area}(\Omega_{e})/\text{Perimeter}(\Omega_{e})$ is somewhat non-standard. It is related not to the diameter of the element but rather to the diameter of the incircle of the element. However, we have found this definition to be critical in ensuring stability in the presence of boundary layer meshes containing high-aspect ratio elements.

We now present a consistency result for our semi-discrete formulation.

Proposition 1 (Consistency)

The semi-discrete HDG method presented in Eqs. (7-10) is consistent provided the exact solution $(\bf{u},p)$ of the incompressible RANS equations is sufficiently smooth. That is, Eqs. (7-10) hold if we replace $\left(\overline{{\bf{u}}}^{h},\hat{\overline{{\bf{u}}}}^{h},\overline{p}^{h},\hat{\overline{p}}^{h}\right)$ with $\left(\overline{\bf{u}},\overline{\bf{u}}|_{\tilde{\Gamma}},\overline{p},\overline{p}|_{\tilde{\Gamma}}\right)$.

Proof. Note that for, a sufficiently smooth exact solution $({\bf{u}},p)$, Eq. (4) are satisfied in a point-wise manner. Using this, we show that each of Eqs. (7-10) hold if we replace $\left(\overline{{\bf{u}}}^{h},\hat{\overline{{\bf{u}}}}^{h},\overline{p}^{h},\hat{\overline{p}}^{h}\right)$ with $\left(\overline{\bf{u}},\overline{\bf{u}}|_{\tilde{\Gamma}},\overline{p},\overline{p}|_{\tilde{\Gamma}}\right)$. We begin with the continuity equation, that is, Eq. (7). This holds trivially since $\overline{{\bf{u}}}$ is divergence-free and hence:

$\displaystyle\sum_{e}\int_{\Omega_{e}}\frac{1}{\rho}\nabla\cdot\overline{{\bf{u}}}q^{h}d\Omega=0\ \hskip 28.45274pt\forall q^{h}\in Q^{h}$

(13)

We next consider the continuity conservativity equation, Eq. (8). This holds since:

$\displaystyle\sum_{e}\sum_{i}\int_{\Gamma_{e_{i}}}\frac{1}{\rho}\overline{{\bf{u}}}\cdot{\bf{n}}\hat{q}^{h}d\Gamma-\int_{\partial\Omega}\frac{1}{\rho}\overline{{\bf{u}}}|_{\tilde{\Gamma}}\cdot{\bf{n}}\hat{q}^{h}d\Gamma=\int_{\tilde{\Gamma}\backslash\partial\Omega}\frac{1}{\rho}\llbracket\overline{{\bf{u}}}\rrbracket\hat{q}^{h}d\Gamma=0\ \hskip 28.45274pt\forall\hat{q}^{h}\in\hat{Q}^{h}$

(14)

as $\llbracket\overline{{\bf{u}}}\rrbracket=0$ on each interior facet $F\in\mathcal{F}_{int}$. We now continue with the momentum equation, Eq. (9). Through reverse integration by parts, we have:

		$\displaystyle\int_{\Omega}\frac{\partial\overline{{\bf{u}}}}{\partial t}\cdot{\bf{v}}^{h}d\Omega-\sum_{e}\int_{\Omega_{e}}(\overline{{\bf{u}}}\otimes\overline{{\bf{u}}}):\nabla{\bf{v}}^{h}d\Omega-\sum_{e}\int_{\Omega_{e}}\frac{1}{\rho}\overline{p}\textbf{I}:\nabla{\bf{v}}^{h}d\Omega+\sum_{e}\int_{\Omega_{e}}2(\nu+\nu_{T})\nabla^{S}\overline{{\bf{u}}}:\nabla^{S}{{\bf{v}}}^{h}d\Omega$		(15)
		$\displaystyle+\sum_{e}\sum_{i}\int_{\Gamma_{e_{i}}}(\overline{{\bf{u}}}\otimes\overline{{\bf{u}}}):({\bf{v}}^{h}\otimes{\bf{n}})d\Gamma+\sum_{e}\sum_{i}\int_{\Gamma_{e_{i}}}(\overline{{\bf{u}}}|_{\tilde{\Gamma}}-\overline{{\bf{u}}})\otimes\lambda\overline{{\bf{u}}}:({\bf{v}}^{h}\otimes{\bf{n}})d\Gamma+\sum_{e}\sum_{i}\int_{\Gamma_{e_{i}}}\frac{1}{\rho}\overline{p}|_{\tilde{\Gamma}}\textbf{I}:({\bf{v}}^{h}\otimes{\bf{n}})d\Gamma$
		$\displaystyle-\sum_{e}\sum_{i}\int_{\Gamma_{e_{i}}}2(\nu+\nu_{T})\nabla^{S}\overline{{\bf{u}}}:({\bf{v}}^{h}\otimes{\bf{n}})d\Gamma-\sum_{e}\sum_{i}\int_{\Gamma_{e_{i}}}{\frac{2C_{pen}}{h_{e}}(\nu+\nu_{T})}(\overline{{\bf{u}}}|_{\tilde{\Gamma}}-\overline{{\bf{u}}})\otimes{\bf{n}}:({\bf{v}}^{h}\otimes{\bf{n}})d\Gamma$
		$\displaystyle+\sum_{e}\sum_{i}\int_{\Gamma_{e_{i}}}2(\nu+\nu_{T})[(\overline{{\bf{u}}}|_{\tilde{\Gamma}}-\overline{{\bf{u}}})\otimes{\bf{n}}]:\nabla^{S}{{\bf{v}}}^{h}d\Gamma-\int_{\Omega}{\bf{f}}\cdot{\bf{v}}^{h}d\Omega$
		$\displaystyle=\sum_{e}\int_{\Omega_{e}}\left[\frac{\partial\overline{{\bf{u}}}}{\partial t}+\nabla\cdot(\overline{{\bf{u}}}\otimes\overline{{\bf{u}}})+\frac{1}{\rho}\nabla\overline{p}-\nabla\cdot(2(\nu+\nu_{T})\nabla^{S}\overline{{\bf{u}}})-{\bf{f}}\right]\cdot{\bf{v}}^{h}d\Omega=\vec{0}\hskip 28.45274pt\forall{\bf{v}}^{h}\in{V}^{h}$

since:

$\displaystyle\frac{\partial\overline{{\bf{u}}}}{\partial t}+\nabla\cdot(\overline{{\bf{u}}}\otimes\overline{{\bf{u}}})+\frac{1}{\rho}\nabla\overline{p}-\nabla\cdot(2(\nu+\nu_{T})\nabla^{S}\overline{{\bf{u}}})-{\bf{f}}={\bf{0}}$

$\displaystyle\text{ in }\Omega\times(0,\infty)$

(16)

We finish with the momentum conservativity equation, Eq. (10). This holds since:

		$\displaystyle-\sum_{e}\sum_{i}\int_{\Gamma_{e_{i}}}(\overline{{\bf{u}}}\otimes\overline{{\bf{u}}}):(\hat{{{\bf{v}}}}^{h}\otimes{\bf{n}})d\Gamma-\sum_{e}\sum_{i}\int_{\Gamma_{e_{i}}}(\overline{{\bf{u}}}|_{\tilde{\Gamma}}-\overline{{\bf{u}}})\otimes\lambda\overline{{\bf{u}}}^{h}:(\hat{{{\bf{v}}}}^{h}\otimes{\bf{n}})d\Gamma$		(17)
		$\displaystyle-\sum_{e}\sum_{i}\int_{\Gamma_{e_{i}}}\frac{1}{\rho}\overline{p}|_{\tilde{\Gamma}}\textbf{I}:(\hat{{{\bf{v}}}}^{h}\otimes{\bf{n}})d\Gamma+\sum_{e}\sum_{i}\int_{\Gamma_{e_{i}}}2(\nu+\nu_{T})\nabla^{S}\overline{{\bf{u}}}:(\hat{{{\bf{v}}}}^{h}\otimes{\bf{n}})d\Gamma$
		$\displaystyle+\sum_{e}\sum_{i}\int_{\Gamma_{e_{i}}}{\frac{2C_{pen}}{h_{e}}(\nu+\nu_{T})}(\overline{{\bf{u}}}|_{\tilde{\Gamma}}-\overline{{\bf{u}}})\otimes{\bf{n}}:(\hat{{{\bf{v}}}}^{h}\otimes{\bf{n}})d\Gamma$
		$\displaystyle+\int_{\Gamma_{N}}(1-\lambda)(\overline{{\bf{u}}}|_{\tilde{\Gamma}}\cdot{\bf{n}})(\overline{{\bf{u}}}|_{\tilde{\Gamma}}\cdot\hat{{{\bf{v}}}}^{h})d\Gamma+\int_{\Gamma_{N}}{\bf{h}}\cdot\hat{{{\bf{v}}}}^{h}d\Gamma$
		$\displaystyle=\int_{\tilde{\Gamma}/\partial\Omega}\left(-\left(\overline{\bf{u}}\cdot\hat{\bf{v}}^{h}\right)\llbracket\overline{\bf{u}}\rrbracket+\llbracket 2(\nu+\nu_{T})\nabla^{s}\bar{\bf{u}}\rrbracket\cdot\hat{\bf{v}}^{h}\right)d\Gamma$
		$\displaystyle-\int_{\Gamma_{N}}\left(\left[-\frac{1}{\rho}\overline{p}\textbf{I}+2(\nu+\nu_{T})\nabla^{S}\overline{\bf{u}}\right]\cdot{\bf{n}}-\text{min}(\overline{\bf{u}}\cdot{\bf{n}},0)\overline{\bf{u}}-{\bf{h}}\right)\cdot\hat{{\bf{v}}}^{h}d\Gamma=0\hskip 54.06023pt\forall\hat{{{\bf{v}}}}^{h}\in\hat{{V}}^{h}$

as $\llbracket\overline{{\bf{u}}}\rrbracket=0$ and $\llbracket 2(\nu+\nu_{T})\nabla^{S}\overline{{\bf{u}}}\rrbracket={\bf{0}}$ on each interior facet $F\in\mathcal{F}_{int}$ and:

$\displaystyle\left[-\frac{1}{\rho}\overline{p}\textbf{I}+2(\nu+\nu_{T})\nabla^{S}\overline{\bf{u}}\right]\cdot{\bf{n}}-\text{min}(\overline{\bf{u}}\cdot{\bf{n}},0)\overline{\bf{u}}-{\bf{h}}$

$\displaystyle={\bf{0}}$

$\displaystyle\text{ on }\Gamma_{N}\times(0,\infty)$

(18)

This completes the proof.∎

We next present a result illustrating our semi-discrete formulation conserves mass in a pointwise manner.

Proposition 2 (Pointwise Mass Conservation)

If $\overline{{\bf{u}}}^{h}\in V^{h}$ and $\hat{\overline{{\bf{u}}}}^{h}\in\hat{V}^{h}$ satisfy Eqs. (7) and (8), with $V_{h}$ and $\hat{V}^{h}$ defined in Eq .(6), then:

$$\nabla\cdot\overline{{\bf{u}}}^{h}=0\hskip 14.22636pt\forall{\bf{x}}\in\Omega_{e},\hskip 14.22636pt\forall\Omega_{e}\in\mathcal{T}$$

(19)

and:

	$\displaystyle\llbracket\overline{{\bf{u}}}^{h}\rrbracket=0$	$\displaystyle\forall{\bf{x}}\in F,\hskip 14.22636pt\forall F\in\mathcal{F}_{\text{int}}$		(20)
	$\displaystyle\overline{{\bf{u}}}^{h}\cdot{\bf{n}}=\hat{\overline{{\bf{u}}}}^{h}\cdot{\bf{n}}$	$\displaystyle\forall{\bf{x}}\in F,\hskip 14.22636pt\forall F\in\mathcal{F}_{\text{bdy}}$

Proof. The proof follows in the same manner as the proof of Proposition 1 in [21], though we review the proof as pointwise mass conservation is a key feature of our HDG method. From Eq. (7) it follows that:

$$0=\int_{\Omega_{e}}q^{h}\nabla\cdot\overline{{\bf{u}}}^{h}d\Omega\hskip 14.22636pt\forall q^{h}\in P_{k-1}(\Omega_{e}),\hskip 14.22636pt\forall\Omega_{e}\in\mathcal{T}$$

(21)

Since $\nabla\cdot\overline{{\bf{u}}}^{h}\in P_{k-1}(\Omega_{e})$ for every $\Omega_{e}\in\mathcal{T}$, we can take $q^{h}=\nabla\cdot\overline{{\bf{u}}}^{h}$ in the above equation, yielding $\int_{\Omega_{e}}(\nabla\cdot\overline{{\bf{u}}}^{h})^{2}d\Omega=0$ for every $\Omega_{e}\in\mathcal{T}$. Thus $\nabla\cdot\overline{{\bf{u}}}^{h}\equiv 0$ in $\Omega$.

From Eq. (8) it follows that:

$$0=\int_{F}\llbracket\overline{{\bf{u}}}^{h}\rrbracket\hat{q}^{h}dS\hskip 14.22636pt\forall\hat{q}^{h}\in P_{k}(F),\hskip 14.22636pt\forall F\in\mathcal{F}_{\text{int}}$$

(22)

Since $\llbracket\overline{{\bf{u}}}^{h}\rrbracket\in P_{k}(F)$ for all $F\in\mathcal{F}_{\text{int}}$, we can choose $\hat{q}^{h}=\llbracket\overline{{\bf{u}}}^{h}\rrbracket$, yielding $\int_{F}\llbracket\overline{{\bf{u}}}^{h}\rrbracket^{2}d\Gamma=0$ for all $F\in\mathcal{F}_{\text{int}}$. Thus $\llbracket\overline{{\bf{u}}}^{h}\rrbracket\equiv 0$ on $\tilde{\Gamma}/\partial\Omega$.

From Eq. (8) it also follows that:

$$0=\int_{F}\left(\overline{{\bf{u}}}^{h}-\hat{\overline{{\bf{u}}}}^{h}\right)\cdot{\bf{n}}\hat{q}^{h}dS\hskip 14.22636pt\forall\hat{q}^{h}\in P_{k}(F),\hskip 14.22636pt\forall F\in\mathcal{F}_{\text{bdy}}$$

(23)

Since $\left(\overline{{\bf{u}}}^{h}-\hat{\overline{{\bf{u}}}}^{h}\right)\cdot{\bf{n}}\in P_{k}(F)$ for all $F\in\mathcal{F}_{\text{bdy}}$, we can choose $\hat{q}^{h}=\left(\overline{{\bf{u}}}^{h}-\hat{\overline{{\bf{u}}}}^{h}\right)\cdot{\bf{n}}$, yielding $\int_{F}\left(\left(\overline{{\bf{u}}}^{h}-\hat{\overline{{\bf{u}}}}^{h}\right)\cdot{\bf{n}}\right)^{2}d\Gamma=0$ for all $F\in\mathcal{F}_{\text{bdy}}$. Thus $\left(\overline{{\bf{u}}}^{h}-\hat{\overline{{\bf{u}}}}^{h}\right)\cdot{\bf{n}}\equiv 0$ on $\partial\Omega$. ∎

Proposition 2 is a stronger statement of mass conservation than general DG or HDG finite element methods, in which mass conservation is normally satisfied locally (element-wise) in an integral sense only. There are two fundamental reasons why our semi-discrete formulation conserves mass in a pointwise manner. First of all, the element-wise divergence operator maps the velocity space $V^{h}$ into the pressure space $Q^{h}$. Second of all, the element-wise trace operator maps the velocity space $V^{h}$ into the pressure trace space $\hat{Q}^{h}$. Any choice of velocity, pressure, velocity trace, and velocity pressure spaces that satisfy these two properties will result in a divergence-conforming HDG method. It can easily be seen that the choice of approximation spaces presented here corresponds to a hybridization of a divergence-conforming DG method employing a Brezzi-Douglas-Marini velocity-pressure pair [34].

The next result gives a global momentum balance law for our semi-discrete method.

Proposition 3 (Global Momentum Balance)

Let $(\overline{{\bf{u}}}^{h},\hat{\overline{{\bf{u}}}}^{h},\overline{p}^{h},\hat{\overline{p}}^{h})\in V^{h}\times\hat{V}^{h}\times Q^{h}\times\hat{Q}^{h}$ satisfy Eq.(7 -10). If $\Gamma_{D}=\emptyset$:

$\displaystyle\frac{d}{dt}$

$\displaystyle\int_{\Omega}\overline{{\bf{u}}}^{h}d\Omega=\int_{\Omega}{\bf{f}}d\Omega-\int_{\Gamma_{N}}(1-\lambda)(\hat{\overline{{\bf{u}}}}^{h}\cdot{\bf{n}})\hat{\overline{{\bf{u}}}}^{h}d\Gamma-\int_{\Gamma_{N}}{\bf{h}}d\Gamma$

(24)

Proof. The proof follows in the same manner as the proof of Proposition 2 in [21]. ∎

The next result states that our semi-discrete method is globally energy stable.

Proposition 4 (Global Energy Stability)

If $(\overline{{\bf{u}}}^{h},\hat{\overline{{\bf{u}}}}^{h},\overline{p}^{h},\hat{\overline{p}}^{h})\in V^{h}\times\hat{V}^{h}\times Q^{h}\times\hat{Q}^{h}$ satisfy Eqs. (7 - 10), ${\bf{f}}={\bf{0}}$, ${\bf{g}}={\bf{0}}$, ${\bf{h}}={\bf{0}}$, and $\nu_{T}\geq 0$:

$$\frac{d}{dt}\sum_{e}\int_{\Omega_{e}}|\overline{{\bf{u}}}^{h}|^{2}d\Omega_{e}\leq 0$$

(25)

Proof. The proof follows in the same manner as the proof of Proposition 3 in [21], though we review the proof to illustrate the role of the penalty constant $C_{pen}$. Setting ${\bf{v}}^{h}=\overline{{\bf{u}}}^{h}$ and $\hat{{\bf{v}}}^{h}=\hat{\overline{{\bf{u}}}}^{h}$ in Eqs. (7-10), summing the two expressions together, recognizing that $\nabla\cdot\overline{\bf{u}}^{h}=0$ in $\Omega$ and $\overline{{\bf{u}}}^{h}\cdot{\bf{n}}=\hat{\overline{{\bf{u}}}}^{h}\cdot{\bf{n}}$ for all $F\in\mathcal{F}$, and exploiting that $\lambda\overline{{\bf{u}}}^{h}\cdot{\bf{n}}=(\overline{{\bf{u}}}^{h}\cdot{\bf{n}}-|\overline{{\bf{u}}}^{h}\cdot{\bf{n}}|)/2$, we obtain:

	$\displaystyle\sum_{e}\frac{1}{2}$	$\displaystyle\int_{\Omega_{e}}\frac{\partial|\overline{{\bf{u}}}^{h}|^{2}}{\partial t}d\Omega+\sum_{e}\sum_{i}\frac{1}{2}\int_{\Gamma_{e_{i}}}(\overline{{\bf{u}}}^{h}\cdot{\bf{n}})|\overline{{\bf{u}}}^{h}|^{2}d\Gamma$		(26)
	$\displaystyle-\sum_{e}\sum_{i}\frac{1}{2}$	$\displaystyle\int_{\Gamma_{e_{i}}}(\overline{{\bf{u}}}^{h}\cdot{\bf{n}})|\hat{\overline{{\bf{u}}}}^{h}|^{2}d\Gamma+\sum_{e}\sum_{i}\frac{1}{2}\int_{\Gamma_{e_{i}}}|\overline{{\bf{u}}}^{h}\cdot{\bf{n}}||\overline{{\bf{u}}}^{h}-\hat{\overline{{\bf{u}}}}^{h}|^{2}d\Gamma$
	$\displaystyle+\sum_{e}$	$\displaystyle\int_{\Omega_{e}}2(\nu+\nu_{T})|\nabla^{S}\overline{{\bf{u}}}^{h}|^{2}d\Omega+\sum_{e}\sum_{i}\int_{\Gamma_{e_{i}}}\frac{2C_{pen}}{h_{e}}(\nu+\nu_{T})|\hat{\overline{{\bf{u}}}}^{h}-\overline{{\bf{u}}}^{h}|^{2}d\Gamma$
	$\displaystyle+4\sum_{e}\sum_{i}$	$\displaystyle\int_{\Gamma_{e_{i}}}(\nu+\nu_{T})(\nabla^{S}\overline{{\bf{u}}}^{h}\cdot{\bf{n}})\cdot(\hat{\overline{{\bf{u}}}}^{h}-\overline{{\bf{u}}}^{h})d\Gamma+\int_{\Gamma_{N}}(1-\lambda)(\hat{\overline{{\bf{u}}}}^{h}\cdot{\bf{n}})|\hat{\overline{{\bf{u}}}}^{h}|^{2}d\Gamma$
	$\displaystyle-\sum_{e}$	$\displaystyle\int_{\Omega_{e}}(\overline{{\bf{u}}}^{h}\otimes\overline{{\bf{u}}}^{h}):\nabla\overline{{\bf{u}}}^{h}d\Omega=0$

Since $\hat{\overline{{\bf{u}}}}^{h}$ is continuous across facets, $\llbracket\overline{{\bf{u}}}^{h}\rrbracket=0$ on interior facets, and $\hat{\overline{{\bf{u}}}}^{h}\cdot{\bf{n}}=\overline{{\bf{u}}}^{h}\cdot{\bf{n}}$ on the domain boundary, the third integral on the left hand side of Eq.(26) can be simplified as follows:

$$-\sum_{e}\sum_{i}\frac{1}{2}\int_{\Gamma_{e_{i}}}(\overline{{\bf{u}}}^{h}\cdot{\bf{n}})|\hat{\overline{{\bf{u}}}}^{h}|^{2}d\Gamma=-\frac{1}{2}\int_{\Gamma_{N}}(\hat{\overline{{\bf{u}}}}^{h}\cdot{\bf{n}})|\hat{\overline{{\bf{u}}}}^{h}|^{2}d\Gamma$$

(27)

By the product rule, it holds that $-\overline{{\bf{u}}}^{h}\otimes\overline{{\bf{u}}}^{h}:\nabla\overline{{\bf{u}}}^{h}=-\nabla\cdot((\overline{{\bf{u}}}^{h}\otimes\overline{{\bf{u}}}^{h})\cdot\overline{{\bf{u}}}^{h})/2$ since $\nabla\cdot\overline{{\bf{u}}}^{h}=0$ on each element $\Omega_{e}$. By the divergence theorem, it then holds that the last term on the left hand side of Eq. (26) is equal to:

$$-\sum_{e}\int_{\Omega_{e}}(\overline{{\bf{u}}}^{h}\otimes\overline{{\bf{u}}}^{h}):\nabla\overline{{\bf{u}}}^{h}d\Omega=-\frac{1}{2}\sum_{e}\sum_{i}\int_{\Gamma_{e_{i}}}(\overline{{\bf{u}}}^{h}\cdot{\bf{n}})|\overline{{\bf{u}}}^{h}|^{2}d\Gamma$$

(28)

Combining Eqs. (26-28) yields:

	$\displaystyle\sum_{e}\frac{1}{2}\int_{\Omega_{e}}\frac{\partial|\overline{{\bf{u}}}^{h}|^{2}}{\partial t}d\Omega=-\sum_{e}\sum_{i}\frac{1}{2}\int_{\Gamma_{e_{i}}}|\overline{{\bf{u}}}^{h}\cdot{\bf{n}}||\overline{{\bf{u}}}^{h}-\hat{\overline{{\bf{u}}}}^{h}|^{2}d\Gamma$		(29)
	$\displaystyle-\sum_{e}\int_{\Omega_{e}}2(\nu+\nu_{T})|\nabla^{S}\overline{{\bf{u}}}^{h}|^{2}d\Omega-\sum_{e}\sum_{i}\int_{\Gamma_{e_{i}}}\frac{2C_{pen}}{h_{e}}(\nu+\nu_{T})|\hat{\overline{{\bf{u}}}}^{h}-\overline{{\bf{u}}}^{h}|^{2}d\Gamma$
	$\displaystyle-4\sum_{e}\sum_{i}\int_{\Gamma_{e_{i}}}(\nu+\nu_{T})(\nabla^{S}\overline{{\bf{u}}}^{h}\cdot{\bf{n}})\cdot(\hat{\overline{{\bf{u}}}}^{h}-\overline{{\bf{u}}}^{h})d\Gamma-\frac{1}{2}\sum_{e}\sum_{i}\int_{\Gamma_{e_{i}}}(\hat{\overline{{\bf{u}}}}^{h}\cdot{\bf{n}})|\hat{\overline{{\bf{u}}}}^{h}|^{2}d\Gamma$

where we have used that:

$$\int_{\Gamma_{N}}(1-\lambda)(\hat{\overline{{\bf{u}}}}^{h}\cdot{\bf{n}})|\hat{\overline{{\bf{u}}}}^{h}|^{2}d\Gamma-\frac{1}{2}\int_{\Gamma_{N}}(\hat{\overline{{\bf{u}}}}^{h}\cdot{\bf{n}})|\hat{\overline{{\bf{u}}}}^{h}|^{2}d\Gamma=\frac{1}{2}\sum_{e}\sum_{i}\int_{\Gamma_{e_{i}}}(\hat{\overline{{\bf{u}}}}^{h}\cdot{\bf{n}})|\hat{\overline{{\bf{u}}}}^{h}|^{2}d\Gamma$$

(30)

By the Cauchy-Schwarz inequality and Young’s inequality, it holds that:

		$\displaystyle-4\sum_{e}\sum_{i}\int_{\Gamma_{e_{i}}}(\nu+\nu_{T})(\nabla^{S}\overline{{\bf{u}}}^{h}\cdot{\bf{n}})\cdot(\hat{\overline{{\bf{u}}}}^{h}-\overline{{\bf{u}}}^{h})d\Gamma\leq$		(31)
		$\displaystyle\sum_{e}\int_{\partial\Omega_{e}}\frac{2h_{e}}{C_{pen}}(\nu+\nu_{T})|\nabla^{S}\overline{{\bf{u}}}^{h}|^{2}d\Gamma+\sum_{e}\sum_{i}\int_{\Gamma_{e_{i}}}\frac{2C_{pen}}{h_{e}}(\nu+\nu_{T})|\hat{\overline{{\bf{u}}}}^{h}-\overline{{\bf{u}}}^{h}|^{2}d\Gamma$

and by the trace inequality, it holds that

$$\int_{\partial\Omega_{e}}(\nu+\nu_{T})|\nabla^{S}\overline{{\bf{u}}}^{h}|^{2}d\Gamma\leq\frac{(k+1)(k+2)}{h_{e}}\int_{\Omega_{e}}(\nu+\nu_{T})|\nabla^{S}\overline{{\bf{u}}}^{h}|^{2}d\Omega$$

(32)

for each element $\Omega_{e}$. Combining Eqs. (29-32) and exploiting the fact that $C_{pen}\geq(k+1)(k+2)$ yields:

$\displaystyle\sum_{e}\frac{1}{2}\int_{\Omega_{e}}\frac{\partial|\overline{{\bf{u}}}^{h}|^{2}}{\partial t}d\Omega\leq-\sum_{e}\sum_{i}\frac{1}{2}\int_{\Gamma_{e_{i}}}|\overline{{\bf{u}}}^{h}\cdot{\bf{n}}||\overline{{\bf{u}}}^{h}-\hat{\overline{{\bf{u}}}}^{h}|^{2}d\Gamma-\frac{1}{2}\sum_{e}\sum_{i}\int_{\Gamma_{e_{i}}}(\hat{\overline{{\bf{u}}}}^{h}\cdot{\bf{n}})|\hat{\overline{{\bf{u}}}}^{h}|^{2}d\Gamma\leq 0$

(33)

This completes the proof. ∎