High-order DG solvers for under-resolved turbulent incompressible flows: A comparison of $L^2$ and $H$(div) methods

TL;DR

This paper compares high-order DG methods based on $L^2$ and $H(div)$ approaches for simulating turbulent incompressible flows, focusing on their mass conservation, energy stability, and performance in under-resolved regimes.

Contribution

It provides a comparative analysis of $L^2$ and $H(div)$ DG discretisations, highlighting their similarities, differences, and suitability for turbulent flow simulations.

Findings

Both methods are promising for under-resolved turbulence simulations.

$H(div)$ methods ensure exact mass conservation and energy stability.

$L^2$ methods achieve similar robustness with stabilization techniques.

Abstract

The accurate numerical simulation of turbulent incompressible flows is a challenging topic in computational fluid dynamics. For discretisation methods to be robust in the under-resolved regime, mass conservation as well as energy stability are key ingredients to obtain robust and accurate discretisations. Recently, two approaches have been proposed in the context of high-order discontinuous Galerkin (DG) discretisations that address these aspects differently. On the one hand, standard $L^2$-based DG discretisations enforce mass conservation and energy stability weakly by the use of additional stabilisation terms. On the other hand, pointwise divergence-free $H(\operatorname{div})$-conforming approaches ensure exact mass conservation and energy stability by the use of tailored finite element function spaces. The present work raises the question whether and to which extent these two approaches are equivalent when applied to under-resolved turbulent flows. This comparative study highlights similarities and differences of these two approaches. The numerical results emphasise that both discretisation strategies are promising for under-resolved simulations of turbulent flows due to their inherent dissipation mechanisms.

Figures

Equations

$${\partial }_t\mathbf{u}-\nu \Delta \mathbf{u}+\left(\mathbf{u}\text{vbox}\text{scalebox}0.635\bullet \nabla \right)\mathbf{u}+\nabla p$$

$$\nabla \text{vbox}\text{scalebox}0.635\bullet \mathbf{u}$$

$$\mathbf{u}\left(0,\cdot \right)$$

$$\displaystyle{\left\llbracket\phi\right\rrbracket}_{F}=\phi^{+}-\phi^{-}\quad\text{and}\quad{{\big{\{}\big{\{}}\phi{\big{\}}\big{\}}}}_{F}=\frac{1}{2}{\left(\phi^{+}+\phi^{-}\right)}.$$

$${\mathbf{H}}^m\left({\mathcal{T}}_h\right)=\left\{\mathbf{v}\in {\mathbf{L}}^2\left(\Omega \right):{\mathbf{v}|}_K\in {\mathbf{H}}^m\left(K\right),\forall K\in {\mathcal{T}}_h\right\},$$

$${\mathbb{Q}_{k}{}}\coloneqq\bigoplus_{K\in{\mathcal{T}_{h}}}{\mathbb{Q}_{k}{(K)}},\quad{\ThisStyle{\ooalign{$\SavedStyle\mathbb{Q}$\cr\kern-0.22pt$\SavedStyle\mathbb{Q}$\cr\kern 0.22pt$\SavedStyle\mathbb{Q}$}}_{k}{}}\coloneqq\bigoplus_{K\in{\mathcal{T}_{h}}}{\ThisStyle{\ooalign{$\SavedStyle\mathbb{Q}$\cr\kern-0.22pt$\SavedStyle\mathbb{Q}$\cr\kern 0.22pt$\SavedStyle\mathbb{Q}$}}_{k}{(K)}}.$$

$${\mathbf{H}}_0\left(\mathrm{div};\Omega \right)=\left\{\mathbf{v}\in {\mathbf{L}}^2\left(\Omega \right):\nabla \text{vbox}\text{scalebox}0.635\bullet \mathbf{v}\in L^2\left(\Omega \right);{\mathbf{v}\text{vbox}\text{scalebox}0.635\bullet \mathbf{n}|}_{{\Gamma }_0}=0\right\}.$$

$${\ThisStyle{\ooalign{$\SavedStyle\mathbb{RT}$\cr\kern-0.22pt$\SavedStyle\mathbb{RT}$\cr\kern 0.22pt$\SavedStyle\mathbb{RT}$}}_{{\left[k\right]}}{{\left(K\right)}}}={\left({\mathbb{P}_{k+1,k,k}{{\left(K\right)}}},{\mathbb{P}_{k,k+1,k}{{\left(K\right)}}},{\mathbb{P}_{k,k,k+1}{{\left(K\right)}}}\right)}^{\dagger},\quad{\mathbb{P}_{k_{1},k_{2},k_{3}}{{\left(K\right)}}}={\mathbb{P}_{k_{1}}{(K)}}\otimes{\mathbb{P}_{k_{2}}{(K)}}\otimes{\mathbb{P}_{k_{3}}{(K)}}.$$

$${\ThisStyle{\ooalign{$\SavedStyle\mathbb{RT}$\cr\kern-0.22pt$\SavedStyle\mathbb{RT}$\cr\kern 0.22pt$\SavedStyle\mathbb{RT}$}}_{{\left[k\right]}}{}}\coloneqq{\bm{H}_{0}{{\left({\mathrm{div}};\Omega\right)}}}\cap\bigoplus_{K\in{\mathcal{T}_{h}}}{\ThisStyle{\ooalign{$\SavedStyle\mathbb{RT}$\cr\kern-0.22pt$\SavedStyle\mathbb{RT}$\cr\kern 0.22pt$\SavedStyle\mathbb{RT}$}}_{{\left[k\right]}}{(K)}}.$$

$${\bm{V}}_{h}={\ThisStyle{\ooalign{$\SavedStyle\mathbb{Q}$\cr\kern-0.22pt$\SavedStyle\mathbb{Q}$\cr\kern 0.22pt$\SavedStyle\mathbb{Q}$}}_{k}{}},\quad{Q}_{h}={\mathbb{Q}_{k-1}{}}\cap{L_{0}^{2}{{{\left(\Omega\right)}}}}.$$

$$\text{Find }\left({\mathbf{u}}_h,p_h\right):\left(0,T\right]\to {\mathbf{V}}_h\times Q_h\text{ with }{\mathbf{u}}_h\left(0\right)={\mathbf{u}}_{0h}\text{ s.t., }\forall \left({\mathbf{v}}_h,q_h\right)\in {\mathbf{V}}_h\times Q_h,$$

$$\left({\partial }_t{\mathbf{u}}_h,{\mathbf{v}}_h\right)+\nu a_h\left({\mathbf{u}}_h,{\mathbf{v}}_h\right)+c_h\left({\mathbf{u}}_h;{\mathbf{u}}_h,{\mathbf{v}}_h\right)+j_h\left({\mathbf{u}}_h,{\mathbf{v}}_h\right)+g_h\left(p_h,{\mathbf{v}}_h\right)+d_h\left({\mathbf{u}}_h,q_h\right)=\left(\mathbf{f},{\mathbf{v}}_h\right).$$

$$a_{h}{\left({\bm{u}}_{h},{\bm{v}}_{h}\right)}=\int_{\Omega}\nabla_{h}{\bm{u}}_{h}\operatorname{{\bm{\mathbin{\vbox{\hbox{$\ldotp$}\hbox{$\ldotp$}}}}}}\nabla_{h}{\bm{v}}_{h}{\,{\mathrm{d}}{\bm{x}}}-\sum_{F\in{\mathcal{F}_{h}}}\int_{F}\left({{\big{\{}\big{\{}}\nabla{\bm{u}}_{h}{\big{\}}\big{\}}}}\operatorname{{\mathchoice{\mathbin{\vbox{\hbox{\scalebox{0.635}{$\displaystyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.635}{$\textstyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.635}{$\scriptstyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.635}{$\scriptscriptstyle\bullet$}}}}}}}{\bm{n}}\right)\operatorname{{\mathchoice{\mathbin{\vbox{\hbox{\scalebox{0.635}{$\displaystyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.635}{$\textstyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.635}{$\scriptstyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.635}{$\scriptscriptstyle\bullet$}}}}}}}{\left\llbracket{\bm{v}}_{h}\right\rrbracket}{\,{\mathrm{d}}\bm{s}}-\sum_{F\in{\mathcal{F}_{h}}}\int_{F}{\left\llbracket{\bm{u}}_{h}\right\rrbracket}\operatorname{{\mathchoice{\mathbin{\vbox{\hbox{\scalebox{0.635}{$\displaystyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.635}{$\textstyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.635}{$\scriptstyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.635}{$\scriptscriptstyle\bullet$}}}}}}}\left({{\big{\{}\big{\{}}\nabla{\bm{v}}_{h}{\big{\}}\big{\}}}}\operatorname{{\mathchoice{\mathbin{\vbox{\hbox{\scalebox{0.635}{$\displaystyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.635}{$\textstyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.635}{$\scriptstyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.635}{$\scriptscriptstyle\bullet$}}}}}}}{\bm{n}}\right){\,{\mathrm{d}}\bm{s}}+\sum_{F\in{\mathcal{F}_{h}}}\int_{F}\tau{\left\llbracket{\bm{u}}_{h}\right\rrbracket}\operatorname{{\mathchoice{\mathbin{\vbox{\hbox{\scalebox{0.635}{$\displaystyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.635}{$\textstyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.635}{$\scriptstyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.635}{$\scriptscriptstyle\bullet$}}}}}}}{\left\llbracket{\bm{v}}_{h}\right\rrbracket}{\,{\mathrm{d}}\bm{s}},$$

$$\tau =\{\begin{array}{ll}\max \left({\tau }_{K^{-}},{\tau }_{K^{+}}\right)& \text{if facet }F\in {\mathcal{F}}_h^i,\\ {\tau }_K& \text{if facet }F\in {\mathcal{F}}_h^{\partial },\end{array}$$

$${\tau }_K=(k+1{)}^2\frac{{\left|\partial K\cap {\mathcal{F}}_h^i\right|}_{d-1}/2+{\left|\partial K\cap {\mathcal{F}}_h^{\partial }\right|}_{d-1}}{{\left|K\right|}_d}.$$

$$c_{h}{\left({\bm{w}}_{h};{\bm{u}}_{h},{\bm{v}}_{h}\right)}=-\int_{\Omega}{\left({\bm{u}}_{h}\otimes{\bm{w}}_{h}\right)}\operatorname{{\bm{\mathbin{\vbox{\hbox{$\ldotp$}\hbox{$\ldotp$}}}}}}\nabla_{h}{\bm{v}}_{h}{\,{\mathrm{d}}{\bm{x}}}+\sum_{F\in{\mathcal{F}_{h}}}\int_{F}{\left({{\big{\{}\big{\{}}{\bm{u}}_{h}\otimes{\bm{w}}_{h}{\big{\}}\big{\}}}}\operatorname{{\mathchoice{\mathbin{\vbox{\hbox{\scalebox{0.635}{$\displaystyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.635}{$\textstyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.635}{$\scriptstyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.635}{$\scriptscriptstyle\bullet$}}}}}}}{\bm{n}}+\frac{\Lambda({\bm{w}}_{h})}{2}{\left\llbracket{\bm{u}}_{n}\right\rrbracket}\right)}\operatorname{{\mathchoice{\mathbin{\vbox{\hbox{\scalebox{0.635}{$\displaystyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.635}{$\textstyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.635}{$\scriptstyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.635}{$\scriptscriptstyle\bullet$}}}}}}}{\left\llbracket{\bm{v}}_{h}\right\rrbracket}{\,{\mathrm{d}}\bm{s}},$$

$$g_h\left(p_h,{\mathbf{v}}_h\right)$$

$$d_h\left({\mathbf{u}}_h,q_h\right)$$

$$j_{\mathrm{div},h}\left({\mathbf{u}}_h,{\mathbf{v}}_h\right)$$

$$j_{\mathrm{conti},h}\left({\mathbf{u}}_h,{\mathbf{v}}_h\right)$$

$${\bm{V}}_{h}={\ThisStyle{\ooalign{$\SavedStyle\mathbb{RT}$\cr\kern-0.22pt$\SavedStyle\mathbb{RT}$\cr\kern 0.22pt$\SavedStyle\mathbb{RT}$}}_{{\left[k\right]}}{}},\quad{Q}_{h}={\mathbb{Q}_{k}{}}\cap{L_{0}^{2}{{{\left(\Omega\right)}}}},\quad{\bm{\widehat{V}}}_{h}={\left\{{\bm{\widehat{v}}}_{h}\in{\bm{L}^{2}{{\left({\mathcal{F}_{h}}\right)}}}\colon{\left.\kern-1.2pt{\bm{\widehat{v}}}_{h}\vphantom{|}\right|_{F}}\in{\ThisStyle{\ooalign{$\SavedStyle\mathbb{Q}$\cr\kern-0.22pt$\SavedStyle\mathbb{Q}$\cr\kern 0.22pt$\SavedStyle\mathbb{Q}$}}_{k-1}{{\left(F\right)}}};~{}{\left.\kern-1.2pt{\bm{\widehat{v}}}_{h}\operatorname{{\mathchoice{\mathbin{\vbox{\hbox{\scalebox{0.635}{$\displaystyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.635}{$\textstyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.635}{$\scriptstyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.635}{$\scriptscriptstyle\bullet$}}}}}}}{\bm{n}}\vphantom{|}\right|_{F}}=0,~{}\forall\,F\in{\mathcal{F}_{h}};~{}{\left.\kern-1.2pt{\bm{\widehat{v}}}_{h}\vphantom{|}\right|_{\Gamma_{0}}}={\bm{0}}\right\}},$$

$$\text{Find }\left({\mathbf{u}}_h,{\hat{\mathbf{u}}}_h,p_h\right):\left(0,T\right]\to {\mathbf{V}}_h\times {\hat{\mathbf{V}}}_h\times Q_h\text{ with }{\mathbf{u}}_h\left(0\right)={\mathbf{u}}_{0h}\text{ s.t., }\forall \left({\mathbf{v}}_h,{\hat{\mathbf{v}}}_h,q_h\right)\in {\mathbf{V}}_h\times {\hat{\mathbf{V}}}_h\times Q_h,$$

$$\left({\partial }_t{\mathbf{u}}_h,{\mathbf{v}}_h\right)+\nu a_h\left(\left({\mathbf{u}}_h,{\hat{\mathbf{u}}}_h\right),\left({\mathbf{v}}_h,{\hat{\mathbf{v}}}_h\right)\right)+c_h\left({\mathbf{u}}_h;{\mathbf{u}}_h,{\mathbf{v}}_h\right)+b_h\left({\mathbf{v}}_h,p_h\right)+b_h\left({\mathbf{u}}_h,q_h\right)=\left(\mathbf{f},{\mathbf{v}}_h\right).$$

$$a_h\left(\left({\mathbf{u}}_h,{\hat{\mathbf{u}}}_h\right),\left({\mathbf{v}}_h,{\hat{\mathbf{v}}}_h\right)\right)=\sum _{K\in {\mathcal{T}}_h}\left[{\int }_K\nabla {\mathbf{u}}_h\text{vbox}..\nabla {\mathbf{v}}_h\mathrm{d}\mathbf{x}-{\int }_{\partial K}\left(\nabla {\mathbf{u}}_h\text{vbox}\text{scalebox}0.635\bullet \mathbf{n}\right)\text{vbox}\text{scalebox}0.635\bullet [[[[{\mathbf{v}}_h]]]]\mathrm{d}\mathbf{s}-{\int }_{\partial K}[[[[{\mathbf{u}}_h]]]]\text{vbox}\text{scalebox}0.635\bullet \left(\nabla {\mathbf{v}}_h\text{vbox}\text{scalebox}0.635\bullet \mathbf{n}\right)\mathrm{d}\mathbf{s}+{\int }_{\partial K}\frac{{\sigma }_K}{h_K\left(\mathbf{s}\right)}[[[[{\mathbf{u}}_h]]]]\text{vbox}\text{scalebox}0.635\bullet [[[[{\mathbf{v}}_h]]]]\mathrm{d}\mathbf{s}\right],$$

$$[[[[\mathbf{v}]]]]={\left({\mathbf{\Pi }}_{{\mathcal{F}}_h}^{k-1}\mathbf{v}-\hat{\mathbf{v}}\right)}_{\mathbf{t}}={\mathbf{\Pi }}_{{\mathcal{F}}_h}^{k-1}{\mathbf{v}}_{\mathbf{t}}-\hat{\mathbf{v}},$$

$$b_h\left({\mathbf{u}}_h,q_h\right)=-{\int }_{\Omega }{q}_h\left(\nabla \text{vbox}\text{scalebox}0.635\bullet {\mathbf{u}}_h\right)\mathrm{d}\mathbf{x}.$$

$$-{\int }_{\Omega }{q}_h\left(\nabla \text{vbox}\text{scalebox}0.635\bullet {\mathbf{u}}_h\right)\mathrm{d}\mathbf{x}=0,\forall q_h\in Q_h⟹{\int }_{\Omega }{\left|\nabla \text{vbox}\text{scalebox}0.635\bullet {\mathbf{u}}_h\right|}^2\mathrm{d}\mathbf{x}=0⟹\nabla \text{vbox}\text{scalebox}0.635\bullet {\mathbf{u}}_h=0.$$

$$c_{h}{\left({\bm{w}}_{h};{\bm{u}}_{h},{\bm{v}}_{h}\right)}=\int_{\Omega}{\left({\bm{w}}_{h}\operatorname{{\mathchoice{\mathbin{\vbox{\hbox{\scalebox{0.635}{$\displaystyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.635}{$\textstyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.635}{$\scriptstyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.635}{$\scriptscriptstyle\bullet$}}}}}}}\nabla_{h}\right)}{\bm{u}}_{h}\operatorname{{\mathchoice{\mathbin{\vbox{\hbox{\scalebox{0.635}{$\displaystyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.635}{$\textstyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.635}{$\scriptstyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.635}{$\scriptscriptstyle\bullet$}}}}}}}{\bm{v}}_{h}{\,{\mathrm{d}}{\bm{x}}}-\sum_{F\in{\mathcal{F}_{h}^{i}}}\int_{F}{\left({\bm{w}}_{h}\operatorname{{\mathchoice{\mathbin{\vbox{\hbox{\scalebox{0.635}{$\displaystyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.635}{$\textstyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.635}{$\scriptstyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.635}{$\scriptscriptstyle\bullet$}}}}}}}{\bm{n}}\right)}{\left\llbracket{\bm{u}}_{h}\right\rrbracket}\operatorname{{\mathchoice{\mathbin{\vbox{\hbox{\scalebox{0.635}{$\displaystyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.635}{$\textstyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.635}{$\scriptstyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.635}{$\scriptscriptstyle\bullet$}}}}}}}{{\big{\{}\big{\{}}{\bm{v}}_{h}{\big{\}}\big{\}}}}{\,{\mathrm{d}}\bm{s}}+\sum_{F\in{\mathcal{F}_{h}^{i}}}\int_{F}\frac{\theta}{2}{\left\lvert{\bm{w}}_{h}\operatorname{{\mathchoice{\mathbin{\vbox{\hbox{\scalebox{0.635}{$\displaystyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.635}{$\textstyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.635}{$\scriptstyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.635}{$\scriptscriptstyle\bullet$}}}}}}}{\bm{n}}\right\rvert}{\left\llbracket{\bm{u}}_{h}\right\rrbracket}\operatorname{{\mathchoice{\mathbin{\vbox{\hbox{\scalebox{0.635}{$\displaystyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.635}{$\textstyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.635}{$\scriptstyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.635}{$\scriptscriptstyle\bullet$}}}}}}}{\left\llbracket{\bm{v}}_{h}\right\rrbracket}{\,{\mathrm{d}}\bm{s}},$$

$$N_{\mathrm{DOFs},\mathbf{u}}(N,k)=N^{d}d{\left(k+1\right)}^d.$$

$$N_{\mathrm{DOFs},\mathbf{u}}^{\text{L2,c}}(N,k)=N^d\left((d-1){\left(k+1\right)}^d-d(k+1{)}^{d-1}+1\right),$$

$$N_{\mathrm{DOFs},\mathbf{u}}^{\text{L2,c}}(N,k)<N_{\mathrm{DOFs},\mathbf{u}}^{\text{Hdiv,c}}(N,k)=N^d(d-1){\left(k+1\right)}^d<N_{\mathrm{DOFs},\mathbf{u}}^{\text{L2,c}}(N,k+1).$$

$$d_h\left({\mathbf{u}}_h,q_h\right)$$

$$\underset{\text{dt\_ekin}}{\underbrace{\frac{1}{2}\frac{\mathrm{d}}{\mathrm{d}t}{\Vert {\mathbf{u}}_h\Vert }_{{\mathbf{L}}^2}^2}}+\underset{\text{visc\_diss}}{\underbrace{\nu a_h\left({\mathbf{u}}_h,{\mathbf{u}}_h\right)}}+\underset{\text{num\_diss\_conv}}{\underbrace{c_h\left({\mathbf{u}}_h;{\mathbf{u}}_h,{\mathbf{u}}_h\right)}}+\underset{\text{num\_diss\_div}+\text{num\_diss\_conti}}{\underbrace{j_h\left({\mathbf{u}}_h,{\mathbf{u}}_h\right)}}=0.$$

Full text

††articletype: Research article

1 Introduction

1.1 Motivation

**Simulating turbulent flows is still a challenging undertaking, even on today’s high-performance computing architectures. Discontinuous Galerkin (DG) discretisations are currently investigated in order to develop new discretisation methods with inbuilt stabilisation mechanism rendering these methods robust and accurate when applied to turbulent flow problems. In this contribution, we compare the accuracy of two high-order DG solvers for incompressible flows with a special emphasis on how they perform in the practically relevant situation of only being able to marginally resolve the occurring flow features. While the two ‘extreme cases’ of direct numerical simulation (DNS) on the one hand, and Reynolds-averaged Navier–Stokes (RANS) simulations on the other hand are relatively well-established and well-understood, being able to perform time-dependent turbulent flow simulations with only a limited amount of fine-scale accuracy usually goes by the name large-eddy simulation (LES). This regime is exactly where we are interested in in this work.  **

**LES simulations by now have a rich history in both the engineering and the mathematics literature and there are numerous promising approaches for different flow situations. Our aim is to be able to compute (incompressible) flows without the need of choosing the ‘correct’ turbulence model or set of parameters. Therefore, we rely on what can be called a ‘no-model LES’ in the sense that no explicit turbulence model shall be incorporated in the simulations in this work. Such an approach, in turn, relies heavily on distinctive robustness as well as healthy and controllable dissipation properties of the numerical scheme to approximate the incompressible Navier–Stokes equations. While DG methods have reached a mature state in the field of the compressible Navier–Stokes equations, see for example the works 1, 2, 3, 4, 5, 6, 7 in the context of low- and moderate-Mach number turbulent flows, designing robust DG discretisations of the incompressible Navier–Stokes equations exhibits some subtleties in the context of under-resolved turbulence: applying well-known numerical flux formulations to the discretisation of convective and viscous terms only, see the works 8, 9, 10, 11, 12, 13, has been realised to be not robust enough in under-resolved scenarios, as pointed out in the recent works 14, 15, 16, 17. It is worth mentioning that this lack of robustness is not related to under-integration of nonlinear terms, commonly known as aliasing. Instead, additional techniques are required which are inherently linked to the nature of the incompressible Navier–Stokes equations and in particular the incompressibility constraint. Namely, effort has to be put into ensuring compliance with both the divergence-free and the normal-continuity constraints.  **

Two promising methods have been proposed recently in 18, 19 which rely on different stabilisation concepts. The two candidates under investigation are intelligently stabilised ${\bm{L}^{2}}$-based DG methods 19 and exactly divergence-free ${\bm{H}{{\left({\mathrm{div}}\right)}}}$-conforming methods 18, 20. A first connection between both approaches has been pointed out in 19, 21. However, it remains rather unexplored how these different stabilisation techniques behave regarding their dissipation mechanism in the under-resolved regime, which is particularly relevant for practical LES of incompressible turbulent flows. Therefore, the aim of this work is a comparison of the accuracy of both approaches in under-resolved turbulence simulation, which is examined numerically for the well-known 3D benchmark cases of freely decaying turbulence (Taylor–Green vortex), and wall-bounded turbulence in a channel. The numerical results for the ${\bm{H}{{\left({\mathrm{div}}\right)}}}$-conforming method have been obtained using NGSolve22, the results for the ${\bm{L}^{2}}$-based DG using the open-source finite element library deal.II23.

1.2 Mathematical model

The underlying mathematical model is the incompressible Navier–Stokes problem24, 25, 26** without any additional terms, i.e.**

[TABLE]

In (1), $\Omega\subset\mathbb{R}^{3}$ denotes a connected bounded Lipschitz domain. Concerning boundary conditions (BCs) we assume that $\partial\Omega=\Gamma_{0}\,\dot{\cup}\,\Gamma_{\mathrm{per}}$, where the no-slip BC ${\bm{u}}={\bm{0}}$ is prescribed on $\Gamma_{0}$ and periodic BCs are imposed on $\Gamma_{\mathrm{per}}$. Moreover, ${\bm{u}}\colon{\left(0,{T}\right]}\times\Omega\to\mathbb{R}^{3}$ indicates the velocity field, $p\colon{\left(0,{T}\right]}\times\Omega\to\mathbb{R}$ is the (zero-mean) kinematic pressure, ${\bm{f}}\colon{\left(0,{T}\right]}\times\Omega\to\mathbb{R}^{3}$ represents external body forces and ${\bm{u}}_{0}\colon\Omega\to\mathbb{R}^{3}$ stands for a suitable initial condition for the velocity. The underlying fluid is assumed to be Newtonian with constant (dimensionless) kinematic viscosity $0<\nu\ll 1$.

1.3 Outline

We introduce both DG discretisations under consideration with a focus on the mathematical formulation and discuss properties such as approximation abilities, mass conservation and dissipation mechanisms in Section 2. The accuracy of both methods for the 3D Taylor-Green vortex problem as a prototype problem of decaying homogeneous isotropic turbulence is investigated in Section 3. Apart from a direct comparison in terms of the kinetic energy dissipation rate on different meshes, we also investigate the impact of low-order vs. high-order discretisations and the impact of variations of the DG formulation such as changes in the stabilisation terms, and the flux formulations for convective and viscous terms. To include wall-bounded turbulent flows in our investigations, we consider the turbulent channel flow problem in Section 4. In that section, we compare the behaviour of both methods under mesh refinement and again discuss the influence of method variations including the handling of anisotropic refinements towards the wall boundary with affine and isoparametric finite element meshes. The robustness of the considered approaches and their reliability in correctly predicting turbulent flows is further assessed by parameter studies. The final part, Section 5, summarises the observations made in this work and derives conclusions.

2 Two DG-based solvers for turbulent incompressible flows

2.1 Notation

Let ${\mathcal{T}_{h}}$ be a partition of $\Omega$ into hexahedra with mesh size $h=\max_{K\in{\mathcal{T}_{h}}}h_{K}$, where $h_{K}$ denotes the diameter of the particular element $K\in{\mathcal{T}_{h}}$. The skeleton ${\mathcal{F}_{h}}$ denotes the set of all facets and ${\mathcal{F}_{h}}={\mathcal{F}_{h}^{i}}\cup{\mathcal{F}_{h}^{\partial}}$ where ${\mathcal{F}_{h}^{i}}$ is the subset of interior plus periodic facets and ${\mathcal{F}_{h}^{\partial}}$ collects all Dirichlet boundary facets $F\subset\Gamma_{0}$. To any $F\in{\mathcal{F}_{h}}$ we assign a unit normal vector ${\bm{n}}_{F}$ where, for $F\subset\partial\Omega$, this is the outer unit normal vector ${\bm{n}}$. If $F\in{\mathcal{F}_{h}^{i}}$, there are two adjacent elements $K^{+}$ and $K^{-}$ sharing the facet $F=\partial K^{+}\cap\partial K^{-}$ and ${\bm{n}}_{F}$ points from $K^{+}$ to $K^{-}$. Let $\phi$ be any piecewise smooth (scalar-, vector- or matrix-valued) function with traces from within the interior of $K^{\pm}$ denoted by $\phi^{\pm}$, respectively. Then, we define the jump {\left\llbracket\cdot\right\rrbracket}_{F} and average {{\big{\{}\big{\{}}\cdot{\big{\}}\big{\}}}}_{F} operators across interior and periodic facets $F\in{\mathcal{F}_{h}^{i}}$ by

[TABLE]

**These operators act componentwise for vector- and matrix-valued functions. Frequently, the subscript indicating the particular facet is omitted.  **

As is usual in DG methods, we make use of the broken vector-valued Sobolev space ($m\in\mathbb{N}_{0}$)

[TABLE]

**and define the (elementwise) broken Jacobian $\nabla_{h}\colon{\bm{H}^{1}{{\left({\mathcal{T}_{h}}\right)}}}\to{\bm{L}^{2}{{{\left(\Omega\right)}}}}$ by ${\left.\kern-1.2pt{\left(\nabla_{h}{\bm{v}}\right)}\vphantom{|}\right|_{K}}=\nabla{\left({\left.\kern-1.2pt{\bm{v}}\vphantom{|}\right|_{K}}\right)}$ and the (elementwise) broken divergence \operatorname{\nabla_{\textit{h}}\operatorname{{\mathchoice{\mathbin{\vbox{\hbox{\scalebox{0.635}{\displaystyle\bullet}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.635}{\textstyle\bullet}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.635}{\scriptstyle\bullet}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.635}{\scriptscriptstyle\bullet}}}}}}}}\colon{\bm{H}^{1}{{\left({\mathcal{T}_{h}}\right)}}}\to{L^{2}{{{\left(\Omega\right)}}}} by {\left.\kern-1.2pt{\left(\operatorname{\nabla_{\textit{h}}\operatorname{{\mathchoice{\mathbin{\vbox{\hbox{\scalebox{0.635}{\displaystyle\bullet}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.635}{\textstyle\bullet}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.635}{\scriptstyle\bullet}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.635}{\scriptscriptstyle\bullet}}}}}}}}{\bm{v}}\right)}\vphantom{|}\right|_{K}}=\operatorname{\nabla\operatorname{{\mathchoice{\mathbin{\vbox{\hbox{\scalebox{0.635}{\displaystyle\bullet}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.635}{\textstyle\bullet}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.635}{\scriptstyle\bullet}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.635}{\scriptscriptstyle\bullet}}}}}}}}{\left({\left.\kern-1.2pt{\bm{v}}\vphantom{|}\right|_{K}}\right)} for all $K\in{\mathcal{T}_{h}}$.  **

**In the remainder of this subsection we introduce local and global finite element spaces. For ease of presentation we consider only hexahedral meshes which consist of straight hexahedra $K$ which are aligned with the coordinate axes and are the image of the reference hexahedron $\hat{K}=[0,1]^{3}$ under affine linear mappings. However, the general case of non-affinely mapped hexadra can also be dealt with easily, cf. the discussion in Section 4.1 below.  **

In the following, ${\mathbb{Q}_{k}{{\left(K\right)}}}$ (vector-valued {\ThisStyle{\ooalign{\SavedStyle\mathbb{Q}$\cr\kern-0.22pt$\SavedStyle\mathbb{Q}$\cr\kern 0.22pt$\SavedStyle\mathbb{Q}}}_{k}{{\left(K\right)}}}) denotes the tensor-product space of all polynomials on $K$ with degree less or equal $k$ in each direction separately, ${\mathbb{Q}_{k}{{\left(K\right)}}}={\mathbb{P}_{k}{(K)}}\otimes{\mathbb{P}_{k}{(K)}}\otimes{\mathbb{P}_{k}{(K)}}$. The ${\bm{L}^{2}}$-based method makes use of these local spaces by glueing them together discontinuously, resulting in the global spaces

[TABLE]

For the ${\bm{H}{{\left({\mathrm{div}}\right)}}}$-based method, let us first introduce

[TABLE]

Note that a function in ${\bm{H}_{0}{{\left({\mathrm{div}};\Omega\right)}}}$ at least has to be normal-continuous. Since we are working on hexahedral meshes in 3D in this work, the local velocity space is based on the Raviart–Thomas element 27

[TABLE]

{\ThisStyle{\ooalign{\SavedStyle\mathbb{RT}$\cr\kern-0.22pt$\SavedStyle\mathbb{RT}$\cr\kern 0.22pt$\SavedStyle\mathbb{RT}}}_{{\left[k\right]}}{(}}K)** includes polynomials of order $k+1$ in certain directions (depending on the vector entry), but much less than {\ThisStyle{\ooalign{\SavedStyle\mathbb{Q}$\cr\kern-0.22pt$\SavedStyle\mathbb{Q}$\cr\kern 0.22pt$\SavedStyle\mathbb{Q}}}_{k+1}{(}}K). There holds the inclusion {\ThisStyle{\ooalign{\SavedStyle\mathbb{Q}$\cr\kern-0.22pt$\SavedStyle\mathbb{Q}$\cr\kern 0.22pt$\SavedStyle\mathbb{Q}}}_{k}{(K)}}\subset{\ThisStyle{\ooalign{\SavedStyle\mathbb{RT}$\cr\kern-0.22pt$\SavedStyle\mathbb{RT}$\cr\kern 0.22pt$\SavedStyle\mathbb{RT}}}_{{\left[k\right]}}{(K)}}\subset{\ThisStyle{\ooalign{\SavedStyle\mathbb{Q}$\cr\kern-0.22pt$\SavedStyle\mathbb{Q}$\cr\kern 0.22pt$\SavedStyle\mathbb{Q}}}_{k+1}{(K)}}. We further note that the highest order polynomials in {\ThisStyle{\ooalign{\SavedStyle\mathbb{RT}$\cr\kern-0.22pt$\SavedStyle\mathbb{RT}$\cr\kern 0.22pt$\SavedStyle\mathbb{RT}}}_{{\left[k\right]}}{(}}K) do not improve the approximation order of the local polynomial space, but are added on top of {\ThisStyle{\ooalign{\SavedStyle\mathbb{Q}$\cr\kern-0.22pt$\SavedStyle\mathbb{Q}$\cr\kern 0.22pt$\SavedStyle\mathbb{Q}}}_{k}{(K)}} for the purpose of improving the approximation of the divergence operator as there holds \operatorname{\nabla\operatorname{{\mathchoice{\mathbin{\vbox{\hbox{\scalebox{0.635}{\displaystyle\bullet}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.635}{\textstyle\bullet}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.635}{\scriptstyle\bullet}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.635}{\scriptscriptstyle\bullet}}}}}}}}{\ThisStyle{\ooalign{\SavedStyle\mathbb{RT}$\cr\kern-0.22pt$\SavedStyle\mathbb{RT}$\cr\kern 0.22pt$\SavedStyle\mathbb{RT}}}_{{\left[k\right]}}{(K)}}={\mathbb{Q}_{k}{(K)}}\supset\operatorname{\nabla\operatorname{{\mathchoice{\mathbin{\vbox{\hbox{\scalebox{0.635}{\displaystyle\bullet}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.635}{\textstyle\bullet}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.635}{\scriptstyle\bullet}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.635}{\scriptscriptstyle\bullet}}}}}}}}{\ThisStyle{\ooalign{\SavedStyle\mathbb{Q}$\cr\kern-0.22pt$\SavedStyle\mathbb{Q}$\cr\kern 0.22pt$\SavedStyle\mathbb{Q}}}_{k}{(K)}}. We define the global finite element spaces for the ${\bm{H}{{\left({\mathrm{div}}\right)}}}$-based method as**

[TABLE]

2.2 L2-DG method with consistent stabilisation terms

The first DG method investigated in this work is an ${\bm{L}^{2}}$-based DG incompressible Navier–Stokes solver characterised by consistent divergence and continuity penalty terms acting as stabilisation terms for improved mass conservation and energy stability. This method is based on the formulation proposed in recent works of Fehn et al. 19, 28. Both velocity and pressure are approximated by discontinuous functions, where the polynomial space for the velocity is one degree higher than for the pressure for reasons of discrete inf-sup stability:

[TABLE]

To obtain the contribution of boundary face integrals from the weak forms given below, equation (2) holds where the exterior weighting function is simply set to zero on boundary facets (${\bm{v}}_{h}^{+}=\bm{0}$, \nabla{\bm{v}}_{h}^{+}\operatorname{{\mathchoice{\mathbin{\vbox{\hbox{\scalebox{0.635}{\displaystyle\bullet}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.635}{\textstyle\bullet}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.635}{\scriptstyle\bullet}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.635}{\scriptscriptstyle\bullet}}}}}}}{\bm{n}}=\bm{0},q_{h}^{+}=0). Homogeneous Dirichlet boundary conditions are prescribed according to the mirror principle ${\bm{u}}_{h}^{+}=-{\bm{u}}_{h}^{-}$ (with \nabla{\bm{u}}_{h}^{+}\operatorname{{\mathchoice{\mathbin{\vbox{\hbox{\scalebox{0.635}{\displaystyle\bullet}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.635}{\textstyle\bullet}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.635}{\scriptstyle\bullet}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.635}{\scriptscriptstyle\bullet}}}}}}}{\bm{n}}=\nabla{\bm{u}}_{h}^{-}\operatorname{{\mathchoice{\mathbin{\vbox{\hbox{\scalebox{0.635}{\displaystyle\bullet}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.635}{\textstyle\bullet}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.635}{\scriptstyle\bullet}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.635}{\scriptscriptstyle\bullet}}}}}}}{\bm{n}} for the velocity gradient and $p_{h}^{+}=p_{h}^{-}$ for the pressure). The space-semidiscrete weak formulation of the incompressible Navier–Stokes equations (1) reads as follows:

[TABLE]

The DG discretisation of the viscous term is based on the symmetric interior penalty Galerkin (SIPG) method

[TABLE]

where the penalty parameter $\tau$ is

[TABLE]

using the definition according to Hillewaert 29 for hexahedral elements

[TABLE]

Here, $|\cdot|_{d}$ denotes the $d$-dimensional Lebesgue measure. For the convective term, the local Lax–Friedrichs flux is applied

[TABLE]

where \Lambda({\bm{w}}_{h})=\max{\left(2|{\bm{w}}^{-}_{h}\operatorname{{\mathchoice{\mathbin{\vbox{\hbox{\scalebox{0.635}{\displaystyle\bullet}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.635}{\textstyle\bullet}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.635}{\scriptstyle\bullet}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.635}{\scriptscriptstyle\bullet}}}}}}}{\bm{n}}|,2|{\bm{w}}^{+}_{h}\operatorname{{\mathchoice{\mathbin{\vbox{\hbox{\scalebox{0.635}{\displaystyle\bullet}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.635}{\textstyle\bullet}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.635}{\scriptstyle\bullet}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.635}{\scriptscriptstyle\bullet}}}}}}}{\bm{n}}|\right)}. Integration by parts of the velocity-pressure coupling terms along with central numerical flux functions yields

[TABLE]

The pressure gradient term and velocity divergence term are implemented in the so-called weak formulation of DG methods, cf. Ref. 8, and it holds $g_{h}{\left(p_{h},{\bm{u}}_{h}\right)}=d_{h}{\left({\bm{u}}_{h},p_{h}\right)}$ in case of exact numerical quadrature, cf. also the paragraph on numerical quadrature below. The consistent stabilisation term $j_{h}=j_{\mathrm{div},h}+j_{\mathrm{conti},h}$ is composed of a divergence penalty term $j_{\mathrm{div},h}$ and a continuity penalty term $j_{\mathrm{conti},h}$, which are given as 19

[TABLE]

where the penalty parameters of the divergence and continuity penalty terms are defined as $\tau_{\mathrm{D}}=\zeta\overline{\|{\bm{u}}^{\mathrm{ex}}_{h}\|}\frac{h}{k+1}$ and $\tau_{\mathrm{C}}=\zeta\overline{\|{\bm{u}}^{\mathrm{ex}}_{h}\|}$, respectively. Here, $\overline{\|{\bm{u}}^{\mathrm{ex}}_{h}\|}$ is the magnitude of the velocity vector averaged over the element volume and $\left(\cdot\right)^{\mathrm{ex}}$ indicates that an extrapolation of the velocity field from previous instants of time is used in the time-discrete setting, resulting in a weak formulation $j_{h}$ that is linear in ${\bm{u}}_{h}$. This definition of the penalty parameters originates from a dimensional analysis ensuring that all terms in the weak formulation (5) have consistent physical units. It was shown by means of numerical investigation in 19 that with this definition of the penalty parameters, a default value of $\zeta=1$ ensures robustness for under-resolved flows in the sense that stability has been demonstrated for a wide range of Reynolds numbers (and also in the inviscid limit), as well as for a wide range of the spatial resolution parameters $h,k$.

Numerical quadrature

The terms in the weak formulation are integrated numerically with Gaussian quadrature and $k+1$ quadrature points per coordinate direction, except for the convective term $c_{h}$ where an over-integration strategy with $\lceil\frac{3k+1}{2}\rceil$ quadrature points is used due to the quadratic nonlinearities of the convective term. On affine element geometries this setup ensures exact integration. Exact integration allows to apply integration by parts to the integrals in $c_{h}(\cdot,\cdot),~{}d_{h}(\cdot,\cdot),~{}g_{h}(\cdot,\cdot)$ which yields $g_{h}{\left(p_{h},{\bm{u}}_{h}\right)}=d_{h}{\left({\bm{u}}_{h},p_{h}\right)}$ and $c_{h}({\bm{w}}_{h};{\bm{u}}_{h},{\bm{u}}_{h})\geq 0$ for divergence-free and normal-continuous velocities ${\bm{w}}_{h}$ (and suitable boundary conditions) and arbitrary ${\bm{u}}_{h}\in{\bm{V}}_{h}$, $p_{h}\in{Q}_{h}$.

Time integration and iterative solution of linear systems of equations

Discretisation in time is based on the backward differentiation formula (BDF) time integration method using a second order accurate formulation (BDF2). The convective term is formulated explicitly in time (using a second order accurate extrapolation scheme) which results in a CFL-type restriction of the time step size, see 19 for details on the time discretisation. For the fully discrete problem, an unsteady Stokes problem has to be solved in each time step. This coupled system of equations including the stabilisation terms is solved iteratively using the GMRES method with block preconditioning. Other solution strategies for time integration such as projection methods can be used as well to obtain computationally efficient solution algorithms and have been described in 19 in the context of the present ${\bm{L}^{2}}$-based DG discretisation. The solver is stopped once the $l_{2}$ norm of the residual has been reduced by $10^{6}$ compared to the initial residual norm (where the initial solution is an extrapolation of the solution from previous time steps), or if the absolute value of the discrete $l_{2}$ norm of the residual goes below $10^{-12}$. All solver components rely on fast matrix-free operator evaluation (exploiting the so-called sum-factorisation technique) to achieve optimal computational complexity, and we refer to 30 for a documentation of the computational efficiency of this ${\bm{L}^{2}}$-based DG solver.

2.3 Divergence-free H(div)-HDG method

For the ${\bm{H}{{\left({\mathrm{div}}\right)}}}$-based method we basically rely on the method in Ref. 18; the underlying FE spaces are given by

[TABLE]

where ${\bm{V}}_{h}$ and ${Q}_{h}$ are the discrete velocity and pressure spaces, respectively, and ${\bm{\widehat{V}}}_{h}$ is the hybrid facet space which contains discrete tangential velocities on the skeleton of ${\mathcal{T}_{h}}$. Note that, thus, the global space ${\bm{V}}_{h}$ has order $k$ approximation properties but in general, locally, the discrete velocity can even be a polynomial of order $k+1$. To obtain the contribution of boundary face integrals on Dirichlet boundary facets $F\in{\mathcal{F}_{h}^{\partial}}$ with homogeneous boundary values, we set {{\big{\{}\big{\{}}\phi{\big{\}}\big{\}}}}_{F}=\phi and {\left\llbracket\phi\right\rrbracket}_{F}=\phi for the average and jump operators introduced in equation (2).

Example 2.1** (Remark 2.0).**

In order to remove some degrees of freedom (DOFs) for the velocity and all pressure unknowns except for piecewise constants, we exploit the a priori knowledge that \operatorname{\nabla\operatorname{{\mathchoice{\mathbin{\vbox{\hbox{\scalebox{0.635}{\displaystyle\bullet}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.635}{\textstyle\bullet}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.635}{\scriptstyle\bullet}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.635}{\scriptscriptstyle\bullet}}}}}}}}{\bm{u}}_{h}=0; cf. Remark 118 and Sec. 2.2.4.231. This can be achieved with a smart choice of basis functions for ${\bm{V}}_{h}$ based on an exact sequence property (De Rham complex) on the discrete level; cf. Refs. 32, 33. $\hfill\blacktriangle$

The space-semidiscrete weak formulation of the hybridised ${\bm{H}{{\left({\mathrm{div}}\right)}}}$-DG method for (1) reads as follows:

[TABLE]

Let us now specify the various terms in (14). In order to discretise the viscous term, a hybridised variant of the symmetric interior penalty (SIP) method is used. The corresponding bilinear form, with $\sigma_{K}=6{\left(k+1\right)}^{2}$, is given by

[TABLE]

where $h_{K}{\left({\bm{s}}\right)}={\left\lvert K\right\rvert}_{d}/{\left\lvert F\right\rvert}_{d-1}$ with ${\bm{s}}\in F\subset\partial K$. Here, in contrast to {\left\llbracket\cdot\right\rrbracket}, the operator ${{\llbracket}{\llbracket}\cdot{\rrbracket}{\rrbracket}}$ represents the tangential jump between cell and facet velocity projected onto the space {\ThisStyle{\ooalign{\SavedStyle\mathbb{Q}$\cr\kern-0.22pt$\SavedStyle\mathbb{Q}$\cr\kern 0.22pt$\SavedStyle\mathbb{Q}}}_{k-1}{(F)}}; that is, using the projection onto the tangent plane {\bm{v}}_{\bm{t}}={\bm{v}}-{\left({\bm{v}}\operatorname{{\mathchoice{\mathbin{\vbox{\hbox{\scalebox{0.635}{\displaystyle\bullet}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.635}{\textstyle\bullet}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.635}{\scriptstyle\bullet}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.635}{\scriptscriptstyle\bullet}}}}}}}{\bm{n}}\right)}{\bm{n}},

[TABLE]

where $\bm{\Pi}_{\mathcal{F}_{h}}^{k-1}$ is the facet-wise ${\bm{L}^{2}}$-orthogonal projection onto {\ThisStyle{\ooalign{\SavedStyle\mathbb{Q}$\cr\kern-0.22pt$\SavedStyle\mathbb{Q}$\cr\kern 0.22pt$\SavedStyle\mathbb{Q}}}_{k-1}{(F)}}. The choice $\sigma_{K}=6{\left(k+1\right)}^{2}$ for the SIP stabilisation parameter renders the SIP penalty scaling comparable to that of the ${\bm{L}^{2}}$-based method. Due to ${\bm{H}{{\left({\mathrm{div}}\right)}}}$-conformity, the pressure-velocity coupling is simply given by

[TABLE]

Here, the property \operatorname{\nabla\operatorname{{\mathchoice{\mathbin{\vbox{\hbox{\scalebox{0.635}{\displaystyle\bullet}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.635}{\textstyle\bullet}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.635}{\scriptstyle\bullet}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.635}{\scriptscriptstyle\bullet}}}}}}}}{\bm{V}}_{h}={Q}_{h} allows to test with q_{h}=\operatorname{\nabla\operatorname{{\mathchoice{\mathbin{\vbox{\hbox{\scalebox{0.635}{\displaystyle\bullet}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.635}{\textstyle\bullet}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.635}{\scriptstyle\bullet}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.635}{\scriptscriptstyle\bullet}}}}}}}}{\bm{u}}_{h} which results in the fact that we are obtaining an exactly divergence-free discrete velocity solution since discretely divergence-free velocities are actually exactly divergence-free; i.e.

[TABLE]

**Finally, for the nonlinear convection term, we use the form **

[TABLE]

where the choice $\theta=1$ leads to an upwind-stabilised formulation and $\theta=0$ corresponds to a method without any kind of convection stabilisation. Note that in the convective form (17), instead of the hybrid jump ${{\llbracket}{\llbracket}\cdot{\rrbracket}{\rrbracket}}$, the ‘normal’ DG jump {\left\llbracket\cdot\right\rrbracket}{} between neighbouring cell velocities is used.

Numerical quadrature

In the realisation of the bilinear forms, numerical quadrature is applied so that on affine linear elements, numerical integration is performed exactly. Due to polynomials of degree $k+1$ in {\ThisStyle{\ooalign{\SavedStyle\mathbb{RT}$\cr\kern-0.22pt$\SavedStyle\mathbb{RT}$\cr\kern 0.22pt$\SavedStyle\mathbb{RT}}}_{{\left[k\right]}}{}} this requires $k+2$ Gaussian quadrature points in each direction for the integrals $a_{h}(\cdot,\cdot)$; for $b_{h}(\cdot,\cdot)$ $k+1$ Gaussian quadrature points suffice as \operatorname{\nabla\operatorname{{\mathchoice{\mathbin{\vbox{\hbox{\scalebox{0.635}{\displaystyle\bullet}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.635}{\textstyle\bullet}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.635}{\scriptstyle\bullet}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.635}{\scriptscriptstyle\bullet}}}}}}}}{\ThisStyle{\ooalign{\SavedStyle\mathbb{RT}$\cr\kern-0.22pt$\SavedStyle\mathbb{RT}$\cr\kern 0.22pt$\SavedStyle\mathbb{RT}}}_{{\left[k\right]}}{}}={\mathbb{Q}_{k}{}} while, due to the nonlinearity in $c_{h}(\cdot;\cdot,\cdot)$, $\lceil\frac{3k+3}{2}\rceil$ points are required.

Time integration and linear solvers

Concerning time integration, we use the second-order Runge–Kutta variant ARS(2,2,2) of the implicit-explicit (IMEX) method introduced in Ref. 34. Here, the pressure-velocity coupling $b_{h}{\left(\cdot\right)}$ and the viscosity term $a_{h}{\left(\cdot\right)}$ is always treated implicitly in order to maintain the exactly divergence-free property of the method, and the convection term $c_{h}{\left(\cdot\right)}$ is always treated explicitly. Therefore, only linear systems have to be solved in each time step. We use static condensation in order to eliminate element-local unknowns and solve the linear systems involving the Schur complement with a BDDC-preconditioned35** CG solver.**

2.4 Comparison of approximation spaces under the divergence constraint

For both the {\ThisStyle{\ooalign{\SavedStyle\mathbb{Q}$\cr\kern-0.22pt$\SavedStyle\mathbb{Q}$\cr\kern 0.22pt$\SavedStyle\mathbb{Q}}}_{k}{}}- and the {\ThisStyle{\ooalign{\SavedStyle\mathbb{RT}$\cr\kern-0.22pt$\SavedStyle\mathbb{RT}$\cr\kern 0.22pt$\SavedStyle\mathbb{RT}}}_{{\left[k\right]}}{}}-based discretisations, the global number of velocity unknowns on a periodic Cartesian mesh with $N$ elements per coordinate direction in $d$ dimensions is given as

[TABLE]

**We note that the global number of unknowns for both velocity spaces is the same despite the fact that the local spaces are different according to (3) and (4).  **

In order to characterise the approximation properties of the two methods, it is useful to incorporate the constraints acting on the velocity space. For the ${\bm{H}{{\left({\mathrm{div}}\right)}}}$ method, the divergence-free condition $b_{h}{\left({\bm{u}}_{h},q_{h}\right)}=0$ constrains $(k+1)^{d}$ out of the $d(k+1)^{d}$ velocity degrees of freedom per element. For the ${\bm{L}^{2}}$-based method, the limit $\tau_{\mathrm{C}},\tau_{\mathrm{D}}\to\infty$ implies continuity in the normal direction and pointwise divergence-free solutions. Both imply constraints on the velocity. In Appendix A we give a derivation and characterisation of the limit method that the solution of the ${\bm{L}^{2}}$-based method converges to as $\tau_{\mathrm{C}},\tau_{\mathrm{D}}\to\infty$. Normal-continuity ($\tau_{\mathrm{C}}\to\infty$) corresponds to $(k+1)^{d-1}$ DOFs per coordinate direction, i.e. $d(k+1)^{(d-1)}$ DOFs per element, whereas the pointwise divergence-free condition ($\tau_{\mathrm{D}}\to\infty$) corresponds to $(k+1)^{d}-1$ DOFs (cf. $\dim(Q_{h}^{\text{vol}})$ in Appendix A). Further, it can be easily seen that pointwise divergence-free velocities imply the weak divergence-free condition, rendering the weak divergence-free condition $d_{h}{\left({\bm{u}}_{h},q_{h}\right)}=0$ in (5) superfluous. Taking these results together and denoting by $N_{\mathrm{DOFs},{\bm{u}}}^{\text{L2,c}}$ and $N_{\mathrm{DOFs},{\bm{u}}}^{\text{Hdiv,c}}$ the DOFs that remain after taking the constraints into account, we obtain

[TABLE]

for the DG method with infinite stabilisation. Moreover, one obtains the relation

[TABLE]

As a consequence, it makes sense to consider the bracket $\{k,k+1\}$ as polynomial degrees for the ${\bm{L}^{2}}$-DG method to compare against {\ThisStyle{\ooalign{\SavedStyle\mathbb{RT}$\cr\kern-0.22pt$\SavedStyle\mathbb{RT}$\cr\kern 0.22pt$\SavedStyle\mathbb{RT}}}_{{\left[k\right]}}{}}. In the setting $\tau_{\mathrm{C}}<\infty$ and $\tau_{\mathrm{D}}<\infty$, the velocity approximation has some freedom left that could in principle increase the solution quality and move the ${\bm{L}^{2}}$-DG method of degree $k$ closer to {\ThisStyle{\ooalign{\SavedStyle\mathbb{RT}$\cr\kern-0.22pt$\SavedStyle\mathbb{RT}$\cr\kern 0.22pt$\SavedStyle\mathbb{RT}}}_{{\left[k\right]}}{}} (and beyond). This behaviour again motivates a closer study of the two adjacent polynomial degrees with respect to Raviart–Thomas. We argue that the added solution quality by $\tau_{\mathrm{C}}<\infty$ and $\tau_{\mathrm{D}}<\infty$ can be expected to be minor because the penalty terms are necessary to render the approximation stable, dominating over “spurious” contributions.

2.5 Mass conservation

To highlight the similarities and differences of the two discretisation approaches in terms of discrete mass conservation, it is instructive to reformulate the discrete continuity equation (10) for the ${\bm{L}^{2}}$-based approach

[TABLE]

For the ${\bm{H}{{\left({\mathrm{div}}\right)}}}$-conforming approach, the velocity is normal-continuous across interior facets so that the second term in the above equation vanishes, see equation (15). By using Raviart–Thomas function spaces the velocity is exactly divergence-free for the ${\bm{H}{{\left({\mathrm{div}}\right)}}}$-conforming approach, see equation (16). Accordingly, it can be seen as a specialisation of the ${\bm{L}^{2}}$-based approach with additional restrictions for the function spaces of velocity/pressure. In case of the ${\bm{L}^{2}}$-based approach, these restrictions are not set explicitly; instead, these restrictions are enforced weakly by the use of standard function spaces along with additional stabilisation terms, see equations (11) and (12).

2.6 Energy balance and dissipation mechanisms

**For a unified presentation below one can formally write ${\bm{\widehat{u}}}_{h}={\bm{\widehat{u}}}_{h}{\left({\bm{u}}_{h}\right)}$ for the HDG method and obtain a corresponding DG bilinear form $a_{h}{\left({\left({\bm{u}}_{h},{\bm{\widehat{u}}}_{h}\right)},{\left({\bm{u}}_{h},{\bm{\widehat{u}}}_{h}\right)}\right)}\to a_{h}{\left({\bm{u}}_{h},{\bm{u}}_{h}\right)}$ (with abuse of notation) which only depends on the usual velocity variable.  **

Assuming that ${\bm{f}}\equiv{\bm{0}}$, testing symmetrically in (5) (respectively in (14) with ${\bm{\widehat{u}}}_{h}={\bm{\widehat{u}}}_{h}{\left({\bm{u}}_{h}\right)}$) leads to the discrete energy balance

[TABLE]

**Whenever $c_{h}{\left({\bm{u}}_{h};{\bm{u}}_{h},{\bm{u}}_{h}\right)}>0$ or $j_{h}{\left({\bm{u}}_{h},{\bm{u}}_{h}\right)}>0$, the corresponding non-viscosity related mechanism can clearly be characterised as artificial or numerical dissipation since there is no counterpart in the continuous Navier–Stokes model (1).  **

It may occur for some discretisations that $c_{h}{\left({\bm{u}}_{h};{\bm{u}}_{h},{\bm{u}}_{h}\right)}$ is sign-indefinite, which may lead to an artificial increase of kinetic energy. This is often compensated by additional artificial or numerical dissipation in $a_{h}{\left({\bm{u}}_{h},{\bm{u}}_{h}\right)}$ and/or $j_{h}{\left({\bm{u}}_{h},{\bm{u}}_{h}\right)}$. To illustrate this, let us explicitly state the discrete energy balance for the ${\bm{L}^{2}}$-based approach with local Lax–Friedrichs flux under the assumption of periodic boundary conditions and vanishing body forces (see the work of Fehn et al. 19 for details)

[TABLE]

**The above equation highlights the second main difference between the ${\bm{L}^{2}}$-based and ${\bm{H}{{\left({\mathrm{div}}\right)}}}$-based approaches considered in this work. The last two terms in the first row of the above equation (belonging to the convective term $c_{h}{\left({\bm{u}}_{h};{\bm{u}}_{h},{\bm{u}}_{h}\right)}$) are sign-indefinite for the ${\bm{L}^{2}}$-based approach. Note that similar terms would also arise for other formulations of the convective term, i.e., the convective formulation with upwind flux discretisation. In case of ${\bm{H}{{\left({\mathrm{div}}\right)}}}$-conforming function spaces with exactly divergence-free velocity, these two terms add up to zero by definition and only the sign-definite stabilisation term of the convective term contributes to the discrete energy evolution. In case of the stabilised ${\bm{L}^{2}}$-based approach, the aim of the divergence and continuity penalty terms (which are positive-semidefinite, i.e., they exhibit a purely dissipative character) is to compensate sign-indefinite terms of the convective term.  **

Let us also mention that in DG discretisations, $a_{h}{\left({\bm{u}}_{h},{\bm{u}}_{h}\right)}$ consists of two types of viscous dissipation: physical and numerical dissipation. Despite its frequent use in the literature, the broken gradient norm $\left\lVert\nabla_{h}{\bm{u}}_{h}\right\rVert_{\bm{L}^{2}{}}^{2}$ is not a good measure for the physical dissipation in under-resolved situations as the corresponding numerical dissipation $a_{h}{\left({\bm{u}}_{h},{\bm{u}}_{h}\right)}-\left\lVert\nabla_{h}{\bm{u}}_{h}\right\rVert_{\bm{L}^{2}{}}^{2}$ is in general not sign-definite (even if the viscosity bilinear form is coercive). We refer to Ref. 36 for an alternative evaluation of physical and numerical dissipation for DG methods. In the remainder of this work, we circumvent this problem by only considering the sum of physical and numerical viscous dissipation.

3 Decaying homogeneous isotropic turbulence:

The 3D Taylor–Green vortex problem

In the box $\Omega={\left(0,2\pi\right)}^{3}$, equipped with periodic boundary conditions on all faces, consider the case ${\bm{f}}\equiv{\bm{0}}$ and the space-periodic initial condition 37, 38

[TABLE]

This initial condition is imposed and the resulting flow is monitored over time. The resulting Taylor–Green vortex (TGV) problem is possibly the easiest flow system for which one can observe the key physical mechanisms inherent to turbulence: transition, vortex roll-up, 3D vortex stretching and interaction, and finally, in the viscous case, molecular energy dissipation. As done frequently in the literature, simulations are run for $\mbox{{Re}}=\nu^{-1}=$1600$$ until ${T}=20$; see Refs. 3, 2, 17. All simulations are performed on structured hexahedral meshes consisting of $N^{3}$ cubes and polynomial order $k$. The accuracy of the under-resolved simulations performed here is assessed by comparing the results to accurate DNS reference data from Ref. 30 with $k=7$ and $N=128$.

3.1 Comparison of methods

**At first, we want to compare the results obtained by the ${\bm{L}^{2}}$ method (5) and the ${\bm{H}{{\left({\mathrm{div}}\right)}}}$ method (14) on the basis of the evolution of kinetic energy ${\mathcal{K}}{\left({\bm{u}}_{h}\right)}=\frac{1}{2{\left\lvert\Omega\right\rvert}}\left\lVert{\bm{u}}_{h}\right\rVert_{\bm{L}^{2}{}}^{2}$ and (negative) total kinetic energy dissipation rate $-\partial_{t}{\mathcal{K}}{\left({\bm{u}}_{h}\right)}$ for the high-order choice $k=8$; see also (21). In Fig. 1 (top) both quantities are displayed on a sequence of meshes $N\in{\left\{4,8,16\right\}}$, which correspond to strongly under-resolved, moderately resolved and essentially resolved simulations. The most important observation is that the results are, in general, not significantly different with respect to the particular quantity. While on the coarse meshes, slight differences between ${\bm{L}^{2}}$ and ${\bm{H}{{\left({\mathrm{div}}\right)}}}$ results can be observed, they are basically identical on the finest mesh. Moreover, it is comparably easy to have a good approximation for the kinetic energy, even in a strongly under-resolved situation, but in order to capture the total dissipation rate, significantly more resolution is necessary. In the remainder of this section, we will thus not discuss the evolution of kinetic energy anymore, but focus on the more demanding energy dissipation rate.  **

**Let $E{\left(\kappa\right)}$ be defined as the amount of kinetic energy concentrated in the wavenumber vector ${\bm{\kappa}}\in\mathbb{R}^{n}$ with wavenumber $\kappa\coloneqq{\left\lvert{\bm{\kappa}}\right\rvert}$, the energy spectrum. Concerning the distribution of $E{\left(\kappa\right)}$ over different wavenumbers $\kappa$, Fig 1 (bottom) shows the kinetic energy spectrum at $t=10$ (shortly after the dissipation peak) for both methods under $h$-refinement. The most important observation is that the ${\bm{L}^{2}}$- and the ${\bm{H}{{\left({\mathrm{div}}\right)}}}$-based results are very similar to each other. For the considered Reynolds number and spatial resolutions, the flow exhibits only a comparably small inertial range where the Kolmogorov $-5/3$ rule can be observed. However, the inertial range becomes larger for finer resolutions for which viscous effects are better resolved and, hence, more energy gets dissipated in the fine scales. On the coarse mesh, the resolution is not sufficient to dissipate the fine-scale features on the molecular level.  **

Next, we refine the above analysis by comparing the ${\bm{H}{{\left({\mathrm{div}}\right)}}}$-conforming method of degree $k$ with the ${\bm{L}^{2}}$-conforming method of degrees $k$ and $k+1$, so that the ${\bm{H}{{\left({\mathrm{div}}\right)}}}$-conforming method is located in between the ${\bm{L}^{2}}$-conforming method in terms of degrees of freedom, see Sec. 2. We show results of this comparison in Fig. 2for $k=4$ in the top row and for $k=8$ in the bottom row. For $k=4$, the ${\bm{L}^{2}}$-conforming method of degree $k$ appears to be of similar accuracy or slightly more accurate than the ${\bm{H}{{\left({\mathrm{div}}\right)}}}$-conforming method of degree $k$. For degree $k=8$, the ${\bm{L}^{2}}$-method and ${\bm{H}{{\left({\mathrm{div}}\right)}}}$-method of degree $k$ also show similar accuracy, where the ${\bm{H}{{\left({\mathrm{div}}\right)}}}$-method tends to be slightly more accurate for the $N=4$ case. The ${\bm{L}^{2}}$-conforming method of degree $k+1$ seems to be the most accurate one in this comparison but we want to emphasise that it is difficult to draw precise conclusions due to the wriggling behaviour of the kinetic energy dissipation rate, especially for the coarser resolutions shown in the left part of the figure. Overall, these results suggest that both methods provide a similar level of accuracy for the Taylor–Green vortex problem investigated here, if the function spaces for the individual methods are chosen in a way that they offer comparable resolution capabilities in terms of DOFs. For all following considerations, comparisons between the two methods are made on the basis of the same polynomial degree $k$, which also results in a fairly comparable number of degrees of freedom according to Sec. 2.

3.2 High-order dissipation mechanisms under $h$-refinement

In Fig. 3 we consider a fixed order ($k=8$**) under $h$-refinement and display the sum of different contributions to the dissipation. Note that the plot is stacked in the sense that the solid line on the top of the orange area represents the sum of all the contributions below. For the ${\bm{L}^{2}}$-based method (top row), there are four contributions, cf. (21), while there are only two (non-zero) contributions for the ${\bm{H}{{\left({\mathrm{div}}\right)}}}$-based method. One observes that the sign-indefinite convection parts in the ${\bm{L}^{2}}$-based methods (where the corresponding part is negative, the area is coloured purple) are compensated by the stabilisation terms. For both methods, one observes that the dominant contribution in the dissipation stems from the viscosity, which eventually is the only relevant dissipation mechanism once the flow is sufficiently resolved, e.g. $(k,N)=(8,16)$.  **

3.3 High-order vs. low-order dissipation mechanisms for fixed under-resolution

**Up to now, we only considered a fixed polynomial order of $k=8$ and did $h$-refinement procedures. This section is concerned with the question whether high-order computations are actually superior to lower-order ones, given that a comparably fair environment is provided. We already saw that whenever the resolution is sufficiently large, no differences between different methods can be observed anymore with the naked eye. Therefore, here, we restrict ourselves to the strongly under-resolved situation $(k,N)=(2,16),(4,8),(8,4)$ while we note that $N_{\mathrm{DOFs},{\bm{u}}}(2,16)>N_{\mathrm{DOFs},{\bm{u}}}(4,8)>N_{\mathrm{DOFs},{\bm{u}}}(8,4)$ , cf. (18), i.e. the low-order methods have slightly more unknowns. In Fig. 4, the evolution of the kinetic energy dissipation rate (-dt_ekin) shows the differences between high-order or low-order methods. We observe that the resolution improves when going to higher-order (although using slightly fewer unknowns). However, one can also see that the total dissipation of high-order methods is concentrated in the viscous term (visc_diss); cf. (21). Consequently, in the strongly under-resolved situation, low-order methods fundamentally rely on the numerical dissipation provided by stabilisation.  **

Summarising, we believe that high-order methods more accurate than low-order methods for freely decaying turbulence problems such as the Taylor–Green vortex. Therefore, in the following, we will restrict ourselves to $k=8$ and investigate the influence of different treatments of convection, viscosity, and divergence-conformity.

3.4 Influence of different treatments of convection, viscosity and divergence-conformity

In this section, we will investigate the influence of small changes to the considered methods. Namely, we compare different fluxes for the discrete convection term for both ${\bm{L}^{2}}$ and ${\bm{H}{{\left({\mathrm{div}}\right)}}}$ method, investigate the impact of different viscosity treatments for the ${\bm{H}{{\left({\mathrm{div}}\right)}}}$ method and vary the penalisation parameter for the ${\bm{L}^{2}}$ method. As the particular differences decrease with an increasing resolution, we will only consider the strongly and mildly under-resolved cases $(k,N)=(8,4),(8,8)$.

3.4.1 Different fluxes for the convection term

**Let us briefly explain the different choices of the discrete convective term $c_{h}{\left(\cdot;\cdot,\cdot\right)}$ which are to be considered and compared here. The basic variants (also used for all computations above) are the Lax–Friedrichs form (8) for ${\bm{L}^{2}}$ and the upwind form (17) with $\theta=1$ for ${\bm{H}{{\left({\mathrm{div}}\right)}}}$. Now, we will also use the upwind form (17) with $\theta=1$ and the central flux form (17) with $\theta=0$ for the ${\bm{L}^{2}}$ method, where a single-valued normal flux is ensured via {\bm{w}}_{h}\operatorname{{\mathchoice{\mathbin{\vbox{\hbox{\scalebox{0.635}{\displaystyle\bullet}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.635}{\textstyle\bullet}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.635}{\scriptstyle\bullet}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.635}{\scriptscriptstyle\bullet}}}}}}}{\bm{n}}\to{{\big{\{}\big{\{}}{\bm{w}}_{h}{\big{\}}\big{\}}}}\operatorname{{\mathchoice{\mathbin{\vbox{\hbox{\scalebox{0.635}{\displaystyle\bullet}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.635}{\textstyle\bullet}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.635}{\scriptstyle\bullet}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.635}{\scriptscriptstyle\bullet}}}}}}}{\bm{n}}. For the ${\bm{H}{{\left({\mathrm{div}}\right)}}}$ method, we additionally regard the performance of the central flux form (17) with $\theta=0$, and the Lax–Friedrichs (LF) form (8), where \Lambda({\bm{w}}_{h})=2{\left\lvert{\bm{w}}_{h}\operatorname{{\mathchoice{\mathbin{\vbox{\hbox{\scalebox{0.635}{\displaystyle\bullet}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.635}{\textstyle\bullet}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.635}{\scriptstyle\bullet}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.635}{\scriptscriptstyle\bullet}}}}}}}{\bm{n}}\right\rvert}. Note that for exactly divergence-free ${\bm{H}{{\left({\mathrm{div}}\right)}}}$ methods (and under exact numerical integration), the LF form simply emerges from (17) by choosing $\theta=2$ and performing integration by parts of the corresponding volume term. Put differently, the difference between LF and upwinding for the ${\bm{H}{{\left({\mathrm{div}}\right)}}}$ method is solely a factor two in the positive semi-definite facet stabilisation term. On the other hand, for the ${\bm{L}^{2}}$ method, integration by parts using a not exactly divergence-free convective velocity leads to a conceptually different discrete convection term.  **

**The corresponding results for the stacked dissipation rates for $k=4$, $N=8$ simulations can be seen in Fig. 5. Regarding the total dissipation rate (-dt_ekin), it does not seem to be possible to draw a conclusion as to whether one particular convective term is superior to the others. In fact, all variants result in a comparably accurate solution even in this strongly under-resolved situation. Moreover, regarding the LF and upwinding ${\bm{H}{{\left({\mathrm{div}}\right)}}}$ results, one can observe that the amount of convective dissipation remains approximately the same regardless of which convection stabilisation is used. In particular, the additional factor two in the LF form does not lead to more numerical dissipation. Surprisingly, the third column shows that both the ${\bm{L}^{2}}$- and the ${\bm{H}{{\left({\mathrm{div}}\right)}}}$-based method can even manage the case without convection stabilisation in form of upwinding or LF for the TGV problem. We conjecture that choosing a seemingly less dissipative convection term is simply counterbalanced by other (numerical) dissipation mechanisms. Summarising, the particular choice of the discretisation of the convection term does not have a large impact on the results for this TGV problem.  **

3.4.2 Different viscosity treatment for the H(div)-HDG method

**In addition to considering variations in the convection term in the last subsection, let us now take a look at what happens when we do not use the popular SIP method for the treatment of the viscosity effects. The motivation for investigating a different method for Laplacian-like terms originates in the frequently heard criticism that the SIP parameter choice can lead to overly large penalties and, therefore, a seemingly large amount of numerical viscous dissipation. To this end, we use a more subtle mechanism to ensure discrete coercivity of the viscosity bilinear form through a so-called lifting technique. The resulting idea/concept goes back to Bassi & Rebay 39 and we will use it only in the ${\bm{H}{{\left({\mathrm{div}}\right)}}}$-HDG context; see also Ref. 40, 35.  **

**The usual penalty term in the SIP formulation is responsible for controlling the skeleton integrals involving non-quadratic forms. To obtain this control with the penalty integral, the penalty parameter has to be chosen sufficiently large, depending on a constant which, in turn, depends on polynomial trace inverse inequalities. Lifting methods tackle the problem differently by reformulating the problematic skeleton integrals as volume integrals. From this a stabilisation bilinear form can be characterised that guarantees non-negativity of the resulting bilinear form without the usual SIP jump penalty. This additional stabilisation bilinear form relies on a new variable, the discrete lifting, whose DOFs can all be eliminated locally by static condensation. Therefore, their efficiency is comparable to the corresponding method without lifting. Similarly, in the context of hybrid mixed DG methods, the additional (dual) diffusive flux variable for approximating $\nabla{\bm{u}}$ can also be eliminated locally by static condensation, thereby resulting in a comparable situation in terms of global free DOFs. In fact, many mixed formulations can be expressed as primal formulations involving lifting operations 41.  **

We define the discrete (${\bm{H}^{1}{}}$) lifting of a vector-valued function living on the boundary of each element. {thmDef} For all $K\in{\mathcal{T}_{h}}$ and ${\bm{\varphi}}\in{\bm{L}^{2}{{\left(\partial K\right)}}}$, the discrete (${\bm{H}^{1}{}}$) lifting ${\mathcal{L}{\left({\bm{\varphi}}\right)}}\in{\left.\kern-1.2pt{\bm{V}}_{h}\vphantom{|}\right|_{K}}\cap{\bm{L}_{0}^{2}{{\left(K\right)}}}$ is defined as

[TABLE]

Note that the space for the discrete lifting is chosen such that ${\mathcal{L}{\left(\cdot\right)}}$ is unique and allows to rewrite -\int_{\partial K}{{\llbracket}{\llbracket}{\bm{w}}_{h}{\rrbracket}{\rrbracket}}\operatorname{{\mathchoice{\mathbin{\vbox{\hbox{\scalebox{0.635}{\displaystyle\bullet}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.635}{\textstyle\bullet}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.635}{\scriptstyle\bullet}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.635}{\scriptscriptstyle\bullet}}}}}}}{\left(\nabla{\bm{v}}_{h}\right)}{\bm{n}}{\,{\mathrm{d}}\bm{s}}=\int_{K}\nabla{\mathcal{L}{\left({{\llbracket}{\llbracket}{\bm{w}}_{h}{\rrbracket}{\rrbracket}}\right)}}\operatorname{{\bm{\mathbin{\vbox{\hbox{\ldotp}\hbox{\ldotp}}}}}}\nabla{\bm{v}}_{h}{\,{\mathrm{d}}{\bm{x}}} in the definition of $a_{h}{\left(\cdot,\cdot\right)}$. Our ${\bm{H}{{\left({\mathrm{div}}\right)}}}$-HDG method with (${\bm{H}^{1}{}}$) lifting results from adding the form

[TABLE]

to the left-hand side of (14). One easily checks that this additional term ensures non-negativity of the bilinear form $a_{h}{\left(\cdot,\cdot\right)}$ for any $\sigma_{\ell}\geqslant 1$. Hence, we do not require the penalty parameter $\sigma_{K}$ to obey a largeness constraint (usually related to the constant in discrete trace inequalities). In fact, we simply use $\sigma_{K}=k$ and $\sigma_{\ell}=2$ and thus the (${\bm{H}^{1}{}}$**) lifting method enforces weak (tangential) continuity of the discrete ${\bm{H}{{\left({\mathrm{div}}\right)}}}$ solution in a (much) weaker sense compared to the SIP method.  **

Fig. 6 shows the dissipation rates for the ${\bm{H}{{\left({\mathrm{div}}\right)}}}$ method for $(k,N)=(8,4)$ for our basic choice SIP with upwinding, compared against the (${\bm{H}^{1}{}}$) lifting variant with upwinding, Lax–Friedrichs, and the central flux choice for the discrete convection term. Comparing the two results with upwinding, one can observe that convection dissipation plays a larger role when the (${\bm{H}^{1}{}}$) lifting method is chosen. The (${\bm{H}^{1}{}}$) lifting method is less dissipative in the viscous term, but the convection stabilisation takes over this role such that the total dissipation rate again is not different. In agreement with the last subsection, the same observation holds true for the Lax–Friedrichs convection term. When the (${\bm{H}^{1}{}}$) lifting method is used without convection stabilisation, surprisingly, it still works and indeed gives comparably accurate results. We conclude that the particular choice of the discretisation of the viscous term does not have a large impact on the results for the TGV problem.

3.4.3 Variation of penalty factors for the L2-DG method

**Next, we investigate the influence of the penalty factor $\zeta$ for the ${\bm{L}^{2}}$-based method. Fig. 7 shows results for values of $\zeta=0.1,1,10,\infty$ (from left to right and top to bottom) considering a spatial resolution of $k=8,N=8$. Note that the limit method is characterised in Appendix A. For a penalty factor to $\zeta=0.1$, the convective term exhibits a negative dissipation rate which is due to sign-indefinite terms in the energy balance, see equation (22), and instabilities might occur for smaller penalty factors $\zeta\rightarrow 0$. For $\zeta=1$, only minor undershoots can be observed where the dissipation of the convective term takes a negative value. For a penalty factor of $\zeta=10$, the numerical dissipation of the convective term is non-negative for all times and the numerical dissipation of the convective term is larger than that of the penalty terms. Finally, the results for the limit method $\zeta=\infty$ show that, by construction, the solution will be normal-continuous and exactly divergence-free. As a consequence, the only dissipation that is not related to viscosity stems from the convection stabilisation (Lax–Friedrichs). Moreover, note that for this problem, the total amount of numerical dissipation is similar for the ${\bm{L}^{2}}$-method with $\zeta=1,10$ compared to the limit method with $\zeta\to\infty$.  **

This parameter study nicely demonstrates the impact of the penalty terms on the dissipation rate of the convective term. However, let us explicitly emphasise that a potential conclusion like “the convective term stabilises the discretisation scheme by providing the required numerical dissipation for under-resolved problems” drawn from the $\zeta=10$ results is misleading. In fact, the convective term would not be able to stabilise the scheme without further action. As demonstrated in 19, the ${\bm{L}^{2}}$-based discretisation scheme without additional penalty terms does not lead to a robust discretisation in the under-resolved regime potentially leading to a blow-up of the solution. Similarly, the fact that the numerical dissipation of the continuity penalty term is small compared to the overall dissipation rate for the example considered here does not imply that the continuity penalty term is not required to obtain a robust discretisation scheme. As shown in 19, the continuity penalty term is indeed an essential ingredient for the robustness of the ${\bm{L}^{2}}$-based method. Finally, we mention that the penalty factor is not intended to be a parameter that can be used to adjust the discretisation scheme. Instead, our goal is a parameter-free turbulent flow solver and designing the penalty parameter in a way to obtain a robust discretisation scheme for a constant penalty factor that is chosen once and for all (the default value is $\zeta=1$ unless specified otherwise).

4 Turbulent channel flow

The 3D turbulent channel flow problem is a frequently used benchmark problem for assessing the ability of flow solvers to deal with wall-bounded turbulence, see Refs. 42, 43, 44. As the domain for all channel flows we consider the rectangular cuboid $\Omega={\left(0,L_{x}\right)}\times{\left(0,L_{y}\right)}\times{\left(0,L_{z}\right)}$ with $L_{x}=2\pi\delta_{c}$, $L_{y}=2\delta_{c}$, $L_{z}=\pi\delta_{c}$ and channel half-width $\delta_{c}=1$. In $x$-direction (streamwise) and $z$-direction (spanwise) periodic boundary conditions are prescribed whereas for $y\in{\left\{0,L_{y}\right\}}$ the no-slip condition ${\bm{u}}={\bm{0}}$ is imposed. Due to sharp velocity gradients, it is common practice to stretch the mesh in $y$-direction (wall normal direction). We choose the stretching function $\Phi\colon{\left[0,1\right]}\to{\left[0,L_{y}\right]}$,

[TABLE]

**with a constant $C=1.8$ for all our simulations. Therefore, given a mesh with $N$ elements in each direction, the resulting meshes for the channel flow problem consist of $N^{3}$ hexahedra which will be equidistant in $x$- and $z$-directions and stretched in $y$-direction. Usually, the whole motion is driven solely by a constant pressure gradient source term ${\bm{f}}={\left(f_{p},0,0\right)}^{\dagger}$ acting in the streamwise direction; cf., for example, Ref. 45 (Sec. 13.4).  **

**In order to distinguish different turbulent channel flow situations, the friction Reynolds number $\mbox{{Re}}_{\tau}$ is considered most frequently. It is defined as $\mbox{{Re}}_{\tau}=u_{\tau}\delta_{c}/\nu$ where $u_{\tau}$ denotes the so-called wall friction velocity which, in turn, depends on the wall shear stress $\tau_{w}$ as $u_{\tau}^{2}=\tau_{w}/\rho$. Under the assumption of a statistically steady state flow and with $\rho=1$, one obtains $\tau_{w}=f_{p}\delta_{c}$; see Ref. 46 (Sec. 2.3). Then, by choosing $f_{p}=1$, one obtains $u_{\tau}=1$ and hence, friction Reynolds number and viscosity are connected by the simple relation $\nu=1/\mbox{{Re}}_{\tau}$.  **

The relevant quantities of interest for the channel flow involve averaging in time ${\langle\cdot\rangle}_{t}$ and in spatial directions of homogeneity ${\langle\cdot\rangle}_{s}$, where the abbreviation ${\langle\cdot\rangle}={\langle{\langle\cdot\rangle}_{s}\rangle}_{t}$ is used. Then, we consider the following normalised quantities: mean velocities ${\langle u_{i}\rangle}^{+}={\langle u_{i}\rangle}/u_{\tau}$, Reynolds stresses ${\langle u_{i}^{\prime}u_{j}^{\prime}\rangle}^{+}={\langle u_{i}u_{j}\rangle}/{u_{\tau}^{2}}-{\langle u_{i}\rangle}{\langle u_{j}\rangle}/{u_{\tau}^{2}}$, and normalised root mean square (rms) velocities $u_{i,{\mathrm{rms}}}^{+}={\left\lvert{\langle u_{i}^{\prime}u_{i}^{\prime}\rangle}^{+}\right\rvert}^{\nicefrac{{1}}{{2}}}$. Furthermore, so-called wall units $y^{+}=yu_{\tau}/\nu$ for $y\in{\left[0,\delta_{c}\right]}$ are used. When normalising the numerical results, $u_{\tau}$ is the friction velocity obtained as a result of the numerical simulation and not the theoretical value $u_{\tau}=1$.

4.1 Affine vs. isoparametric finite element meshes

**In Sec. 4.2 below we will compare results based on finite elements meshes using either affine or isoparametric mappings. We use this section to explain the possible impact of the choice of finite element meshes.  **

In the case of affine linear mappings, the stretching function $\Phi$ in (24) is approximated by a continuous piecewise linear approximation $\Phi_{h}^{1}$ on an equidistant mesh for $[0,1]$ with $N$ elements while in the isoparametric case the approximation is done with $\Phi_{h}^{k}$, a continuous piecewise polynomial of degree $k$.

**In Fig. 8 both approximations and the stretching function (24) for $N=2$ and $N=4$ are sketched (for $k=3$) for illustration purposes. One can observe that for $N=2$ there is no stretching in the affine case and also for $N=4$ the difference between an affine stretching and the isoparametric stretching is significant. From this simple observation, we realise that while the analytical stretching function is used for the purpose of accumulating resolution close to the boundary, the realisation of the corresponding stretching function by the finite element mesh is essential. Note that this aspect is particularly relevant for high-order finite element discretisations considered here as opposed to low-order discretisations for which the number of elements is much larger.  **

The choice of affine or isoparametric finite element meshes also influences the choice of finite element spaces. In Section 2 we introduced the finite element spaces based on the assumption of affine linear mappings. If this assumption is no longer fulfilled and the element mappings $\Phi_{K}:\hat{K}\to K$, where $\hat{K}=[0,1]^{3}$ is the reference hexahedron, are not affine linear the local finite element spaces are adapted accordingly. Let ${\hat{\mathbb{Q}}_{k}{{\left(\hat{K}\right)}}}$ be the tensor-product polynomial space as introduced (without the hat notation) in Section 2 but w.r.t. $\hat{K}=[0,1]^{3}$ and {\hat{\ThisStyle{\ooalign{\SavedStyle\mathbb{RT}$\cr\kern-0.22pt$\SavedStyle\mathbb{RT}$\cr\kern 0.22pt$\SavedStyle\mathbb{RT}}}}_{{\left[k\right]}}{{\left(\hat{K}\right)}}} accordingly. Then, we define the mapped polynomial spaces as

[TABLE]

**where the latter is the well-known Piola transformation. Let us stress that the resulting spaces contain non-polynomial functions. From here on, we want to focus on the practical impact that a decision between affine linear and isoparametric finite element meshes can have.  **

Therefore, let us consider a simple 1D configuration with the underlying question of how the approximation properties of both approaches differ for approximating a typical (averaged) turbulent velocity profile. The following simple (explicit) formulation of a corresponding law of the wall has been derived in Reichardt 47:

[TABLE]

**The accuracy of this formula shall suffice for the intended purpose of showing the different abilities to approximate functions with sharp gradients in a boundary layer as $y\to 0$.  **

In Fig. 9, a simple 1D example demonstrates the different approaches by showing the ${L^{2}}$-best approximation of ${\langle u\rangle}$ on the inverval $y^{+}\in(0,2000)$. Here, the mesh stretching (24) is approximated with an affine (left) and an isoparametric mapping (right) of order $k=6$ where globally discontinuous polynomials of element-wise order $k=6$ are used to approximate ${\langle u\rangle}^{+}$. The resulting element boundaries are indicated by vertical lines, at least when there are more than one. One can observe that the isoparametric approach roughly allows to use a once less refined mesh in this setting. In the following section, the different approaches will be compared based on the actual 3D turbulent channel flow problem. Note that the corresponding combination of polynomial order and meshes will be identical to the ones chosen in Fig. 9.

4.2 Mesh convergence study for $\mbox{{Re}}_{\tau}=$ 395

**We investigate the accuracy of the two DG methods by comparing statistical quantities of the turbulent channel flow to accurate DNS reference data from Ref. 44. Two different polynomial degrees, $k=3$ and $k=6$, are studied for three different meshes with $N=4,8,16$ elements for $k=3$ and $N=2,4,8$ elements for $k=6$. The results are shown in Fig. 10 for $k=3$ and in Fig. 11 for $k=6$, where the ${\bm{L}^{2}}$-based discretisation and the ${\bm{H}{{\left({\mathrm{div}}\right)}}}$-based discretisation are compared, considering both affine mappings (solid lines) and isoparametric mapping (dashed lines) for each discretisation.  **

**For both polynomial degrees, the numerical results converge towards the DNS reference data under mesh refinement with accurate predictions of the mean velocity profile and Reynolds stresses on the finest meshes. Comparing the $k=3$ results to those for $k=6$, the solution quality is comparable for the coarse mesh and the intermediate mesh. Note that on the coarsest mesh, the domain is discretised by only two elements in wall normal direction for each half of the channel for $k=3$, and by only one element for $k=6$. Remarkably, even this strongly under-resolved simulations deliver an at least meaningful prediction of the involved statistical quantities, e.g., the Reynolds stress tensor. On the fine mesh, the different methods analysed here show less variations in the results for $k=6$ than for $k=3$, which can be seen as an indication of an improved accuracy of the high-order simulations with $k=6$ compared to those with $k=3$ once the flow becomes resolved, i.e., the results are less sensitive to changes in the discretisation scheme for the higher order method (note also that the number of degrees of freedom is slightly lower for the $k=6$ simulations). In contrast, no clear improvement of the $k=6$ simulations over the $k=3$ simulations can be observed in the limit of under-resolution.  **

**Comparing the results for the ${\bm{L}^{2}}$-based discretisation to the ${\bm{H}{{\left({\mathrm{div}}\right)}}}$-based discretisation, no noticeable differences between the two methods can be observed in terms of accuracy. Furthermore, the results demonstrate rigorously that the spatial resolution described by $k$ and $h$ (or $N$) is the relevant parameter that drives the accuracy of the results and that shifts the results towards the DNS data, having a much higher influence on the solution quality than the particular DG discretisation method itself.  **

**Regarding the polynomial degree used for the mapping of the geometry, a rather clear trend can be observed on the coarsest meshes: The isoparametric mapping shows more accurate results than the affine mapping. The results for the affine mapping on the coarse mesh appear too inaccurate with the profile for the mean velocity being significantly higher than the DNS reference results and the Reynolds stresses showing larger oscillations compared to the results with isoparametric mapping. Given the under-resolution of the flow field on the coarsest mesh, the mean velocity profile is predicted very accurately in case of the isoparametric mapping. However, let us point out that a seemingly exact prediction of the mean velocity profile in terms of the mean velocity in the middle of the channel for these coarse meshes might be coincidence and it cannot be guaranteed that this result is robust under variations of the DG discretisation scheme, e.g., when changing the numerical fluxes or the discretisation parameters for the convective and viscous terms. In Sec. 4.3 below, we explicitly demonstrate the impact of the DG discretisation parameters on the accuracy of the results for the coarsest mesh. With increasing resolution, the differences between the affine and isoparametric mappings diminish. This is in agreement with Fig. 9 showing that the affine mapping is of course also able to resolve sharp gradients if the mesh is fine enough.  **

Due to the above observations, in the following we restrict ourselves to showing results only for the isoparametric setting. Furthermore, as $k=6$ does not show a clear advantage compared to $k=3$ for the coarsest spatial resolutions, all following investigations are performed only for $k=3$. However, we do not expect different results (qualitatively) for different polynomial orders.

4.3 Influence of the SIP parameter and convection stabilisation

**Having analysed the convergence behaviour of the considered discretisation schemes towards the reference data under mesh refinement, it remains to demonstrate the sensitivity of the results for a given spatial resolution with respect to the parameters of the discretisation scheme. Instead of tuning the parameters to improve the solution quality for a given mesh, we want to demonstrate the “guaranteed” accuracy of the discretisation schemes, defined by the least accurate results achieved for “unsuitable” or “non-optimised” parameters. For these investigations, we focus on the coarsest mesh for degree $k=3$ since the variations in the results (i.e., the discretisation error) can be expected to be largest for the lowest spatial resolution. On the one hand, we compare results for “upwind” fluxes to central fluxes for the convective term to investigate its impact for wall-bounded turbulent flows. On the other hand, we study the impact of the SIP penalty factor, which can be expected to be the second main parameter that determines the dissipation properties of the methods. Results for these method variations are shown in Fig. 12 for the ${\bm{L}^{2}}$-method and in Fig. 13 for the ${\bm{H}{{\left({\mathrm{div}}\right)}}}$-method.  **

**For the ${\bm{L}^{2}}$-method, we multiply our basic SIP penalty parameter described above with factors of $1,2,4,8,16$ and investigate both the “standard” Lax–Friedrichs (LF) flux formulation (top row in Fig. 12) and the convective formulation of the convective term with central flux as described in Sec. 3.4.1 (bottom row in Fig. 12). For the LF formulation, the SIP penalty factor has a rather small impact on the results. The values of the mean velocity profile increase for increasing penalty factors and the curves move closer to the DNS reference data. Moreover, the wall shear stress increases slightly with increasing penalty factor. For the central flux formulation, the SIP penalty factor has a larger impact. Without upwind stabilisation in the convective term, a SIP penalty factor $\geq 2$ was required to render the method stable, which can be explained by the fact that the upwind stabilisation term and the SIP penalty term are the only terms in the discretisation scheme that weakly enforce tangential continuity in the velocity field and that the “portion” of the convective term has to be taken over by the SIP term if the upwind stabilisation is switched off. For a small SIP penalty factor of $2$, the prediction of the wall shear stress is still inaccurate, but for increasing penalty factors of $4,8,16$ the accuracy achieved with the LF formulation is retained. Analogously, switching off the upwind-like stabilisation term in the LF formulation leads to qualitatively similar results as those shown in Fig. 12 for the convective formulation with central flux. Since the LF formulation produces better results with the SIP penalty factor as defined in Sec. 2.2, the LF formulation is the preferred choice here, but we mention that a purely central convective flux can be used as well even though it remains unclear whether this extends to more complex problems. Since we want to avoid parameter adjustments for the present discretisation schemes, the LF formulation with minimal SIP penalty factor is advantageous given the fact that the iteration counts of iterative solvers will deteriorate with increasing penalty factors and given the fact that parameter tuning is not a reliable or preferable tool in under-resolved turbulence simulations in our opinion.  **

For the ${\bm{H}{{\left({\mathrm{div}}\right)}}}$ method, our basic SIP penalty parameter described above is multiplied with factors $1,2,4,16$, where $8$ is omitted only in order to improve the visibility of the results. Fig. 13 shows the resulting channel profiles with different SIP parameters and with and without upwinding (convection stabilisation). The first conclusion can be drawn from comparing SIP1-upw with SIP16-upw (dashed lines), which show that once upwinding is used, the particular parameter for the SIP mechanism does not play a major role anymore. On the other hand, SIP1 without upwinding (central) obviously leads to worse (but still stable) results. Thus, we conclude also for the ${\bm{H}{{\left({\mathrm{div}}\right)}}}$ method that a certain amount of numerical dissipation is beneficial for channel flow problems. Interestingly, when the SIP parameter is increased for the central flux situation, the curves approach the SIP1-upw curves in a continuous fashion up to the point that SIP16-central is very similar to SIP16-upw. From this behaviour we conclude that for the ${\bm{H}{{\left({\mathrm{div}}\right)}}}$ method and for this problem, the particular source for numerical dissipation, be it from SIP or from upwinding, is not decisive. A possible explanation for this phenomenon can be derived from regarding how upwinding works; cf. (17): In principle, upwinding is an anisotropic mechanism acting only on facets where in- or outflow occurs. However, in turbulent channel flow simulations, the turbulent mixing properties of the instantaneous flow basically result in the fact that upwinding acts isotropically. Thus, SIP and upwinding are penalising the tangential discontinuity in a very similar way here.

5 Summary and conclusions

**Two conceptually different discontinuous Galerkin discretisation methods for the numerical simulation of incompressible turbulent flow problems have been analysed and compared in this work. Fulfilling the incompressibility constraint as well as inter-element continuity of the velocity field in the direction normal to element interfaces are of particular importance regarding the numerical robustness and accuracy of the discretisation schemes. Remarkably, the robustness of the considered methods allows to obtain meaningful results even in significantly under-resolved settings.  **

**The present work shows that these aspects can be addressed successfully by the use of either specialised function spaces in an ${\bm{H}{{\left({\mathrm{div}}\right)}}}$-context fulfilling these properties exactly, or by plain ${\bm{L}^{2}}$-conforming spaces equipped with additional stabilisation terms enforcing the above requirements in a weak sense. To assess the properties of the numerical discretisation schemes, the three-dimensional Taylor–Green vortex problem has been investigated as a representative of transitional flows and decaying turbulence as well as the turbulent channel flow problem as a representative of wall-bounded turbulent flows. The comparison of the ${\bm{L}^{2}}$-based and ${\bm{H}{{\left({\mathrm{div}}\right)}}}$-based methods, as well as studying the impact of different discretisation variants for the convective and viscous terms and the influence of certain discretisation parameters has several important implications. Overall, variations of the discretisation scheme were found to have only a small or moderate influence on the accuracy of the results. This insensitivity can be seen as an advantage of the methods analysed here in the sense that accurate results are not only achieved for a certain set of parameters but that this accuracy can be expected in general. At the same time and as shown in this work, this implies that the possibilities to optimise high-order DG discretisations by choosing one or the other flux formulation appear to be limited as compared to the influence other parameters – most importantly the spatial resolution parameters $h$ and $k$ – have on the accuracy of the results.  **

Interestingly, it is commonly believed that the discretisation of the convective term plays a crucial role in DG discretisations rendering the method stable for convection-dominated problems by the choice of suitable flux functions such as upwind fluxes. Let us emphasise that even in under-resolved settings, basically independent of the chosen convective flux, we still have been able to obtain accurate results with respect to e.g. predictions of statistical quantities such as the Reynolds stress tensor. In this context, a key message of the present work is that – for DG discretisations of the incompressible Navier–Stokes equations – measures enforcing mass conservation and energy stability in a weak sense or in a local (exact) sense are more important to achieve a robust discretisation scheme. Indeed, the present work shows that convection stabilisation (upwind or Lax–Friedrichs) is neither necessary to ensure energy-stability of the overall method nor crucial for the robustness of the simulations. Whereas the ${\bm{L}^{2}}$ method relies on divergence and normal-continuity stabilisation for energy-stability, the divergence-free ${\bm{H}{{\left({\mathrm{div}}\right)}}}$ method is designed to be energy stable without any additional stabilisation. Thus, concerning robustness for incompressible flows, the crucial ingredient turns out to be a meaningful fulfilment of the divergence-free constraint and normal-continuity of the velocity, realised weakly (${\bm{L}^{2}}$) or strongly (${\bm{H}{{\left({\mathrm{div}}\right)}}}$**) in this work. Nonetheless, the use of upwind-like fluxes for the convective term was found advantageous, especially for the wall-bounded turbulent channel flow problem, since alternative central flux formulations have shown a larger sensitivity with respect to the discretisation parameters of the viscous term.  **

We conclude that both ${\bm{L}^{2}}$-based and ${\bm{H}{{\left({\mathrm{div}}\right)}}}$-based discretisations show convincing results in highly under-resolved scenarios and are promising candidates as generic turbulent flow solvers. Preconceptions like “high-order methods are not robust” or “no-model LES is not physical” are widespread and the results shown here – albeit promising – are certainly not enough to conclusively demonstrate the opposite. Therefore, given the maturity of the results shown here, the applicability of these novel high-order DG discretisations to complex engineering applications should be considered as part of future work.

Appendix A Characterisation of a limit method

The ${\bm{L}^{2}}$-DG method is defined by the velocity and pressure spaces

[TABLE]

**and the variational formulation (5) which depends on the stabilisation parameters $\tau_{D}$ and $\tau_{C}$, cf. (11) and (12).  **

One may ask the question if the limit $\tau_{D},\tau_{C}\to\infty$ can be characterised as a meaningful discretisation. To this end, we first realise that in the limits $\tau_{D},\tau_{C}\to\infty$, the solution ${\bm{u}}_{h}$ has to fulfil $j_{{\mathrm{div}},h}{\left({\bm{u}}_{h},{\bm{v}}_{h}\right)}=0$ and $j_{\mathrm{conti},h}{\left({\bm{u}}_{h},{\bm{v}}_{h}\right)}=0$ for all ${\bm{v}}_{h}\in{\bm{V}}_{h}$, as both bilinear forms are positive semi-definite. The subspace of functions in ${\bm{V}}_{h}$ that fulfil these constraints are pointwise divergence-free and normal-continuous across element interfaces. In this sense, the limit shares the properties of the ${\bm{H}{{\left({\mathrm{div}}\right)}}}$-based method. However, the finite element space ${\bm{V}}_{h}$ and hence the pointwise divergence-free and normal-continuous subspace

[TABLE]

**are different from the ${\bm{H}{{\left({\mathrm{div}}\right)}}}$-based method (on non-simplex meshes).  **

**Let us – for ease of presentation – assume that the mesh consists of straight, affine linearly mapped hexahedral elements. In this case, if ${\bm{u}}_{h}\in{\bm{V}}_{h}^{0}$, the weak divergence-constraint $d_{h}{\left({\bm{u}}_{h},q_{h}\right)}=0$ for all $q_{h}\in{Q}_{h}$, cf. (10), can be rewritten as $d_{h}({\bm{u}}_{h},q_{h})=g_{h}(q_{h},{\bm{u}}_{h})=0$ for all $q_{h}\in{Q}_{h}$.  **

To characterise ${\bm{V}}_{h}^{0}$ through constraints by Lagrange multiplier functions we rewrite

[TABLE]

Here $Q_{h}^{\text{skel}}$ is the space of functions defined on the (interior) skeleton of the mesh that is obtained by taking the normal trace of functions in ${\bm{V}}_{h}$. We can explicitly characterise this space as the following space of facet unknowns:

[TABLE]

This space is well-known in the community of hybrid mixed and hybrid DG methods. Less established is the space $Q_{h}^{\text{vol}}$ which, however, can also be characterised explicitly as follows:

[TABLE]

**where $\xi_{i},~{}i=1,2,3$ are the spatial coordinates. We added the zero-mean constraint for uniqueness of solutions in the later given variational formulation. This means $Q_{h}^{\text{vol}}$ is ${\mathbb{Q}_{k}{}}$ except for the element-wise highest-order bubble. Moreover, ${Q}_{h}^{\text{vol}}$ is a superset of $Q_{h}$ (since the pressure is approximated by polynomials of degree $k-1$) so that the weak divergence-constraint $d_{h}{\left({\bm{u}}_{h},q_{h}\right)}=0$ holds for all $q_{h}\in{Q}_{h}$, cf. (10), and the corresponding Lagrange-multiplier $p_{h}\in Q_{h}$ become superfluous.  **

With these characterisations, one can define the limit method as

[TABLE]

Let us note that the solution ${\bm{u}}_{h}$ will be in ${\bm{V}}_{h}^{0}$ which has dimension $\dim({\bm{V}}_{h})-\dim(Q_{h}^{\text{vol}})-\dim(Q_{h}^{\text{skel}})$.

Example A.1** (Remark A.0).**

In this section we assumed affine linear element geometries. This assumption can be dropped if the (discontinuous) finite element space ${\bm{V}}_{h}$ is defined through Piola transformations from the reference element to the physical element, cf. Section 4.1. Thereby, the normal traces and the divergence of ${\bm{V}}_{h}$ can be characterised through polynomials again. $\hfill\blacktriangle$