On a linear stability issue of split form schemes for compressible flows

Vikram Singh
The Ocean in the Earth System
Max-Planck Institute for Meteorolgy
Hamburg - 20146, Germany
vikram.singh@mpimet.mpg.de
&Praveen Chandrashekar
Center for Applicable Mathematics
Tata Institute of Fundamental Research
Bangalore – 560065, India
praveen@math.tifrbng.res.in

Abstract

Split form schemes for Euler and Navier-Stokes equations are useful for computation of turbulent flows due to their better robustness. This is because they satisfy additional conservation properties of the governing equations like kinetic energy preservation leading to a reduction in aliasing errors at high orders. Recently, linear stability issues have been pointed out for these schemes for a density wave problem and we investigate this behaviour for some standard split forms. By deriving linearized equations of split form schemes, we show that most existing schemes do not satisfy a perturbation energy equation that holds at the continuous level. A simple modification to the energy flux of some existing schemes is shown to yield a scheme that is consistent with the energy perturbation equation. Numerical tests are given using a discontinuous Galerkin method to demonstrate these results.

Keywords Split forms, summation-by-parts, DG scheme, kinetic energy preservation, linear stability

1 Introduction

High order methods are essential for large eddy simulation (LES) and direct numerical simulation (DNS) of turbulent flow problems but their stability is a major issue, especially in the under-resolved cases. High order methods suffer from aliasing instabilities that can lead to solver blow-up [1, 2, 3, 4], therefore improving robustness of these methods is an active area of research. One approach for improving robustness that has become increasingly popular is carrying over stability properties of the continuous PDE to the discrete setup, using techniques such as Summation-By-Parts operators (SBP) [5]. For example, the compressible Euler equations satisfy kinetic energy preservation (KEP) [6] and entropy preservation [7] properties under some conditions. Many different methods have been developed to extend these properties to the discrete setup using so called split forms, which can increase robustness [8, 9, 10, 6, 11, 12, 13]. However, recent studies [14, 15] have shown linear stability issues with this approach. A simple solution of the Euler equations which involves a density wave propagating with constant speed and constant pressure was considered. The linearization of certain entropy conservative schemes around this density wave solution was shown to exhibit unstable eigenvalues while the standard central scheme does not have this issue¹¹1If the semi-discrete scheme is $\frac{\textrm{d}U}{\textrm{d}t}=R(U)$ , its linearization around some state $\bar{U}$ is given by $\frac{\textrm{d}U^{\prime}}{\textrm{d}t}=R^{\prime}(\bar{U})U^{\prime}$ . For stability, all the eigenvalues of $R^{\prime}(\bar{U})$ must have negative real parts.. Many split form schemes which are designed to be kinetic energy preserving also suffer from linear instability in the sense of eigenvalues.

In this study, we analyze several split form schemes for their linear stability for the density wave problem. First, we derive linearized Euler equations around the density wave solution considered as the base solution, and show that small perturbations around the base solution are not strictly bounded by the initial data but can admit some growth. Since the base solution has small values of density, the perturbations can drive it to still smaller values, possibly negative, and make the solution to be unstable in a non-linear sense, if the underlying schemes are not monotone. While there does not exist an energy norm for the perturbations which is non-increasing, there exists a reduced energy that is conserved, but it does not bound the density perturbations.

We also linearize the split form schemes around the density wave solution, and investigate their stability in terms of conserving the reduced energy. We find that even the central scheme does not satisfy reduced energy conservation property present in the linearized PDE even though it is stable in the sense of eigenvalues. We propose a modified KEP scheme that mimics the energy principles of the linearized PDE. Finally, we present numerical results using nodal discontinuous Galerkin methods and compare some split form schemes in terms of their stability and ability to conserve quantities of interest.

2 DG split form schemes

Consider the Euler equations in conservation form which can be written as

U_{t}+F_{x}+G_{y}+H_{z}=0

(1)

where $U$ is the vector of conserved variables and $F,G,H$ are the Cartesian components of the flux vectors given by

U=\begin{bmatrix}\rho\\ \rho u\\ \rho v\\ \rho w\\ E\end{bmatrix},\qquad F=\begin{bmatrix}\rho u\\ p+\rho u^{2}\\ \rho uv\\ \rho uw\\ (E+p)u\end{bmatrix},\qquad G=\begin{bmatrix}\rho v\\ \rho uv\\ p+\rho v^{2}\\ \rho wv\\ (E+p)v\end{bmatrix},\qquad H=\begin{bmatrix}\rho w\\ \rho uw\\ \rho vw\\ p+\rho w^{2}\\ (E+p)w\end{bmatrix}

(2)

The subscripts $t,x,y,z$ in (1) refer to the time derivative and spatial derivatives respectively. In (2), $\rho,p$ are the density and pressure, $u,v,w$ are the Cartesian components of the fluid velocity, $E=p/(\gamma-1)+\rho(u^{2}+v^{2}+w^{2})/2$ is the total energy per unit volume, and $\gamma>1$ is the ratio of specific heats at constant pressure and volume.

We consider the discontinuous Galerkin (DG) schemes using nodal Lagrange polynomials collocated at the Gauss-Lobatto-Legendre (GLL) nodes to discretize the compressible Euler equations. In the following, we describe the method in 1-D only, and the extension to multiple dimensions can be made using a tensor product approach. Let $N$ be the degree of polynomials used to approximate the solution; each component of the solution inside the $e$ ’th element is written as

U^{e}=\sum_{i=0}^{N}U^{e}_{i}\ell_{i}(\xi)

where $\ell_{i}(\xi)$ are the nodal Lagrange polynomials based on $(N+1)$ -GLL nodes $-1=\xi_{0}<\xi_{1}<\ldots<\xi_{N}=+1$ . Let $\bm{U}^{e}=[U_{0}^{e},\ldots,U_{N}^{e}]^{\top}$ be the nodal values in the $e$ ’th element and the flux at these GLL nodes is $\bm{F}^{e}=[F_{0}^{e},\ldots,F_{N}^{e}]^{\top}$ . Then the semi-discrete scheme in one dimension for each component of $U$ can be written as a collocation scheme given by

\frac{\Delta x_{e}}{2}\frac{\textrm{d}\bm{U}^{e}}{\textrm{d}t}+D\bm{F}^{e}+\bm{\mathcal{M}}^{-1}\bm{\mathcal{R}}^{\top}\bm{\mathcal{B}}(\bm{F}^{*}-\bm{\mathcal{R}}\bm{F}^{e})=0

(3)

where $D\bm{F}_{e}$ approximates the flux derivative at the GLL points, $\bm{\mathcal{M}}$ is the diagonal mass matrix containing the GLL weights $w_{i}$ on its diagonal, $\bm{\mathcal{B}}=\textrm{diag}(-1,+1)$ , $\bm{\mathcal{R}}$ is the boundary restriction operator given by $\bm{\mathcal{R}}\ \bm{U}^{e}=[U_{0}^{e},\ U_{N}^{e}]^{\top}$ , $\bm{F}^{*}=[F_{e-\frac{1}{2}}^{*},F_{e+\frac{1}{2}}^{*}]^{\top}$ is the numerical flux which couples the solution to the neighbouring elements $e-1$ and $e+1$ , respectively, and $\Delta x_{e}$ is the element size. The last term containing the numerical fluxes is also called the simultaneous approximation term (SAT). In the standard DG scheme, the $D\bm{F}^{e}$ term is approximated by a differentiation matrix $\mathcal{D}$ , leading to

\frac{\Delta x_{e}}{2}\frac{\textrm{d}U_{i}^{e}}{\textrm{d}t}+\sum_{j=0}^{N}\mathcal{D}_{ij}F_{j}^{e}-\frac{\delta_{i0}}{w_{0}}(F_{e-\frac{1}{2}}^{*}-F_{0}^{e})+\frac{\delta_{iN}}{w_{N}}(F_{e+\frac{1}{2}}^{*}-F_{N}^{e})=0,\qquad 0\leq i\leq N

(4)

where $\mathcal{D}_{ij}$ are elements of the differentiation matrix which is explained below and $\delta_{ij}$ is the Kronecker symbol. For non-linear equations, we cannot prove any type of energy/entropy property because we cannot apply SBP on the $D\bm{F}^{e}$ terms. Following [13], we can modify the scheme (4) into the so called flux differencing form given by

\frac{\Delta x_{e}}{2}\frac{\textrm{d}U_{i}^{e}}{\textrm{d}t}+2\sum_{j=0}^{N}\mathcal{D}_{ij}F^{\#}_{ij}-\frac{\delta_{i0}}{w_{0}}(F_{e-\frac{1}{2}}^{*}-F_{0}^{e})+\frac{\delta_{iN}}{w_{N}}(F_{e+\frac{1}{2}}^{*}-F_{N}^{e})=0,\qquad 0\leq i\leq N

(5)

where $F^{\#}_{ij}=F^{\#}(U_{i}^{e},U_{j}^{e})$ and $F^{\#}$ is a symmetric flux function. The proper choice of the volumetric flux $F^{\#}$ and the surface flux $F^{*}$ can lead to schemes which satisfy certain energy/entropy preservation property.

The scheme (5) can also be written as a split form derivative scheme [13], where the derivatives $D\bm{F}$ of the non-linear flux in (3) are computed in a special way to achieve SBP property. Define the spatial derivatives of linear ( $a$ ), quadratic ( $ab$ ) and cubic ( $abc$ ) terms by the following formulae,

$\displaystyle D_{1}a$	$\displaystyle=$	$\displaystyle a_{x}$
$\displaystyle D_{2}(a\cdot b)$	$\displaystyle=$	$\displaystyle\frac{1}{2}[(ab)_{x}+ab_{x}+ba_{x}]$
$\displaystyle D_{3}(a\cdot b\cdot c)$	$\displaystyle=$	$\displaystyle\frac{1}{4}[(abc)_{x}+a(bc)_{x}+b(ac)_{x}+c(ab)_{x}+bca_{x}+acb_{x}+abc_{x}]$

The derivatives on the right hand side are approximated at the GLL nodes using a differentiation matrix $\bm{\mathcal{D}}$ ,

a_{x}=\bm{\mathcal{D}}\bm{a},\qquad(ab)_{x}=\bm{\mathcal{D}}\bm{a}\bm{b},\qquad(abc)_{x}=\bm{\mathcal{D}}\bm{a}\bm{b}\bm{c}

where

\bm{\mathcal{D}}_{ij}=\ell_{j}^{\prime}(\xi_{i}),\qquad 0\leq i,j\leq N

In the above equations, the product of vectors, e.g., $\bm{a}\bm{b}$ , must be interpreted componentwise.

We now list some standard split form schemes for the compressible Euler equations from the literature along with their flux form $F^{\#}$ followed by our proposed split form. The flux $F^{\#}$ depends on two states $U_{l},U_{r}$ , and we use curly braces to denote the arithmetic average

\{\phi\}=\{\phi\}_{lr}=\frac{1}{2}(\phi_{l}+\phi_{r})

The subscript $lr$ may be suppressed in some of the following formulae for clarity.

Central scheme

DF=\begin{bmatrix}D_{1}(\rho u)\\ D_{1}(p+\rho uu)\\ D_{1}(\rho uv)\\ D_{1}(\rho uw)\\ D_{1}(pu+Eu)\end{bmatrix},\qquad F^{\#}=\begin{bmatrix}\{\rho u\}\\ \{p+\rho uu\}\\ \{\rho uv\}\\ \{\rho uw\}\\ \{pu+Eu\}\end{bmatrix}

This is the canonical approach where the flux vector is differentiated and the numerical flux is just the arithmetic average of the flux at the cell boundaries. With this flux, scheme (4) and (5) are identical; but it does not have SBP property for non-linear equations, and it does not satisfy KEP or any other energy boundedness property.

Kennedy and Gruber (KG) [8]

DF=\begin{bmatrix}D_{2}(\rho\cdot u)\\ D_{1}p+D_{3}(\rho\cdot u\cdot u)\\ D_{3}(\rho\cdot u\cdot v)\\ D_{3}(\rho\cdot u\cdot w)\\ D_{2}(p\cdot u)+D_{3}(\rho\cdot e\cdot u)\end{bmatrix},\qquad F^{\#}=\begin{bmatrix}\{\rho\}\{u\}\\ \{p\}+\{\rho\}\{u\}\{u\}\\ \{\rho\}\{u\}\{v\}\\ \{\rho\}\{u\}\{w\}\\ \{p\}\{u\}+\{\rho\}\{e\}\{u\}\end{bmatrix},\qquad e=\frac{E}{\rho}

Ducros [11]

DF=\begin{bmatrix}D_{2}(\rho\cdot u)\\ D_{1}p+D_{2}(\rho u\cdot u)\\ D_{2}(\rho u\cdot v)\\ D_{2}(\rho u\cdot w)\\ D_{2}(p\cdot u)+D_{2}(E\cdot u)\end{bmatrix},\qquad F^{\#}=\begin{bmatrix}\{\rho\}\{u\}\\ \{p\}+\{\rho u\}\{u\}\\ \{\rho u\}\{v\}\\ \{\rho u\}\{w\}\\ \{p\}\{u\}+\{E\}\{u\}\end{bmatrix}

Shima et al. (KEEP-PE ) [16]

DF=\begin{bmatrix}D_{2}(\rho\cdot u)\\ D_{1}p+D_{3}(\rho\cdot u\cdot u)\\ D_{3}(\rho\cdot u\cdot v)\\ D_{3}(\rho\cdot u\cdot w)\\ \frac{1}{\gamma-1}D_{2}(p\cdot u)+T\end{bmatrix},\qquad F^{\#}=\begin{bmatrix}\{\rho\}\{u\}\\ \{p\}+\{\rho\}\{u\}\{u\}\\ \{\rho\}\{u\}\{v\}\\ \{\rho\}\{u\}\{w\}\\ \frac{1}{\gamma-1}\{p\}\{u\}+\frac{1}{2}\{\rho\}(u_{l}u_{r}+v_{l}v_{r}+w_{l}w_{r})\{u\}+\frac{1}{2}(p_{l}u_{r}+p_{r}u_{l})\end{bmatrix}

where $T$ is the split form for the kinetic energy and the $pu$ terms.

Modified KEP (mKEP)

Finally, our proposed modified KEP scheme is given by

DF=\begin{bmatrix}D_{2}(\rho\cdot u)\\ D_{1}p+D_{3}(\rho\cdot u\cdot u)\\ D_{3}(\rho\cdot u\cdot v)\\ D_{3}(\rho\cdot u\cdot w)\\ \frac{\gamma}{\gamma-1}D_{2}(p\cdot u)+D_{3}(\rho\cdot k\cdot u)\end{bmatrix},\qquad F^{\#}=\begin{bmatrix}\{\rho\}\{u\}\\ \{p\}+\{\rho\}\{u\}\{u\}\\ \{\rho\}\{u\}\{v\}\\ \{\rho\}\{u\}\{w\}\\ \frac{\gamma}{\gamma-1}\{p\}\{u\}+\{\rho\}\{k\}\{u\}\end{bmatrix},\qquad k=\frac{1}{2}(u^{2}+v^{2}+w^{2})

This is similar to the KG flux [8] and Jameson flux [6] in all the components except the energy equation, where we use only primitive variables and treat the flux in quadratic and cubic form.

The KG, KEEP-PE and mKEP schemes are kinetic energy preserving, but the central and Ducros are not.

2.1 SBP property

The differentiation matrix $\mathcal{D}$ satisfies the summation-by-parts (SBP) property given by

\bm{\mathcal{M}}\bm{\mathcal{D}}+\bm{\mathcal{D}}^{\top}\bm{\mathcal{M}}=\textrm{diag}(-1,0,\ldots,0,+1),\qquad\bm{\mathcal{M}}=\textrm{diag}(w_{0},w_{1},\ldots,w_{N})

which can be written in component form as

w_{i}\mathcal{D}_{ij}+w_{j}\mathcal{D}_{ji}=-\delta_{i0}\delta_{j0}+\delta_{iN}\delta_{jN},\qquad 0\leq i,j\leq N

(6)

The SBP property implies the following result

\bm{f}^{\top}\bm{\mathcal{M}}\bm{\mathcal{D}}\bm{g}+\bm{f}^{\top}\bm{\mathcal{D}}^{\top}\bm{\mathcal{M}}\bm{g}=f_{N}g_{N}-f_{0}g_{0}

which is the discrete analogue of the integration by parts formula

\int_{-1}^{+1}(fg_{\xi}+gf_{\xi})\textrm{d}\xi=f(+1)g(+1)-f(-1)g(-1)

Let us also note the following SBP result for later use. If $f_{ij}$ is some symmetric quantity, i.e., $f_{ij}=f_{ji}$ , then it is easy to show that

\sum_{i=0}^{N}\sum_{j=0}^{N}w_{i}\mathcal{D}_{ij}f_{ij}=-f_{00}+f_{NN}

(7)

3 One dimensional density wave

Consider a simple solution to the 1-D Euler equations given by [14]

u(x,t)=V=\textrm{constant},\qquad p(x,t)=P=\textrm{constant}

and where the density $\rho(x,t)=r(x,t)$ solves the linear advection equation

r_{t}+Vr_{x}=0

with periodic boundary conditions. We first check if the split form schemes are able to maintain the constancy of velocity and pressure when applied to solve this problem, by writing the schemes in split derivative form.

3.1 Behaviour of the KG scheme

Consider the initial condition corresponding to the density wave solution. The split flux derivatives at the initial time are given by

D_{2}(\rho\cdot u)=V\bm{\mathcal{D}}\bm{\rho},\qquad D_{1}p=0,\qquad D_{2}(\rho\cdot u\cdot u)=V^{2}\bm{\mathcal{D}}\bm{\rho}

D_{2}(p\cdot u)=0,\qquad D_{3}(\rho\cdot e\cdot u)=\frac{V^{3}}{2}\bm{\mathcal{D}}\bm{\rho}+\frac{PV}{2}\bm{\rho}\bm{\mathcal{D}}(1/\bm{\rho})+\frac{PV}{2}(1/\bm{\rho})\ \bm{\mathcal{D}}\bm{\rho}

so that

$\displaystyle\frac{\textrm{d}\bm{\rho}}{\textrm{d}t}\|_{t=0}$	$\displaystyle=$	$\displaystyle-V\bm{\mathcal{D}}\bm{\rho}+\textrm{SAT}_{\rho}$
$\displaystyle\frac{\textrm{d}\bm{m}}{\textrm{d}t}\|_{t=0}$	$\displaystyle=$	$\displaystyle-V^{2}\bm{\mathcal{D}}\bm{\rho}+\textrm{SAT}_{\rho u}$
$\displaystyle\frac{\textrm{d}\bm{E}}{\textrm{d}t}\|_{t=0}$	$\displaystyle=$	$\displaystyle-\frac{V^{3}}{2}\bm{\mathcal{D}}\bm{\rho}-\frac{PV}{2}\left[\bm{\rho}\bm{\mathcal{D}}(1/\bm{\rho})+(1/\bm{\rho})\bm{\mathcal{D}}\bm{\rho}\right]+\textrm{SAT}_{E},$

where $\textrm{SAT}_{\rho},\textrm{SAT}_{\rho u},\textrm{SAT}_{E}$ denote the SAT terms for the discretization of $\rho,\rho u,E$ respectively. At the continuous level $\rho(1/\rho)_{x}+(1/\rho)\rho_{x}=0$ , but at the discrete level $\bm{\rho}\bm{\mathcal{D}}(1/\bm{\rho})+(1/\bm{\rho})\bm{\mathcal{D}}\bm{\rho}\neq 0$ . The velocity and pressure equation can be obtained from

\frac{\textrm{d}\bm{u}}{\textrm{d}t}|_{t=0}=\frac{1}{\bm{\rho}}\left[\frac{\textrm{d}\bm{m}}{\textrm{d}t}-V\frac{\textrm{d}\bm{\rho}}{\textrm{d}t}\right]_{t=0}=0,\qquad\frac{\textrm{d}\bm{p}}{\textrm{d}t}|_{t=0}=(\gamma-1)\left[\frac{\textrm{d}\bm{E}}{\textrm{d}t}-V\frac{\textrm{d}\bm{m}}{\textrm{d}t}+\frac{1}{2}V^{2}\frac{\textrm{d}\bm{\rho}}{\textrm{d}t}\right]_{t=0}\neq 0

The scheme cannot maintain constancy of pressure, and consequently also the velocity, and it is unstable for this problem in numerical computations.

3.2 Behaviour of the modified KEP scheme

At the initial condition for the density wave solution, we have similar relations as in previous sub-section and

D_{3}(\rho\cdot k\cdot u)=\frac{V^{3}}{2}\bm{\mathcal{D}}\bm{\rho}

so that

\frac{\textrm{d}\bm{\rho}}{\textrm{d}t}|_{t=0}=-V\bm{\mathcal{D}}\bm{\rho}+\textrm{SAT}_{\rho},\qquad\frac{\textrm{d}\bm{m}}{\textrm{d}t}|_{t=0}=-V^{2}\bm{\mathcal{D}}\bm{\rho}+V\ \textrm{SAT}_{\rho},\qquad\frac{\textrm{d}\bm{E}}{\textrm{d}t}|_{t=0}=-\frac{V^{3}}{2}\bm{\mathcal{D}}\bm{\rho}+\frac{V^{3}}{2}\textrm{SAT}_{\rho}

which implies that

\frac{\textrm{d}\bm{u}}{\textrm{d}t}|_{t=0}=0,\qquad\frac{\textrm{d}\bm{p}}{\textrm{d}t}|_{t=0}=0

The scheme maintains constancy of velocity and pressure.

3.3 Behaviour of some other schemes

The central, Ducros and KEEP-PE schemes are also able to maintain constancy of pressure and velocity. However, other schemes like Kennedy and Gruber (KG), Jameson [6], Morinishi[9], Kuya et al. [17], etc., which make use of enthalpy or specific internal energy or specific total energy in the energy flux and thereby cause a coupling of density and pressure, do not maintain constancy of pressure and thereby the velocity. Such schemes are also found to be unstable in numerical computations for this test problem.

3.4 Linearized equations

We now consider a small perturbation around the density wave solution and write the resulting density, velocity and pressure as $r+\rho,V+u,P+p$ . Here $(\rho,u,p)$ are small perturbations which are governed by the linearized Euler equations given by

$\displaystyle\rho_{t}+V\rho_{x}+(ru)_{x}$	$\displaystyle=$	$\displaystyle 0$
$\displaystyle u_{t}+Vu_{x}+\frac{1}{r}p_{x}$	$\displaystyle=$	$\displaystyle 0$	(8)
$\displaystyle p_{t}+Vp_{x}+\gamma Pu_{x}$	$\displaystyle=$	$\displaystyle 0$

These equations admit a reduced energy conservation equation of the form

\mathcal{E}_{t}+(\mathcal{E}V+pu)_{x}=0,\qquad\mathcal{E}=\frac{p^{2}}{2\gamma P}+\frac{1}{2}ru^{2}

(9)

so that we have the conservation of the total energy in the domain

\int\mathcal{E}\textrm{d}x=\textrm{constant}

(10)

under periodic boundary condition. The energy $\mathcal{E}$ does not involve the density perturbations and hence we refer to it as a reduced energy. We are not able to derive an energy conservation equation that includes all three variables. See Appendix A for an alternate approach based on symmetrization where an energy equation is established but this energy is not strictly bounded with respect to time. The energy in the density perturbations is given by the equation

(\rho^{2}/2)_{t}+V(\rho^{2}/2)_{x}+\rho(ru)_{x}=0

so that

\frac{\textrm{d}}{\textrm{d}t}\int\frac{1}{2}\rho^{2}\textrm{d}x=-\int\rho(ru)_{x}\textrm{d}x

(11)

The convective term does not change the total density perturbations in the domain. The density perturbations however can grow with time if the right hand side is positive and we do not have strict control on the density perturbations.

3.5 Linearization of modified KEP scheme

Let $r,V,P$ be solution of the modified KEP scheme with constant velocity and pressure. Note that we are analyzing the semi-discrete scheme where time is still continuous. So $r$ satisfies the ODE

\frac{\Delta x_{e}}{2}\frac{\textrm{d}r_{i}^{e}}{\textrm{d}t}+2\sum_{j=0}^{N}\mathcal{D}_{ij}V\{r\}_{ij}-\frac{\delta_{i0}}{w_{0}}V(\{r\}_{e-\frac{1}{2}}-r_{0}^{e})+\frac{\delta_{iN}}{w_{N}}V(\{r\}_{e+\frac{1}{2}}-r_{N}^{e})=0

(12)

We consider a linearization around this solution. The linearized density equation is

	$\displaystyle\frac{\Delta x_{e}}{2}\frac{\textrm{d}\rho_{i}^{e}}{\textrm{d}t}+2\sum_{j=0}^{N}\mathcal{D}_{ij}[V\{\rho\}_{ij}+\{r\}_{ij}\{u\}_{ij}]-$	$\displaystyle\frac{\delta_{i0}}{w_{0}}[f_{e-\frac{1}{2}}^{\rho}-(ru+\rho V)_{0}^{e}]+\frac{\delta_{iN}}{w_{N}}[f_{e+\frac{1}{2}}^{\rho}-(ru+\rho V)_{N}^{e}]=0$		(13)
	$\displaystyle f^{*\rho}=$	$\displaystyle\ \{r\}\{u\}+V\{\rho\}$		(13)

To first order, the momentum time derivative is given by

\frac{\textrm{d}}{\textrm{d}t}(r_{i}^{e}+\rho_{i}^{e})(V+u_{i}^{e})\approx V\frac{\textrm{d}r_{i}^{e}}{\textrm{d}t}+V\frac{\textrm{d}\rho_{i}^{e}}{\textrm{d}t}+\frac{\textrm{d}}{\textrm{d}t}(r_{i}^{e}u_{i}^{e})

and using (12), the linearized momentum equation can be written as

	$\displaystyle\frac{\Delta x_{e}}{2}\frac{\textrm{d}}{\textrm{d}t}(r_{i}^{e}u_{i}^{e})+2\sum_{j=0}^{N}\mathcal{D}_{ij}[\{p\}_{ij}+V\{r\}_{ij}\{u\}_{ij}]-$	$\displaystyle\frac{\delta_{i0}}{w_{0}}[f_{e-\frac{1}{2}}^{m}-(p+rVu)_{0}^{e}]+\frac{\delta_{iN}}{w_{N}}[f_{e+\frac{1}{2}}^{m}-(p+rVu)_{N}^{e}]=0$		(14)
	$\displaystyle f^{*m}=$	$\displaystyle\ \{p\}+V\{r\}\{u\}$		(14)

To first order, the energy time derivative is

\frac{\textrm{d}E_{i}^{e}}{\textrm{d}t}\approx\frac{V^{2}}{2}\frac{\textrm{d}r_{i}^{e}}{\textrm{d}t}+\frac{V^{2}}{2}\frac{\textrm{d}\rho_{i}^{e}}{\textrm{d}t}+V\frac{\textrm{d}}{\textrm{d}t}(r_{i}^{e}u_{i}^{e})+\frac{1}{\gamma-1}\frac{\textrm{d}p_{i}^{e}}{\textrm{d}t}

Then using (13), (14), the linearized pressure equation is obtained as

	$\displaystyle\frac{\Delta x_{e}}{2}\frac{\textrm{d}p_{i}^{e}}{\textrm{d}t}+2\sum_{j=0}^{N}\mathcal{D}_{ij}[V\{p\}_{ij}+\gamma P\{u\}_{ij}]-$	$\displaystyle\frac{\delta_{i0}}{w_{0}}[f_{e-\frac{1}{2}}^{p}-(pV+\gamma Pu)_{0}^{e}]+\frac{\delta_{iN}}{w_{N}}[f_{e+\frac{1}{2}}^{p}-(pV+\gamma Pu)_{N}^{e}]=0$		(15)
	$\displaystyle f^{*p}=$	$\displaystyle\ V\{p\}+\gamma P\{u\}$		(15)

Then the perturbation energy $\mathcal{E}$ defined in (9) satisfies

\frac{\textrm{d}\mathcal{E}}{\textrm{d}t}=\frac{p}{\gamma P}\frac{\textrm{d}p}{\textrm{d}t}+\frac{\textrm{d}}{\textrm{d}t}\frac{1}{2}ru^{2}=\frac{p}{\gamma P}\frac{\textrm{d}p}{\textrm{d}t}+u\frac{\textrm{d}}{\textrm{d}t}(ru)-\frac{u^{2}}{2}\frac{\textrm{d}r}{\textrm{d}t}

Using the above identity and (12), (14), (15), we can show that the energy at a solution node satisfies

\frac{\Delta x_{e}}{2}\frac{\textrm{d}\mathcal{E}_{i}^{e}}{\textrm{d}t}+2\sum_{j=0}^{N}\mathcal{D}_{ij}\left[\frac{V}{\gamma P}p_{i}\{p\}_{ij}+p_{i}\{u\}_{ij}+u_{i}\{p\}_{ij}+Vu_{i}\{r\}_{ij}\{u\}_{ij}-\frac{1}{2}Vu_{i}^{2}\{r\}_{ij}\right]-A_{i}+B_{i}=0

where $A_{i},B_{i}$ are the interface terms at left and right faces of the element. This equation can be re-written as

\frac{\Delta x_{e}}{2}\frac{\textrm{d}\mathcal{E}_{i}^{e}}{\textrm{d}t}+\sum_{j=0}^{N}\mathcal{D}_{ij}\left[\frac{V}{\gamma P}p_{i}p_{j}+p_{i}u_{j}+u_{i}p_{j}+Vu_{i}u_{j}\{r\}_{ij}\right]-A_{i}+B_{i}=0

We multiply by the quadrature weight $w_{i}$ and sum over all $i$ ; using the symmetry of the quantity inside the square brackets and the SBP property (6), we obtain

	$\displaystyle\frac{\Delta x_{e}}{2}\frac{\textrm{d}}{\textrm{d}t}\sum_{i}w_{i}\mathcal{E}_{i}^{e}-\frac{1}{2}\left[\frac{V}{\gamma P}(p_{0}^{e})^{2}+2p_{0}^{e}u_{0}^{e}+V(u_{0}^{e})^{2}r_{0}^{e}\right]$	$\displaystyle+\frac{1}{2}\left[\frac{V}{\gamma P}(p_{N}^{e})^{2}+2p_{N}^{e}u_{N}^{e}+V(u_{N}^{e})^{2}r_{N}^{e}\right]$
		$\displaystyle-\sum_{i}w_{i}A_{i}+\sum_{i}w_{i}B_{i}=0$

The boundary term can be written as

$\displaystyle\sum_{i}w_{i}A_{i}$	$\displaystyle=$	$\displaystyle\frac{p_{0}^{e}}{\gamma P}f_{e-\frac{1}{2}}^{p}+u_{0}^{e}f_{e-\frac{1}{2}}^{m}-\frac{1}{2}(u_{0}^{e})^{2}V\{r\}_{e-\frac{1}{2}}-\left[\frac{p_{0}^{e}}{\gamma P}(p_{0}V+\gamma Pu_{0})+u_{0}^{e}(p_{0}+r_{0}Vu_{0})-\frac{1}{2}(u_{0}^{e})^{2}Vr_{0}^{e}\right]$
	$\displaystyle=$	$\displaystyle\left[\frac{V}{\gamma P}p_{0}^{e}\{p\}_{e-\frac{1}{2}}+p_{0}^{e}\{u\}_{e-\frac{1}{2}}+u_{0}^{e}\{p\}_{e-\frac{1}{2}}+\frac{1}{2}Vu_{N}^{e-1}u_{0}^{e}\{r\}_{e-\frac{1}{2}}\right]$
		$\displaystyle-\left[\frac{V}{\gamma P}(p_{0}^{e})^{2}+2p_{0}^{e}u_{0}^{e}+\frac{1}{2}V(u_{0}^{e})^{2}r_{0}^{e}\right]$

with a similar equation for the other interface term, so that the energy equation becomes

	$\displaystyle\frac{\Delta x_{e}}{2}\frac{\textrm{d}}{\textrm{d}t}\sum_{i}w_{i}\mathcal{E}_{i}^{e}-$	$\displaystyle\frac{1}{2}\left[\frac{V}{\gamma P}p_{N}^{e-1}p_{0}^{e}+p_{0}^{e}u_{N}^{e-1}+p_{N}^{e-1}u_{0}^{e}+Vu_{N}^{e-1}u_{0}^{e}\{r\}_{e-\frac{1}{2}}\right]$
	$\displaystyle+$	$\displaystyle\frac{1}{2}\left[\frac{V}{\gamma P}p_{N}^{e}p_{0}^{e+1}+p_{0}^{e+1}u_{N}^{e}+p_{N}^{e}u_{0}^{e+1}+Vu_{N}^{e}u_{0}^{e+1}\{r\}_{e+\frac{1}{2}}\right]=0$

This has a divergence form; when we add the equations from all the elements, there is a telescopic collapse of the terms inside the square brackets leading to

\frac{\textrm{d}}{\textrm{d}t}\sum_{e}\sum_{i}w_{i}\mathcal{E}_{i}^{e}=0

The total perturbation energy is conserved which is consistent with the continuous case of equation (10).

From (13), we obtain an equation for the energy in the density perturbations, which is given by

\frac{\Delta x_{e}}{2}\sum_{i}w_{i}\rho_{i}^{e}\frac{\textrm{d}\rho_{i}^{e}}{\textrm{d}t}+\sum_{i}\sum_{j=0}^{N}w_{i}\mathcal{D}_{ij}[V\rho_{i}\rho_{j}+2\rho_{i}\{r\}_{ij}\{u\}_{ij}]-\rho_{0}^{e}[f_{e-\frac{1}{2}}^{*\rho}-(ru+\rho V)_{0}^{e}]+\rho_{N}^{e}[f_{e+\frac{1}{2}}^{*\rho}-(ru+\rho V)_{N}^{e}]=0

Using SBP property (7), we get

\sum_{i}\sum_{j}w_{i}\mathcal{D}_{ij}\rho_{i}^{e}\rho_{j}^{e}=-\frac{1}{2}(\rho_{0}^{e})^{2}+\frac{1}{2}(\rho_{N}^{e})^{2}

and adding the equation from all elements, it is easy to show that the convective terms involving $V$ cancel out just as in the continuous analysis.

Remark.

The KEEP-PE flux also leads to the same linearized equations as the mKEP flux. It behaves similarly to the mKEP flux when numerically computing small perturbations as shown in the results.

3.6 Linearization of central scheme

The base flow density $r$ evolves according to (12). The linearized density equation is

\frac{\Delta x_{e}}{2}\frac{\textrm{d}\rho_{i}^{e}}{\textrm{d}t}+2\sum_{j=0}^{N}\mathcal{D}_{ij}[V\{\rho\}_{ij}+\{ru\}_{ij}]-\frac{\delta_{i0}}{w_{0}}[f_{e-\frac{1}{2}}^{*\rho}-(ru+\rho V)_{0}^{e}]+\frac{\delta_{iN}}{w_{N}}[f_{e+\frac{1}{2}}^{*\rho}-(ru+\rho V)_{N}^{e}]=0

(16)

f^{*\rho}=\{ru\}+V\{\rho\}

and the momentum equation is

\frac{\Delta x_{e}}{2}\frac{\textrm{d}}{\textrm{d}t}(r_{i}^{e}u_{i}^{e})+2\sum_{j=0}^{N}\mathcal{D}_{ij}[\{p\}_{ij}+V\{ru\}_{ij}]-\frac{\delta_{i0}}{w_{0}}[f_{e-\frac{1}{2}}^{*m}-(p+rVu)_{0}^{e}]+\frac{\delta_{iN}}{w_{N}}[f_{e+\frac{1}{2}}^{*m}-(p+rVu)_{N}^{e}]=0

(17)

f^{*m}=\{p\}+V\{ru\}

The pressure equation is same as in the case of mKEP scheme and is given by (15). We then obtain the following energy equation

\frac{\Delta x_{e}}{2}\frac{\textrm{d}\mathcal{E}_{i}^{e}}{\textrm{d}t}+\frac{1}{4}V\sum_{j=0}^{N}\mathcal{D}_{ij}u_{i}(r_{j}-r_{i})(u_{j}-u_{i})+\sum_{j=0}^{N}\mathcal{D}_{ij}\left[\frac{V}{\gamma P}p_{i}p_{j}+p_{i}u_{j}+u_{i}p_{j}+Vu_{i}u_{j}\{r\}_{ij}\right]-A_{i}+B_{i}=0

The third term is same as in the previous section and we can apply SBP on it. The second term cannot be simplified using SBP property and leads to a volumetric term with non-determinate sign in the energy equation. This is more clearly shown in the finite difference context in Appendix B. Thus the central scheme is not consistent with the energy equation (10). A similar result is obtained for the case of Ducros flux.

Remark.

Appendix B provides linearized schemes in the finite difference context where it is easier to do the analysis.

4 Numerical tests

In this section, we first present numerical results for the compressible Euler equations for three different test cases. We use the open source code Trixi.jl [18]. Time-stepping is performed with the fourth-order, five-stage, low-storage Runge-Kutta method of [19] and with CFL=0.2. The interface flux $F^{*}$ is taken to be same as the volumetric symmetric flux $F^{\#}$ , i.e., $F^{*}_{e+\frac{1}{2}}=F^{\#}(U_{N}^{e},U_{0}^{e+1})$ , since our purpose here is to investigate the properties of the scheme when there is no added numerical dissipation. All the meshes used for compressible Euler computations are comprised of uniform quadrilateral or hexahedral elements. Next, we present implicit large eddy simulations (ILES) for the compressible Navier Stokes equation for two standard test cases. These simulations are performed using the solver documented in [4, 20].

4.1 2-D density wave

Consider the initial condition given by

\rho=1+0.98\sin(2\pi(x+y)),\qquad u=0.1,\qquad v=0.2,\qquad p=20.0

in the domain $[-1,+1]^{2}$ with periodic boundary conditions. The exact solution is a translation of the density profile at constant velocity. This is an exact, smooth solution for the compressible Euler equations. However, as shown in [14, 15], many popular schemes, including entropy conserving schemes blow-up for this solution due to loss of positivity of density. We run the test until time $T=100$ and monitor the time at which blow-up happens. Table 1 shows the time to blow up using mesh of $4\times 4$ cells and $N=3,4,5$ , while Table 2 shows results for $8\times 8$ mesh. Note that unstable cases were repeated with CFL=0.05 but this did not change the results. All schemes show nearly identical time to blow-up except KG which blows up much earlier, even with increased resolution. These can be attributed to non-preservation of constant velocity and pressure condition by the scheme and the presence of unstable eigenvalues for KG flux as shown in [14].

N	Central	Ducros	KG	KEEP-PE	mKEP
3	0.51	0.51	0.13	0.51	0.51
4	0.49	0.49	0.07	0.49	0.49
5	-	-	0.08	-	-

Table 1: Time to blow-up for density wave test for

4\times 4

mesh. Tests were run until time

T=100

N	Central	Ducros	KG	KEEP-PE	mKEP
3	-	-	0.09	-	-
4	-	-	0.08	-	-
5	-	-	0.08	-	-

Table 2: Time to blow-up for density wave test for

8\times 8

mesh. Tests were run until time

T=100

Next we add a perturbation to the solution, so that the initial velocity is given by

u=0.1+A[\sin(2\pi x)+\sin(2\pi y)],\qquad v=0.2+A[\cos(2\pi x)+\cos(2\pi y)]

where $A$ is the amplitude of the perturbation. Tables 3, 4 and 5 show time to blow-up for $4\ \times\ 4$ mesh with $N=5$ , $8\ \times\ 8$ mesh with $N=3$ , and $8\ \times\ 8$ mesh with $N=4$ , respectively. Stronger perturbations with $A=10^{-3}$ result in small time to blow-up and maybe indicative of nonlinear instabilities. With $A=10^{-4}$ , we can see simulations run for much longer times until blow-up. For this amplitude, mKEP and KEEP-PE always run as long as or longer than Central and Ducros schemes. This is consistent with our linear stability analysis as mKEP and KEEP-PE satisfy the reduced energy conservation equation and are also KEP.

A	Central	Ducros	KG	KEEP-PE	mKEP
1E-3	0.80	7.48	0.08	5.79	5.79
1E-4	40.22	40.22	0.08	49.17	49.17
1E-5	-	-	0.08	-	-

Table 3: Time to blow-up for density wave test with perturbation for

4\times 4

mesh, N=5. Tests were run until time

T=100

A	Central	Ducros	KG	KEEP-PE	mKEP
1E-3	3.90	3.90	0.07	2.47	2.47
1E-4	20.83	20.83	0.10	20.84	20.84
1E-5	-	-	0.10	-	-

Table 4: Time to blow-up for density wave test with perturbation for

8\times 8

mesh, N=3. Tests were run until time

T=100

A	Central	Ducros	KG	KEEP-PE	mKEP
1E-3	2.50	9.13	0.08	9.88	9.88
1E-4	77.47	84.56	0.09	-	-
1E-5	-	-	0.09	-	-

Table 5: Time to blow-up for density wave test with perturbation for

8\times 8

mesh, N=4. Tests were run until time

T=100

4.2 2-D isentropic vortex

The isentropic vortex is another exact solution of the compressible Euler equations, and is a popular test case for evaluating high order methods. Various initial conditions are available in literature [21] for this problem. We use the initial condition given by

\rho=\left[1-\frac{\beta^{2}(\gamma-1)}{8\gamma\pi^{2}}\exp(1-r^{2})\right]^{\frac{1}{\gamma-1}},\qquad u=M\cos\alpha-\frac{\beta(y-y_{c})}{2\pi}\exp\left(\frac{1-r^{2}}{2}\right)

v=M\sin\alpha+\frac{\beta(x-x_{c})}{2\pi}\exp\left(\frac{1-r^{2}}{2}\right),\qquad r^{2}=(x-x_{c})^{2}+(y-y_{c})^{2}

and the pressure is given by $p=\rho^{\gamma}$ . We choose the parameters $\beta=5$ , $M=0.5$ , $\alpha=45^{o}$ , $(x_{c},y_{c})=(0,0)$ and the domain is taken to be $[-10,10]\times[-10,10]$ . We run the computations upto a time of $t=2\sqrt{2}\cdot 20/M=113.13708499$ when the vortex has crossed the domain two times in the diagonal direction. To measure the conservation of total kinetic energy and entropy, we plot the relative change in these quantities defined as

ke(t)=\frac{KE(t)-KE(0)}{KE(0)},\qquad en(t)=\frac{EN(t)-EN(0)}{|EN(0)|},

where $KE(t),EN(t)$ refers to total kinetic energy and entropy in the spatial domain at time $t$ . The entropy refers to the mathematical entropy which is a convex function that should be non-increasing with time. Note that we divide by $|EN(0)|$ since the total entropy can be negative; the non-increasing property of entropy requires that $en(t)$ be a non-increasing function of time.

Refer to caption — Figure 1: Comparison of kinetic energy and entropy variation with time for isentropic vortex on $32\times 32$ mesh with $N=3$ .

The simulations were performed using a $32\ \times\ 32$ uniform quadrilateral mesh with polynomial orders of $N=3,4$ . Figure 1 and 2 show the evolution of normalized kinetic energy and entropy for the various split forms. The central scheme was not stable for these simulations and is not shown in the figures. We observe that KG and mKEP behave very similarly; they are both kinetic energy dissipative and mildly entropy unstable for all simulation cases. In contrast, both Ducros and KEEP-PE schemes are nearly conservative for both kinetic energy and entropy. We observe slight kinetic energy growth for Ducros and slight entropy growth for KEEP-PE . For $N=4$ , the conservation violations are even smaller.

4.3 3-D inviscid Taylor-Green vortex

The density wave and isentropic vortex test cases were 2D cases with smooth solutions, so that turbulence effects on flow stability do not exist in those problems. We use the 3D Taylor-Green vortex case to test the numerical properties of the scheme in the presence of turbulence. Here we use the inviscid version which is useful in studying the stability of the discrete formulation [22, 16] for the Euler equations. The domain is a cube of length $2\pi$ and the initial conditions are given by

\rho=\rho_{0},\qquad u=+U_{0}\sin(x)\cos(y)\cos(z),\qquad v=-U_{0}\cos(x)\sin(y)\cos(z),\qquad w=0

p=P_{0}+\frac{\rho_{0}U_{0}^{2}}{16}(\cos(2x)\cos(2z)+2\cos(2x)+2\cos(2y)+\cos(2y)\cos(2z)),\qquad P_{0}=\frac{\rho_{0}}{M^{2}\gamma}

with $M$ being the Mach number and $\rho_{0}=1,U_{0}=1,M=0.1$ . We run the simulations on a $16\times 16\times 16$ uniform hexahedral mesh with $N=3$ .

Figure 3 shows the evolution of normalized kinetic energy and entropy for the different split forms. As with the isentropic vortex case, central scheme was not stable. All split forms nearly conserve kinetic energy and entropy until about $T=6$ , when the flow starts transitioning and becomes turbulent. The mKEP and KG split forms show nearly identical results and have very low kinetic energy dissipation. The Ducros scheme has slightly higher dissipation levels, while the KEEP-PE scheme has very high dissipation of kinetic energy after the onset of turbulence. Figure 4 shows a zoomed in view of kinetic energy evolution where the advantage of the mKEP and KG schemes over the other two schemes become even more apparent.

In the case of total entropy, KG and mKEP schemes are again nearly identical and are mildly entropy unstable; Ducros is also entropy unstable but less than mKEP, and KEEP-PE is the only scheme that is entropy stable. However, note again that the entropy instability is mild for mKEP and KEEP-PE schemes. This is in contrast to the kinetic energy dissipation which is substantial for the KEEP-PE split form.

4.4 Implicit LES results

High order DG methods and especially their split form versions are suitable for implicit large eddy simulation (ILES) of turbulent flows because of the high wave number dissipation properties of the method, which works as an LES model [23, 24, 25, 4]. To demonstrate the capabilities of the new split form, we present results for two benchmark test cases. In these simulations, we add Lax-Friedrich dissipation to the surface flux between the elements.

4.4.1 3-D viscous Taylor-Green vortex

We use the same initial data as given in Section 4.3 and solve the Navier-Stokes equations at a Reynolds number of 1600. We use a grid consisting of $64^{3}$ uniform hexahedral elements with $N=3$ . We compare the kinetic energy evolution and dissipation rate as a function of time with DNS data from [26] which is shown in Figure 5; we observe that the new split form performs as well as the KG flux and both of them agree well with the DNS data.

4.4.2 Channel flow

To show ILES capabilities for wall bounded flows, we solve the turbulent channel flow problem [27]. This flow has a periodic domain in the streamwise and spanwise directions with walls at the top and bottom where no-slip condition is imposed. The channel has height H and the half height $L=\frac{H}{2}$ decides the Reynolds’ number. The flow is driven by a forcing term to get a constant mass flow rate. The test case is characterized by the parameter $Re_{\tau}$ given by

Re_{\tau}=\frac{u_{\tau}L}{\nu},\qquad u_{\tau}=\sqrt{\frac{\tau_{w}}{\rho}}.

We use a domain length $2\pi$ in the streamwise direction and $\pi$ in the spanwise direction; the channel height is $H=2$ . It is discretized into $32\times 22\times 22$ mesh with parabolic stretching to refine the mesh near the walls. The computations were compared to DNS data from [28] for $Re_{\tau}=180$ . The initial density is taken as $\rho=1$ and the pressure is chosen such that the Mach number for the initial horizontal velocity $u_{0}$ is $M=0.2$ . To this uniform profile, perturbations are added to trigger transition. Figures 6 compares the mean streamwise velocity with DNS data, meanwhile Figure 7 shows the Reynolds stresses. These figures show excellent agreement of the simulations with DNS data for both first order and second order moments, showing the suitability of the new split form for ILES simulations.

5 Summary and conclusions

Kinetic energy and entropy preserving split forms are attractive since they can lead to improved robustness for the computation of turbulent flows, but many such schemes have been shown to be linearly unstable in the sense of eigenvalues when applied to a problem involving the advection of a density wave. The linearized Euler equations are not strongly stable for this problem since the perturbations can possibly grow with time. However, a reduced energy is conserved which we use to distinguish the various split form schemes. We present a new split form for the compressible Euler equations, which we call the modified kinetic energy preserving (mKEP) split form, which conserves the reduced energy. The central scheme which is stable in terms of eigenvalues does not conserve the reduced energy. For non-normal operators which arise by the linearization of the schemes, eigenvalues characterize only the long time behaviour, but can exhibit energy growth at initial times. We show that our new scheme is consistent with the reduced energy equation and behaves similarly to the Kennedy and Gruber flux in most test cases, while being stable for the density wave problem. Note that the KEEP-PE split form also conserves the reduced energy since it has the same linearized equations as the mKEP scheme. To test these schemes, we used a discontinuous Galerkin method which has SBP property.

For the density wave problem, schemes which cannot maintain the constancy of velocity and pressure are always unstable. Schemes which can maintain the constancy like central, Ducros, mKEP and KEEP-PE can also be unstable at coarse resolutions due to loss of positivity of density, even though they can be stable in the sense of eigenvalues. At sufficiently high resolutions, such schemes are found to be stable for long times. When the initial condition is perturbed, mKEP and KEEP-PE schemes behave similarly and are more stable than central and Ducros schemes. If the initial perturbation is sufficiently large, all schemes fail at large times due to loss of positivity of density, since none of them is monotone.

For the 2D isentropic vortex problem, the KEEP-PE and Ducros schemes show better kinetic energy conservation, with Ducros showing a small growth in kinetic energy, while KG and mKEP are sightly more dissipative. The KEEP-PE and Ducros are similarly entropy dissipative while the KG and mKEP produce small amounts of entropy. For all schemes, these errors decrease with increasing resolution.

For the 3D invsicid Taylor-Green problem, the KEEP-PE scheme shows very high kinetic energy dissipation after the transition to turbulence. The mKEP and KG schemes behave similarly and have small levels of kinetic energy dissipation, while the Ducros scheme has higher dissipation but much less than KEEP-PE scheme. The KEEP-PE scheme is entropy dissipative while KG, mKEP and Ducros schemes produce small amounts of entropy.

These observations suggest that the mKEP scheme is a good alternative to existing split form schemes when used in a DG method, since it is more stable for the density wave problem, satisfies the reduced energy conservation and has good kinetic energy properties similar to the KG flux. It also performs well for the two viscous turbulent test cases shown here in an ILES setting. The KEEP-PE scheme has been shown to perform well in finite difference schemes but it seems to have high levels of kinetic energy dissipation when used in DG methods for turbulent flows. These issues have to be further investigated by applying the schemes to more realistic problems.

Acknowledgements

We thank the authors of the Julia code Trixi.jl [18] for making it publicly available. Praveen Chandrashekar was supported by Department of Atomic Energy, Government of India, under project no. 12-R&D-TFR-5.01-0520.

Appendix A Symmetrization of linearized equations for density wave problem

We attempt to derive an energy equation for the linearized Euler equations (8) by following the symmetrization approach of [29]. Let $V=[\rho,u,p]$ denote small perturbations and define the change of variables

W=S^{-1}V,\qquad S^{-1}=\begin{bmatrix}\frac{c}{\sqrt{\gamma}r}&0&0\\ 0&1&0\\ -\frac{c}{r\sqrt{\gamma(\gamma-1)}}&0&\sqrt{\frac{\gamma}{\gamma-1}}\frac{1}{rc}\end{bmatrix},\qquad c=c(x,t)=\sqrt{\frac{\gamma P}{r(x,t)}}

Note that the matrix $S$ is a function of $x,t$ and is not constant. The linearized equations can be written as

W_{t}+\tilde{A}W_{x}+\tilde{B}W=0,\qquad\tilde{A}=\begin{bmatrix}V&\frac{c}{\sqrt{\gamma}}&0\\ \frac{c}{\sqrt{\gamma}}&V&\sqrt{\frac{\gamma-1}{\gamma}}c\\ 0&\sqrt{\frac{\gamma-1}{\gamma}}c&V\end{bmatrix},\qquad\tilde{B}=\begin{bmatrix}0&\frac{c}{r\sqrt{\gamma}}&0\\ \frac{c}{2r\sqrt{\gamma}}&0&\sqrt{\frac{\gamma-1}{\gamma}}\frac{c}{2r}\\ 0&-\frac{c}{r\sqrt{\gamma(\gamma-1)}}&0\end{bmatrix}\frac{\partial r}{\partial x}

If the base flow density is constant, $r=$ constant and $\tilde{B}=0$ , then we obtain a symmetric hyperbolic system with constant coefficients for which the quadratic energy $W^{\top}W$ is conserved. In the general case, $\tilde{A},\tilde{B}$ are functions of $(x,t)$ and we get the energy equation

W^{\top}W_{t}+W^{\top}\tilde{A}W_{x}+W^{\top}\tilde{B}W=0

\frac{1}{2}(W^{\top}W)_{t}+(W^{\top}\tilde{A}W)_{x}+W^{\top}\tilde{C}W=0

where

\tilde{C}=\tilde{B}-\frac{1}{2}\tilde{A}_{x}=\begin{bmatrix}0&\frac{5c}{4r\sqrt{\gamma}}&0\\ \frac{3c}{4r\sqrt{\gamma}}&0&\sqrt{\frac{\gamma-1}{\gamma}}\frac{3c}{4r}\\ 0&\frac{(\gamma-5)c}{4r\sqrt{\gamma(\gamma-1)}}&0\end{bmatrix}\frac{\partial r}{\partial x}

The solution is bounded by

\frac{\textrm{d}}{\textrm{d}t}|W|^{2}\leq 2c|W|^{2}\qquad\implies\qquad|W(t)|\leq\exp(cT)|W(0)|,\qquad 0\leq t\leq T

where $c=\|\tilde{C}\|$ . We are not able to strictly bound the perturbation energy in terms of the initial energy, and the equations do allow for some growth in the solution.

Appendix B Finite difference schemes

In this section, we perform linear stability analysis of some split form fluxes in the finite difference case for the density wave problem. Consider a finite difference semi-discrete scheme for 1-D Euler equations

\frac{\textrm{d}U_{e}}{\textrm{d}t}+\frac{F_{e+\frac{1}{2}}-F_{e-\frac{1}{2}}}{\Delta x}=0

where $F_{e+\frac{1}{2}}=F(Q_{e},Q_{e+1})$ is a numerical flux function. We write the numerical flux in terms of the primitive variables $Q=[\rho,u,p]$ which simplifies our analysis. Some examples of numerical fluxes are given in Section 2.

B.1 Density wave problem in 1-D

Consider the density wave problem in 1-D introduced in Section 3, whose exact solution is of the form $\rho(x,t)=r(x,t)=r(x-Vt,0),u=V,p=P$ . In this case the Euler flux derivative is of the form $F_{x}=[\rho_{x}V,\ \rho_{x}V^{2},\frac{1}{2}\rho_{x}V^{3}]^{\top}$ and all three equations reduce to an advection equation for density. We investigate this property for some numerical schemes.

KG flux

At time $t=0$ , the density and momentum equations reduce to an advection equation

\frac{\textrm{d}\rho_{e}}{\textrm{d}t}|_{t=0}+V\frac{\{\rho\}_{e+\frac{1}{2}}-\{\rho\}_{e-\frac{1}{2}}}{\Delta x}=0

which is just the central difference scheme. The energy flux is

F_{E}=\frac{1}{\gamma-1}\{\rho\}\{p/\rho\}+\{p\}\{u\}+\frac{1}{2}\{\rho\}\{u^{2}\}\{u\}

and the energy equation yields

\frac{\textrm{d}p_{e}}{\textrm{d}t}|_{t=0}+PV\frac{\{\rho\}_{e+\frac{1}{2}}\{1/\rho\}_{e+\frac{1}{2}}-\{\rho\}_{e-\frac{1}{2}}\{1/\rho\}_{e-\frac{1}{2}}}{\Delta x}=0

The second term on the left is not zero; a Taylor expansion around $x=x_{e}$ gives

\frac{\textrm{d}p_{e}}{\textrm{d}t}|_{t=0}+\frac{PV}{2\rho_{e}^{3}}[\rho_{e}\rho^{\prime}_{e}\rho^{\prime\prime}_{e}-(\rho^{\prime}_{e})^{3}]\Delta x^{3}+O(\Delta x^{4})=0

where primes denote spatial derivatives at $t=0$ . While the error in pressure is $O(\Delta x^{2})$ , its magnitude can be large on coarse grids, expecially in regions where the density $\rho_{e}$ is small and the pressure and/or velocity are large. When we evolve the solution in time, the scheme changes the pressure distribution in space and thus creates an artificial pressure gradient out of the density gradient, which then changes the velocity also.

Other fluxes.

A similar effect arises with the KEP flux of Jameson [6] due to the use of enthalpy average, Kuya et al. flux [17], and all other fluxes in which the pressure and density are coupled in some average. These schemes do not keep the pressure constant which then changes both the velocity and pressure profiles. When we compute the solution using such fluxes, the computations break down due to loss of positivity of solution. The Ducros flux and the modified KEP flux maintain constancy of pressure and velocity.

B.2 Linearized schemes for density wave problem

Let $Q_{e}^{0}(t)=[r_{e}(t),\ V,\ P]^{\top}$ be solution of the semi-discrete finite volume scheme with constant velocity and pressure. For a scheme to admit such a solution, we require that the momentum and energy fluxes satisfy the following conditions

F_{m}(Q_{e}^{0},Q_{e+1}^{0})=P+VF_{\rho}(Q_{e}^{0},Q_{e+1}^{0}),\qquad F_{E}(Q_{e}^{0},Q_{e+1}^{0})=\frac{\gamma PV}{\gamma-1}+\frac{1}{2}V^{2}F_{\rho}(Q_{e}^{0},Q_{e+1}^{0})

(18)

in which case, all three equations reduce to the density equation. The central, Ducros and modified KEP fluxes satisfy this condition.

Consider a perturbed solution of the form $Q_{e}^{0}+q_{e}$ where $Q_{e}^{0}$ is the density wave solution and $q_{e}$ is a small perturbation around this. Let us denote the arguments of the numerical flux by $(X,Y)\to F(X,Y)$ and let $F_{X},F_{Y}$ denote jacobians with respect to the two arguments. Then

F(Q_{e}^{0}+q_{e},Q_{e+1}^{0}+q_{e+1})=F(Q_{e}^{0},Q_{e+1}^{0})+F_{X}(Q_{e}^{0},Q_{e+1}^{0})q_{e}+F_{Y}(Q_{e}^{0},Q_{e+1}^{0})q_{e+1}+O(|q|^{2})

The linearized scheme is

	$\displaystyle\Delta x\begin{bmatrix}\frac{\textrm{d}\rho_{e}}{\textrm{d}t}\\ V\frac{\textrm{d}\rho_{e}}{\textrm{d}t}+\frac{\textrm{d}}{\textrm{d}t}(r_{e}u_{e})\\ \frac{V^{2}}{2}\frac{\textrm{d}\rho_{e}}{\textrm{d}t}+V\frac{\textrm{d}}{\textrm{d}t}(r_{e}u_{e})+\frac{1}{\gamma-1}\frac{\textrm{d}p_{e}}{\textrm{d}t}\end{bmatrix}+$	$\displaystyle F_{X}(Q_{e}^{0},Q_{e+1}^{0})q_{e}+F_{Y}(Q_{e}^{0},Q_{e+1}^{0})q_{e+1}$
	$\displaystyle-$	$\displaystyle F_{X}(Q_{e-1}^{0},Q_{e}^{0})q_{e-1}-F_{Y}(Q_{e-1}^{0},Q_{e}^{0})q_{e}=0$

Let us also note the following equations for later use

\frac{\textrm{d}}{\textrm{d}t}\left(\frac{1}{2}r_{e}u_{e}^{2}\right)=u_{e}\frac{\textrm{d}}{\textrm{d}t}(r_{e}u_{e})-\frac{u_{e}^{2}}{2}\frac{\textrm{d}r_{e}}{\textrm{d}t},\qquad\Delta x\frac{\textrm{d}r_{e}}{\textrm{d}t}=-V\Delta\{r\}_{e}

where $\Delta(\cdot)_{e}=(\cdot)_{e+\frac{1}{2}}-(\cdot)_{e-\frac{1}{2}}$ .

B.2.1 Central flux

The flux jacobian is

F_{X}(Q_{e}^{0},Q_{e+1}^{0})=\begin{bmatrix}\frac{1}{2}V&\frac{1}{2}r_{e}&0\\ \\ \frac{1}{2}V^{2}&r_{e}V&\frac{1}{2}\\ \\ \frac{1}{4}V^{3}&\frac{\gamma P}{2(\gamma-1)}+\frac{3}{4}r_{e}V^{2}&\frac{\gamma V}{2(\gamma-1)}\end{bmatrix},\quad F_{Y}(Q_{e}^{0},Q_{e+1}^{0})=\begin{bmatrix}\frac{1}{2}V&\frac{1}{2}\{r\}&0\\ \\ \frac{1}{2}V^{2}&\frac{1}{2}r_{e+1}V+\frac{1}{2}\{r\}V&\frac{1}{2}\\ \\ \frac{1}{4}V^{3}&\frac{\gamma P}{2(\gamma-1)}+\frac{1}{2}r_{e+1}V^{2}+\frac{1}{4}\{r\}V^{2}&\frac{\gamma V}{2(\gamma-1)}\end{bmatrix}

where $\{r\}=\{r\}_{e+\frac{1}{2}}=\{r\}_{e,e+1}$ . The linearized scheme is given by

$\displaystyle\Delta x\frac{\textrm{d}\rho_{e}}{\textrm{d}t}+V\Delta\{\rho\}_{e}+\Delta\{ru\}_{e}$	$\displaystyle=$	$\displaystyle 0$
$\displaystyle\Delta x\frac{\textrm{d}}{\textrm{d}t}(r_{e}u_{e})+V\Delta\{ru\}_{e}+\Delta\{p\}_{e}$	$\displaystyle=$	$\displaystyle 0$
$\displaystyle\Delta x\frac{\textrm{d}p_{e}}{\textrm{d}t}+V\Delta\{p\}_{e}+\gamma P\Delta\{u\}_{e}$	$\displaystyle=$	$\displaystyle 0$

From the above equations, we can derive the following relations

	$\displaystyle\Delta x\frac{\textrm{d}}{\textrm{d}t}\left(\frac{1}{2}r_{e}u_{e}^{2}\right)$	$\displaystyle=$	$\displaystyle-Vu_{e}\Delta\{ru\}_{e}-u_{e}\Delta\{p\}_{e}+\frac{1}{2}u_{e}^{2}V\Delta\{r\}_{e}$
	$\displaystyle\Delta x\frac{\textrm{d}}{\textrm{d}t}\left(\frac{1}{2\gamma P}p_{e}^{2}\right)$	$\displaystyle=$	$\displaystyle-\frac{V}{\gamma P}p_{e}\Delta\{p\}_{e}-p_{e}\Delta\{u\}_{e}$

We add the above two equations and sum over all cells. Note that

	$\displaystyle\sum_{e}p_{e}\Delta\{p\}_{e}$	$\displaystyle=$	$\displaystyle\frac{1}{2}\sum_{e}(p_{e}^{2}-p_{e+1}^{2})=0$
	$\displaystyle\sum_{e}\left[u_{e}\Delta\{p\}_{e}+p_{e}\Delta\{u\}_{e}\right]$	$\displaystyle=$	$\displaystyle\sum_{e}(p_{e}u_{e}-p_{e+1}u_{e+1})=0$

We obtain the energy equation

\Delta x\sum_{e}\frac{\textrm{d}\mathcal{E}_{e}}{\textrm{d}t}=-V\sum_{e}\left[u_{e}\Delta\{ru\}_{e}-\frac{1}{2}u_{e}^{2}\Delta\{r\}_{e}\right]=\frac{V}{2}\sum_{e}(r_{e+1}-r_{e})(u_{e+1}-u_{e})^{2}\neq 0

We cannot conclude from this that the perturbation energy will be bounded with time.

B.2.2 Ducros flux

The flux jacobian is

F_{X}(Q_{e}^{0},Q_{e+1}^{0})=\begin{bmatrix}\frac{1}{2}V&\frac{1}{2}\{r\}&0\\ \\ \frac{1}{2}V^{2}&\frac{1}{2}r_{e}V+\frac{1}{2}\{r\}V&\frac{1}{2}\\ \\ \frac{1}{4}V^{3}&\frac{\gamma P}{2(\gamma-1)}+\frac{1}{2}r_{e}V^{2}+\frac{1}{4}\{r\}V^{2}&\frac{\gamma V}{2(\gamma-1)}\end{bmatrix},\quad F_{Y}(Q_{e}^{0},Q_{e+1}^{0})=\begin{bmatrix}\frac{1}{2}V&\frac{1}{2}\{r\}&0\\ \\ \frac{1}{2}V^{2}&\frac{1}{2}r_{e+1}V+\frac{1}{2}\{r\}V&\frac{1}{2}\\ \\ \frac{1}{4}V^{3}&\frac{\gamma P}{2(\gamma-1)}+\frac{1}{2}r_{e+1}V^{2}+\frac{1}{4}\{r\}V^{2}&\frac{\gamma V}{2(\gamma-1)}\end{bmatrix}

The linearized scheme is given by

$\displaystyle\Delta x\frac{\textrm{d}\rho_{e}}{\textrm{d}t}+V\Delta\{\rho\}_{e}+\Delta(\{r\}\{u\})_{e}$	$\displaystyle=$	$\displaystyle 0$
$\displaystyle\Delta x\frac{\textrm{d}}{\textrm{d}t}(r_{e}u_{e})+V\Delta\{ru\}_{e}+\Delta\{p\}_{e}$	$\displaystyle=$	$\displaystyle 0$
$\displaystyle\Delta x\frac{\textrm{d}p_{e}}{\textrm{d}t}+V\Delta\{p\}_{e}+\gamma P\Delta\{u\}_{e}$	$\displaystyle=$	$\displaystyle 0$

The momentum and pressure equations are same as in the case of central flux and hence we obtain the same equation for perturbation energy as in the case of central flux, and this scheme also does not conserve the perturbation energy $\mathcal{E}$ .

B.2.3 KG flux without advection

The KG flux does not satisfy the condition (18) required to preserve constant velocity and pressure. It does trivially satisfy this condition if there is no advection, i.e., $V=0$ . We will analyze the linear stability in this special case. The flux jacobian is given by

F_{X}(Q_{e}^{0},Q_{e+1}^{0})=F_{Y}(Q_{e}^{0},Q_{e+1}^{0})=\begin{bmatrix}0&\frac{1}{2}\{r\}&0\\ 0&0&\frac{1}{2}\\ 0&\frac{1}{2}\alpha P&0\end{bmatrix},\qquad\alpha_{e+\frac{1}{2}}=1+\frac{\{r\}_{e+\frac{1}{2}}\{1/r\}_{e+\frac{1}{2}}}{\gamma-1}

The linearized equations are given by

$\displaystyle\Delta x\frac{\textrm{d}\rho_{e}}{\textrm{d}t}+\Delta(\{r\}\{u\})_{e}$	$\displaystyle=$	$\displaystyle 0$
$\displaystyle\Delta x\ r_{e}\frac{\textrm{d}u_{e}}{\textrm{d}t}+\Delta\{p\}_{e}$	$\displaystyle=$	$\displaystyle 0$
$\displaystyle\Delta x\frac{\textrm{d}p_{e}}{\textrm{d}t}+(\gamma-1)P\Delta(\alpha\{u\})_{e}$	$\displaystyle=$	$\displaystyle 0$

Let

\{r\}\{1/r\}=1+\beta,\qquad\textrm{then}\qquad\alpha=\frac{\gamma}{\gamma-1}+\frac{\beta}{\gamma-1}

and the pressure equation can be written as

\Delta x\frac{\textrm{d}p_{e}}{\textrm{d}t}+\gamma P\Delta\{u\}_{e}+P\Delta(\beta\{u\})_{e}=0

The first two terms are consistent with the linearized pressure equation but the last term is an error term which goes to zero as $\Delta\to 0$ since $\beta\to 0$ . However, on any finite mesh, the total energy evolves as

\Delta x\sum_{e}\frac{\textrm{d}\mathcal{E}_{e}}{\textrm{d}t}=-\frac{1}{\gamma}\sum_{e}p_{e}\Delta(\beta\{u\})_{e}=\frac{1}{\gamma}\sum_{e}\beta_{e+\frac{1}{2}}\{u\}_{e+\frac{1}{2}}(p_{e+1}-p_{e})\neq 0

Even in this case of no advection, the KG flux is not consistent with the perturbation equation and we cannot conclude stability of this scheme since the sign of the right hand side is indeterminate.

B.2.4 Modified KEP flux

The flux jacobian is

F_{X}(Q_{e}^{0},Q_{e+1}^{0})=F_{Y}(Q_{e}^{0},Q_{e+1}^{0})=\begin{bmatrix}\frac{1}{2}V&\frac{1}{2}\{r\}&0\\ \\ \frac{1}{2}V^{2}&\{r\}V&\frac{1}{2}\\ \\ \frac{1}{4}V^{3}&\frac{\gamma P}{2(\gamma-1)}+\frac{3}{4}\{r\}V^{2}&\frac{\gamma V}{2(\gamma-1)}\end{bmatrix}

and the linearized scheme is given by

$\displaystyle\Delta x\frac{\textrm{d}\rho_{e}}{\textrm{d}t}+V\Delta\{\rho\}_{e}+\Delta(\{r\}\{u\})_{e}$	$\displaystyle=$	$\displaystyle 0$
$\displaystyle\Delta x\frac{\textrm{d}}{\textrm{d}t}(r_{e}u_{e})+V\Delta(\{r\}\{u\})_{e}+\Delta\{p\}_{e}$	$\displaystyle=$	$\displaystyle 0$
$\displaystyle\Delta x\frac{\textrm{d}p_{e}}{\textrm{d}t}+V\Delta\{p\}_{e}+\gamma P\Delta\{u\}_{e}$	$\displaystyle=$	$\displaystyle 0$

These equations lead to

	$\displaystyle\Delta x\frac{\textrm{d}}{\textrm{d}t}\left(\frac{1}{2}r_{e}u_{e}^{2}\right)$	$\displaystyle=$	$\displaystyle-Vu_{e}\Delta(\{r\}\{u\})_{e}-u_{e}\Delta\{p\}_{e}+\frac{1}{2}u_{e}^{2}V\Delta\{r\}_{e}$
	$\displaystyle\Delta x\frac{\textrm{d}}{\textrm{d}t}\left(\frac{1}{2\gamma P}p_{e}^{2}\right)$	$\displaystyle=$	$\displaystyle-\frac{V}{\gamma P}p_{e}\Delta\{p\}_{e}-p_{e}\Delta\{u\}_{e}$

We add the two equations and sum over all cells to obtain the energy equation

\Delta x\sum_{e}\frac{\textrm{d}\mathcal{E}_{e}}{\textrm{d}t}=-V\sum_{e}\left[u_{e}\Delta(\{r\}\{u\})_{e}-\frac{1}{2}u_{e}^{2}\Delta\{r\}_{e}\right]=0

Thus the modified KEP scheme preserves the perturbation energy and is stable in this sense. The density perturbations satisfy the equation

\Delta x\sum_{e}\frac{\textrm{d}}{\textrm{d}t}\frac{1}{2}\rho_{e}^{2}=-V\sum_{e}\rho_{e}\{\rho\}_{e}-\sum_{e}\rho_{e}\Delta(\{r\}\{u\})_{e}=-\sum_{e}\rho_{e}\Delta(\{r\}\{u\})_{e}

which is consistent with the continuous equation (11), i.e., the convective terms do not lead to any accumulation of density perturbations.

References

[1] R.C. Moura, G. Mengaldo, J. Peiró, and S.J. Sherwin. On the eddy-resolving capability of high-order discontinuous Galerkin approaches to implicit LES / under-resolved DNS of Euler turbulence. Journal of Computational Physics, 330:615 – 623, 2017.
[2] Sergio Pirozzoli. Generalized conservative approximations of split convective derivative operators. Journal of Computational Physics, 229(19):7180–7190, September 2010.
[3] Diego Rojas, Radouan Boukharfane, Lisandro Dalcin, David C. Del Rey Fernández, Hendrik Ranocha, David E. Keyes, and Matteo Parsani. On the robustness and performance of entropy stable collocated discontinuous galerkin methods. Journal of Computational Physics, 426:109891, February 2021.
[4] Vikram Singh and Steven Frankel. On the use of split forms and wall modeling to enable accurate high-reynolds number discontinuous galerkin simulations on body-fitted unstructured grids. Computers & Fluids, 208:104616, 2020.
[5] Magnus Svärd and Jan Nordström. Review of summation-by-parts schemes for initial–boundary-value problems. Journal of Computational Physics, 268:17 – 38, 2014.
[6] Antony Jameson. Formulation of Kinetic Energy Preserving Conservative Schemes for Gas Dynamics and Direct Numerical Simulation of One-Dimensional Viscous Compressible Flow in a Shock Tube Using Entropy and Kinetic Energy Preserving Schemes. Journal of Scientific Computing, 34(2):188–208, February 2008.
[7] Philippe G. LeFloch and Christian Rohde. High-Order Schemes, Entropy Inequalities, and Nonclassical Shocks. SIAM Journal on Numerical Analysis, 37(6):2023–2060, January 2000.
[8] Christopher A. Kennedy and Andrea Gruber. Reduced aliasing formulations of the convective terms within the navier–stokes equations for a compressible fluid. Journal of Computational Physics, 227(3):1676 – 1700, 2008.
[9] Yohei Morinishi. Skew-symmetric form of convective terms and fully conservative finite difference schemes for variable density low-mach number flows. Journal of Computational Physics, 229(2):276 – 300, 2010.
[10] G.A. Blaisdell, E.T. Spyropoulos, and J.H. Qin. The effect of the formulation of nonlinear terms on aliasing errors in spectral methods. Applied Numerical Mathematics, 21(3):207–219, July 1996.
[11] F. Ducros, F. Laporte, T. Soulères, V. Guinot, P. Moinat, and B. Caruelle. High-order fluxes for conservative skew-symmetric-like schemes in structured meshes: Application to compressible flows. Journal of Computational Physics, 161(1):114–139, June 2000.
[12] G. Coppola, F. Capuano, S. Pirozzoli, and L. de Luca. Numerically stable formulations of convective terms for turbulent compressible flows. Journal of Computational Physics, 382:86–104, 2019.
[13] Travis C. Fisher, Mark H. Carpenter, Jan Nordström, Nail K. Yamaleev, and Charles Swanson. Discretely conservative finite-difference formulations for nonlinear conservation laws in split form: Theory and boundary conditions. Journal of Computational Physics, 234:353 – 375, 2013.
[14] Gregor J. Gassner, Magnus Svärd, and Florian J. Hindenlang. Stability issues of entropy-stable and/or split-form high-order schemes, 2020.
[15] Hendrik Ranocha and Gregor J Gassner. Preventing pressure oscillations does not fix local linear stability issues of entropy-based split-form high-order schemes, 2020.
[16] Nao Shima, Yuichi Kuya, Yoshiharu Tamaki, and Soshi Kawai. Preventing spurious pressure oscillations in split convective form discretization for compressible flows. Journal of Computational Physics, 427:110060, February 2021. bibtex: Shima2021.
[17] Yuichi Kuya, Kosuke Totani, and Soshi Kawai. Kinetic energy and entropy preserving schemes for compressible flows by split convective forms. Journal of Computational Physics, 375:823–853, December 2018.
[18] Michael Schlottke-Lakemper, Gregor J. Gassner, Hendrik Ranocha, and Andrew R. Winters. Trixi.jl, 2021.
[19] Mark Carpenter and Christopher Kennedy. Fourth-order 2n-storage runge-kutta schemes. Nasa reports TM, 109112,, 07 1994.
[20] Vikram Singh, Steven Frankel, and Jan Nordström. Impact of wall modeling on kinetic energy stability for the compressible navier-stokes equations. Computers & Fluids, 220:104870, April 2021.
[21] Seth C. Spiegel, James R. DeBonis, and H.T. Huynh. Overview of the NASA glenn flux reconstruction based high-order unstructured grid code. In 54th AIAA Aerospace Sciences Meeting. American Institute of Aeronautics and Astronautics, January 2016.
[22] Gregor J. Gassner, Andrew R. Winters, and David A. Kopriva. Split form nodal discontinuous galerkin schemes with summation-by-parts property for the compressible euler equations. Journal of Computational Physics, 327:39–66, December 2016.
[23] Andrea D. Beck, Thomas Bolemann, David Flad, Hannes Frank, Gregor J. Gassner, Florian Hindenlang, and Claus-Dieter Munz. High-order discontinuous galerkin spectral element methods for transitional and turbulent flow simulations. International Journal for Numerical Methods in Fluids, 76(8):522–548, August 2014.
[24] David Flad and Gregor Gassner. On the use of kinetic energy preserving DG-schemes for large eddy simulation. Journal of Computational Physics, 350:782–795, December 2017.
[25] Andrew R. Winters, Rodrigo C. Moura, Gianmarco Mengaldo, Gregor J. Gassner, Stefanie Walch, Joaquim Peiro, and Spencer J. Sherwin. A comparative study on polynomial dealiasing and split form discontinuous galerkin schemes for under-resolved turbulence computations. Journal of Computational Physics, 372:1–21, November 2018.
[26] James DeBonis. Solutions of the Taylor-Green vortex problem using high-resolution explicit finite difference methods. Aerospace Sciences Meetings, Jan 2013.
[27] John Kim, Parviz Moin, and Robert Moser. Turbulence statistics in fully developed channel flow at low reynolds number. Journal of Fluid Mechanics, 177(-1):133, April 1987.
[28] Myoungkyu Lee and Robert D. Moser. Direct numerical simulation of turbulent channel flow up to $\mathit{Re}_{{\it\tau}}\approx 5200$ . Journal of Fluid Mechanics, 774:395–415, 2015.
[29] Saul Abarbanel and David Gottlieb. Optimal time splitting for two- and three-dimensional navier-stokes equations with mixed derivatives. Journal of Computational Physics, 41(1):1 – 33, 1981.