A Kinetic Energy Preserving DG Scheme Based on Gauss–Legendre Points

Abstract

In the context of numerical methods for conservation laws, not only the preservation of the primary conserved quantities can be of interest, but also the balance of secondary ones such kinetic energy in case of the Euler equations of gas dynamics. In this work, we construct a kinetic energy preserving discontinuous Galerkin method on Gauss–Legendre nodes based on the framework of summation-by-parts operators. For a Gauss–Legendre point distribution, boundary terms require special attention. In fact, stability problems will be demonstrated for a combination of skew-symmetric and boundary terms that disagrees with exclusively interior nodal sets. We will theoretically investigate the required form of the corresponding boundary correction terms in the skew-symmetric formulation leading to a conservative and consistent scheme. In numerical experiments, we study the order of convergence for smooth solutions, the kinetic energy balance and the behaviour of different variants of the scheme applied to an acoustic pressure wave and a viscous shock tube. Using Gauss–Legendre nodes results in a more accurate approximation in our numerical experiments for viscous compressible flow. Moreover, for two-dimensional decaying homogeneous turbulence, kinetic energy preservation yields a better representation of the energy spectrum.

Appendix

In the following, the proofs of Lemma 2 in Sect. 3 and Lemma 3 in Sect. 5.1 are given.

1.1 Proof of Lemma 2

For the assertion of conservativity, we multiply the scheme (25) from left by $\underline{1}^T\underline{\underline{M}}$. Thus, the rate of change of the cell mean within a cell $[x_i,x_{i+1}]$ multiplied by the cell size $\varDelta x$ is given by

$$\begin{aligned} \underline{1}^T\underline{\underline{M}}\frac{\varDelta x_i}{2}\,\frac{d}{dt} \underline{u}^i= & {} - \underline{1}^T\underline{\underline{M}}\underline{\underline{D}}\underline{f}^i - \underline{1}^T\underline{\underline{M}}\left[ -\underline{\underline{D}}\underline{\underline{\alpha }}^i\underline{\beta }^i + \underline{\underline{\alpha }}^i\underline{\underline{D}}\underline{\beta }^i + \underline{\underline{\beta }}^i\underline{\underline{D}}\underline{\alpha }^i\right] \\&+\, \underline{1}^T[(f^i_{h} + \underline{\underline{\alpha }}^i\beta ^i_h + k^i_h)\underline{L}]_{-1}^{1} - \underline{1}^T[f^{*,i}\underline{L}]_{-1}^{1} -\underline{1}^T[\underline{\underline{\alpha }}^i\beta ^{*,i}\underline{L}]_{-1}^{1}\\&+\, \underline{1}^T\left( k_i^{*,+}\underline{L}(-1) - k_i^{*,-}\underline{L}(1)\right) \\= & {} - \underline{1}^T(\underline{\underline{B}}-\underline{\underline{D}}^T\underline{\underline{M}}) \underline{f}^i\\&- \,\underline{1}^T\left[ -(\underline{\underline{B}}-\underline{\underline{D}}^T\underline{\underline{M}})\underline{\underline{\alpha }}^i\underline{\beta }^i + \underline{\underline{\alpha }}^i(\underline{\underline{B}}-\underline{\underline{D}}^T\underline{\underline{M}})\underline{\beta }^i + \underline{\underline{\beta }}^i\underline{\underline{M}}\underline{\underline{D}}\underline{\alpha }^i\right] \\&+ \,\underline{1}^T[(f^i_{h} + \underline{\underline{\alpha }}\beta ^i_h + k^i_h)\underline{L}]_{-1}^{1} -\underline{1}^T[f^{*,i}\underline{L}]_{-1}^{1} -\underline{1}^T[\underline{\underline{\alpha }}^i\beta ^{*,i}\underline{L}]_{-1}^{1}\\&+ \,\underline{1}^T\left( k_i^{*,+}\underline{L}(-1) - k_i^{*,-}\underline{L}(1)\right) , \end{aligned}$$

using the SBP property $\underline{\underline{M}}\underline{\underline{D}}= \underline{\underline{B}}- \underline{\underline{D}}^T\underline{\underline{M}}$.

Furthermore, the interpolation of boundary values and the corresponding boundary terms are encoded in the matrix $\underline{\underline{B}}$. Thus, due to Lemma 1 we can replace the corresponding terms by $\underline{\underline{B}}\underline{f}^i = [f^i_{h}\underline{L}]_{-1}^{1},\ \underline{\underline{B}}\underline{\underline{\alpha }}^i\underline{\beta }^i = [(\alpha \beta )^i_{h}\underline{L}]_{-1}^{1}$ and $(\underline{\alpha }^i)^T\underline{\underline{B}}\underline{\beta }^i = \underline{1}^T\underline{\underline{\alpha }}^i[\beta ^i_{h}\underline{L}]_{-1}^{1}$. Since this cancels out the terms containing $f^i_h$ and $\underline{\underline{\alpha }}^i\beta ^i_h$, it holds that

$$\begin{aligned} \underline{1}^T\underline{\underline{M}}\frac{\varDelta x_i}{2}\,\frac{d}{dt} \underline{u}^i= & {} (\underline{\underline{D}}\underline{1})^T\underline{\underline{M}}\underline{f}^i - \left[ (\underline{\underline{D}}\underline{1})^T\underline{\underline{M}}\underline{\underline{\alpha }}^i\underline{\beta }^i - (\underline{\alpha }^i)^T\underline{\underline{D}}^T\underline{\underline{M}}\underline{\beta }^i + (\underline{\beta }^i)^T\underline{\underline{M}}\underline{\underline{D}}\underline{\alpha }^i\right] \\&+\, \underline{1}^T[\left( (\alpha \beta )^i_{h}+k^i_h\right) \underline{L}]_{-1}^{1}\\&- \,\underline{1}^T[f^{*,i}\underline{L}]_{-1}^{1} - \underline{1}^T[\underline{\underline{\alpha }}^i\beta ^{*,i}\underline{L}]_{-1}^{1} + \underline{1}^T\left( k_i^{*,+}\underline{L}(-1) - k_i^{*,-}\underline{L}(1)\right) . \end{aligned}$$

Now, discrete differentiation yields $\underline{\underline{D}}\underline{1} = \underline{0}$ and obviously, $-\underline{\alpha }^T\underline{\underline{D}}^T\underline{\underline{M}}\underline{\beta }+ \underline{\beta }^T\underline{\underline{M}}\underline{\underline{D}}\underline{\alpha }= \underline{0}$ holds. This reduces the rate of change of mass contained in the specific cell to

$$\begin{aligned} \underline{1}^T\underline{\underline{M}}\frac{\varDelta x_i}{2}\,\frac{d}{dt} \underline{u}^i= & {} \underline{1}^T[\left( (\alpha \beta )^i_{h} + k^i_h\right) \underline{L}]_{-1}^{1} -\underline{1}^T[f^{*,i}\underline{L}]_{-1}^{1} -\underline{1}^T[\underline{\underline{\alpha }}^i\beta ^{*,i}\underline{L}]_{-1}^{1}\\&+\, \underline{1}^T\left( k_i^{*,+}\underline{L}(-1) - k_i^{*,-}\underline{L}(1)\right) \end{aligned}$$

For conservativity, the remaining terms may only contain interface contributions. Hence, the balance of flux contributions at cell interfaces has to be investigated. We consider the interface with index $(i-1,i)$ between two cells $[x_{i-1},x_i]$ and $[x_i,x_{i+1}]$. As the interpolation property yields $\underline{1}^T\underline{L}(\xi ) = 1$, we obtain the following fluxes over cell interfaces. Denoting by $C^{*,i-1}_{i-1,i}$ the corresponding flux from left to right based on the values in the left cell, neglecting the term $f_{i-1,i}^*$, we have

$$\begin{aligned} C^{*,i-1}_{i-1,i} = -(\alpha \beta )^{i-1}_{h}(1) - k^{i-1}_h(1) + \alpha ^{i-1}_h(1)\beta _{i-1,i}^*+ k_{i-1}^{*,-}, \end{aligned}$$

whereas the analogous quantity $C^{*,i}_{i-1,i}$ based on the right cell is given by

$$\begin{aligned} C^{*,i}_{i-1,i} = -(\alpha \beta )^i_{h}(-1) - k^i_h(-1) + \alpha _h^i(-1)\beta _{i-1,i}^*+ k_{i}^{*,+} . \end{aligned}$$

Hence, conservativity precisely requires

$$\begin{aligned} C^{*,i-1}_{i-1,i} = C^{*,i}_{i-1,i} =: C^{*}_{i-1,i}, \end{aligned}$$

yielding the first assertion.

As in the second assertion, we now assume that $k^i_h$ only depends on a combination of the interior values $\underline{\alpha }^i,\underline{\beta }^i$ and that $k_i^{*,+}$ only depends on $\beta _{i-1,i}^*$ and $\alpha ^i_h(-1)$. If we then modify any of the input values other than $\underline{\alpha }^i,\underline{\beta }^i,\beta ^*_{i-1,i}$, the value of $C^{*,i}_{i-1,i}$ does not change and thus $C^{*}_{i-1,i}$ remains constant as well. If we furthermore assume that $k_i^{*,-}$ only depends on $\beta _{i,i+1}^*$ and $\alpha ^i_h(1)$, modifying $\underline{\alpha }^i$ or $\underline{\beta }^{i}$ does not change $C^{*,i-1}_{i-1,i}$ and hence does not change $C^{*}_{i-1,i}$. Therefore, $C^{*}_{i-1,i}$ may only depend on $\beta ^*_{i-1,i}$ which proves assertion 2.

Now, consistency in the finite volume sense refers to a consistent numerical flux. With $k^i_h$ and $k_i^{*,\pm }$ set as in the third assertion, the scheme (25) for $N=0$ reduces to

$$\begin{aligned} \varDelta x_i\frac{d}{dt} u^i= & {} [f^{*,i}]_{-1}^{1} -[\alpha ^i\beta ^{*,i}]_{-1}^{1} + \left( k_i^{*,+} - k_i^{*,-}\right) = [f^{*,i}]_{-1}^{1} - \left( C_{i,i+1}^{*} - C_{i-i,i}^{*}\right) . \end{aligned}$$

Since $f^*$ is consistent to f with $f_{i-1,i}^*(u,u) = f(u)$ and $\beta ^*$ is consistent to $\beta$, this above finite volume scheme is consistent if and only if

$$\begin{aligned} f^*_{i-1,i}(u,u) + C^*_{i-1,i}(\beta _{i-1,i}^{*}(u,u)) = f(u) + C^*_{i-1,i}(\beta (u)) = f(u), \end{aligned}$$

for any interface index $(i-1,i)$ and any value of the conserved variable u. Hence, we have the requirement $C_{i-1,i}(\beta (u)) = 0$, which yields the last assertion.

1.2 Proof of Lemma 3

The columns of the matrices $\underline{\underline{B_\xi }}= \underline{\underline{B}}\otimes \underline{\underline{M}}_{1D}$ and $\underline{\underline{B_\eta }}= \underline{\underline{M}}_{1D}\otimes \underline{\underline{B}}$ can be related to the Lagrange polynomials as follows. Let $\nu (i,j)=(i-1)(N+1)+j$. From the proof of Lemma 1 in Sect. 2, we recall that $\underline{\underline{B}}$ has the entries $B_{jk}=[L_jL_k]^1_{-1}$ and $\underline{\underline{M}}_{1D}$ is diagonal with entries $\omega _j$. Therefore, the $\nu$th columns of the above matrices $\underline{\underline{B_\xi }}$ and $\underline{\underline{B_\eta }}$ are given by $\left. \underline{B_\xi }\right. _{\nu (i,j)} = [\underline{L}(\xi )L_i(\xi )]^1_{-1}\otimes \omega _j\underline{\epsilon }_j$ and $\left. \underline{B_\eta }\right. _{\nu (i,j)} = \omega _i\underline{\epsilon }_i\otimes [\underline{L}(\xi )L_j(\xi )]^1_{-1}$, respectively, where $\underline{\epsilon }_j$ denotes the jth unit vector.

Therefore, we obtain

$$\begin{aligned} \underline{\underline{B_\xi }}\underline{g}^\xi= & {} \sum _{i=1}^{N+1}\sum _{j=1}^{N+1}\left. \underline{B_\xi }\right. _{\nu (i,j)} g^\xi _{\nu (i,j)} = \sum _{i=1}^{N+1}\sum _{j=1}^{N+1}\left( [\underline{L}(\xi )L_i(\xi )]^1_{-1}\otimes \omega _j\underline{\epsilon }_j\right) g^\xi _{\nu (i,j)}\\= & {} \left[ \sum _{j=1}^{N+1}\omega _j\left( \sum _{i=1}^{N+1}L_i(\xi ) g^\xi _{\nu (i,j)}\right) \underline{L}(\xi )\otimes \underline{\epsilon }_j\right] ^1_{-1} = \left[ \sum _{j=1}^{N+1}\omega _jg^\xi _h(\xi ,\eta _j)\underline{L}(\xi )\otimes \underline{L}(\eta _j)\right] ^1_{-1}. \end{aligned}$$

Since for the normal vectors and Gauss nodes on edges $e=2$ and $e=4$ we have $n^2_\xi =1, n^4_\xi =-1$ and $\xi ^2_j = 1, \xi ^4_j = -1, j=1,\ldots ,N+1$, it holds that

$$\begin{aligned} \underline{\underline{B_\xi }}\underline{g}^\xi= & {} \sum _{e=2,4}\sum _{j=1}^{N+1}n^e_\xi \omega _jg^\xi _h(\xi ^e_j,\eta ^e_j)\underline{L}(\xi ^e_j)\otimes \underline{L}(\eta ^e_j) \end{aligned}$$

With analogous arguments we obtain

$$\begin{aligned} \underline{\underline{B_\eta }}\underline{g}^\eta= & {} \left[ \sum _{i=1}^{N+1}\omega _i\underline{\epsilon }_i\otimes \underline{L}(\eta )\sum _{j=1}^{N+1}L_j(\eta ) g^\eta _{\nu (i,j)}\right] ^1_{-1}\\= & {} \sum _{e=1,3}\sum _{i=1}^{N+1}n^e_\eta \omega _ig^\eta _h(\xi ^e_i,\eta ^e_i)\underline{L}(\xi _i^e)\otimes \underline{L}(\eta ^e_i). \end{aligned}$$

Summing up $\underline{\underline{B_\xi }}\underline{g}^\xi + \underline{\underline{B_\eta }}\underline{g}^\eta$ and considering that $n^1_\xi =n^2_\eta =n^3_\xi =n^4_\eta =0$ we obtain the assertion.