Stability of Wall Boundary Condition Procedures for Discontinuous Galerkin Spectral Element Approximations of the Compressible Euler Equations

TL;DR

This paper analyzes the stability of wall boundary condition procedures in discontinuous Galerkin spectral element methods for the compressible Euler equations, highlighting conditions for energy and entropy stability.

Contribution

It provides a comprehensive stability analysis of different boundary procedures, including new insights into their dissipation properties and robustness.

Findings

Exact upwind and Lax-Friedrichs solvers are energy dissipative.

Standard Riemann solvers like Lax-Friedrichs, HLL, HLLC are entropy stable.

Roe flux is entropy stable under certain conditions.

Abstract

We perform a linear and entropy stability analysis for wall boundary condition procedures for discontinuous Galerkin spectral element approximations of the compressible Euler equations. Two types of boundary procedures are examined. The first defines a special wall boundary flux that incorporates the boundary condition. The other is the commonly used reflection condition where an external state is specified that has an equal and opposite normal velocity. The internal and external states are then combined through an approximate Riemann solver to weakly impose the boundary condition. We show that with the exact upwind and Lax-Friedrichs solvers the approximations are energy dissipative, with the amount of dissipation proportional to the square of the normal Mach number. Standard approximate Riemann solvers, namely Lax-Friedrichs, HLL, HLLC are entropy stable. The Roe flux is entropy stable under certain conditions. An entropy conserving flux with an entropy stable dissipation term (EC-ES) is also presented. The analysis gives insight into why these boundary conditions are robust in that they introduce large amounts of energy or entropy dissipation when the boundary condition is not accurately satisfied, e.g. due to an impulsive start or under resolution.

Figures

Equations

$${\mathbf{u}}_t+\sum _{i=1}^3\frac{\partial {\mathbf{f}}_i}{\partial x_i}=0.$$

$$\mathbf{u}={\left[\varrho \varrho \text{spacevec}vE\right]}^T={\left[\varrho \varrho v_1\varrho v_2\varrho v_{3}E\right]}^T.$$

$$\mathbf{f}_{1}=\left[{\begin{array}[]{*{20}c}{\varrho v_{1}}\\ {\varrho v_{1}^{2}+p}\\ {\varrho v_{1}\,v_{2}}\\ {\varrho v_{1}\,v_{3}}\\ {(E+p)v_{1}}\\ \end{array}}\right]\quad\mathbf{f}_{2}=\left[{\begin{array}[]{*{20}c}{\varrho v_{2}}\\ {\varrho v_{2}\,v_{1}}\\ {\varrho v_{2}^{2}+p}\\ {\varrho v_{2}\,v_{3}}\\ {(E+p)v_{2}}\\ \end{array}}\right]\quad\mathbf{f}_{3}=\left[{\begin{array}[]{*{20}c}{\varrho v_{3}}\\ {\varrho v_{3}\,v_{1}}\\ {\varrho v_{3}\,v_{2}}\\ {\varrho v_{3}^{2}+p}\\ {(E+p)v_{3}}\\ \end{array}}\right],$$

$$p=(\gamma -1)\left(E-\frac{1}{2}\varrho {\Vert \text{spacevec}v\Vert }^2\right),$$

$$\text{accentset}\leftrightarrow \mathbf{f}={\left[{\mathbf{f}}_1{\mathbf{f}}_2{\mathbf{f}}_3\right]}^T,$$

$${\mathbf{u}}_t+\text{spacevec}{\nabla }_x\cdot \text{accentset}\leftrightarrow \mathbf{f}=0.$$

$$a=\sqrt{\frac{\gamma -1}{\gamma }}\bar{c},b=\frac{\bar{c}}{\sqrt{\gamma }},\bar{c}=\sqrt{\frac{\gamma \bar{p}}{\bar{\varrho }}},$$

$$\mathbf{u}={\left[{\varrho }^{\prime }{v}_1{v}_2{v}_3{p}^{\prime }\right]}^T,$$

$${\varrho }^{\prime }=b\frac{\tilde{\varrho }}{\bar{\varrho }},p^{\prime }=\frac{1}{\bar{\varrho }a}\tilde{p}-\frac{1}{\sqrt{\gamma -1}}{\varrho }^{\prime },$$

$${\mathbf{f}}_i={\underline{A}}_i\mathbf{u},\text{accentset}\leftrightarrow \mathbf{f}=\text{spacevec}\underline{A}\mathbf{u}=\left({\underline{A}}_1\hat{x}+{\underline{A}}_2\hat{y}+{\underline{A}}_3\hat{z}\right)\mathbf{u},$$

$${\underline{A}_{1}}=\left[{\begin{array}[]{*{20}{c}}\bar{v}_{1}&b&0&0&0\\ b&\bar{v}_{1}&0&0&a\\ 0&0&\bar{v}_{1}&0&0\\ 0&0&0&\bar{v}_{1}&0\\ 0&a&0&0&\bar{v}_{1}\end{array}}\right],\quad{\underline{A}_{2}}=\left[{\begin{array}[]{*{20}{c}}\bar{v}_{2}&0&b&0&0\\ 0&\bar{v}_{2}&0&0&0\\ b&0&\bar{v}_{2}&0&a\\ 0&0&0&\bar{v}_{2}&0\\ 0&0&a&0&\bar{v}_{2}\end{array}}\right],\quad{\underline{A}_{3}}=\left[{\begin{array}[]{*{20}{c}}\bar{v}_{3}&0&0&b&0\\ 0&\bar{v}_{3}&0&0&0\\ 0&0&\bar{v}_{3}&0&0\\ b&0&0&\bar{v}_{3}&a\\ 0&0&0&a&\bar{v}_{3}\end{array}}\right]$$

$$⟨\mathbf{v},\mathbf{w}⟩=\underset{\Omega }{\int }{\mathbf{v}}^T\mathbf{w}dxdydz,⟨\text{accentset}\leftrightarrow \mathbf{f},\text{accentset}\leftrightarrow \mathbf{g}⟩=\underset{\Omega }{\int }\sum _{i=1}^3{\mathbf{f}}_i^T{\mathbf{g}}_{i}dxdydz.$$

$$⟨\text{spacevec}{\nabla }_x\cdot \text{accentset}\leftrightarrow \mathbf{f},\mathbf{u}⟩=⟨\text{spacevec}{\nabla }_x\cdot \left(\text{spacevec}\underline{A}\mathbf{u}\right),\mathbf{u}⟩=⟨{\nabla }_x\mathbf{u},\text{accentset}\leftrightarrow \mathbf{f}⟩.$$

$$⟨\text{spacevec}{\nabla }_x\cdot \text{accentset}\leftrightarrow \mathbf{f},\mathbf{u}⟩=\frac{1}{2}{\int }_{\partial \Omega }{\mathbf{u}}^T\text{accentset}\leftrightarrow \mathbf{f}\cdot \text{spacevec}ndS,$$

$$\frac{d}{dt}{\left|\left|\mathbf{u}\right|\right|}^2=-{\int }_{\partial \Omega }{\mathbf{u}}^T\text{accentset}\leftrightarrow \mathbf{f}\cdot \text{spacevec}ndS.$$

$${\mathbf{u}}^T\text{accentset}\leftrightarrow \mathbf{f}\cdot \text{spacevec}n={\mathbf{u}}^T\left(\text{spacevec}\underline{A}\cdot \text{spacevec}n\right)\mathbf{u}=2\left(\varrho b+ap\right)v_n+(\text{spacevec}\bar{v}\cdot \text{spacevec}n)({\varrho }^2+|\text{spacevec}v{|}^2+p^2),$$

$$\frac{d}{dt}{\left|\left|\mathbf{u}\right|\right|}^2=0,$$

$$s(\mathbf{u})=-\frac{\varrho \varsigma }{(\gamma -1)},$$

$$\text{spacevec}{f}^{\varsigma }(\mathbf{u})=\text{spacevec}vs=-\frac{\varrho \varsigma \text{spacevec}v}{(\gamma -1)}.$$

$$\mathbf{w}=\frac{\partial s(\mathbf{u})}{\partial\mathbf{u}}=\left[{\begin{array}[]{*{20}{c}}\frac{\gamma-\varsigma}{\gamma-1}-\beta||\spacevec{v}||^{2},\\ {2\beta\spacevec{v}}\\ {-2\beta}\end{array}}\right]\,,\quad\beta=\frac{\varrho}{2p}\,.$$

$${\mathbf{w}}^T{\mathbf{u}}_t={\left(\frac{\partial s}{\partial \mathbf{u}}\right)}^T{\mathbf{u}}_t=s_t(\mathbf{u}),{\mathbf{w}}^T\text{spacevec}{\nabla }_x\cdot \text{accentset}\leftrightarrow \mathbf{f}=\text{spacevec}{\nabla }_x\cdot \text{spacevec}{f}^{\varsigma }.$$

$$⟨\mathbf{w}(\mathbf{u}),{\mathbf{u}}_t⟩+⟨\mathbf{w}(\mathbf{u}),\text{spacevec}{\nabla }_x\cdot \text{accentset}\leftrightarrow \mathbf{f}⟩=0.$$

$$⟨s_t(\mathbf{u}),1⟩=-⟨\text{spacevec}{\nabla }_x\cdot \text{spacevec}{f}^{\varsigma },1⟩=-\underset{\partial \Omega }{\int }\left(\text{spacevec}{f}^{\varsigma }\cdot \text{spacevec}n\right)\mathrm{dS}$$

$$-\text{spacevec}{f}^{\varsigma }\cdot \text{spacevec}n=\frac{\varrho \varsigma }{(\gamma -1)}{v}_n=0.$$

$$\mathcal{J}{\mathbf{u}}_t+\text{spacevec}{\nabla }_{\xi }\cdot \text{accentset}\leftrightarrow \tilde{\mathbf{f}}=0,$$

$${\mathbf{U}}^T\left\{{\tilde{\mathbf{F}}}^{\ast }-\frac{1}{2}\text{accentset}\leftrightarrow \tilde{\mathbf{F}}\cdot \hat{n}\right\}\ge 0,$$

$$B_L\equiv {\mathbf{U}}^T\left\{{\mathbf{F}}^{\ast }-\frac{1}{2}\text{accentset}\leftrightarrow \mathbf{F}\cdot \text{spacevec}n\right\}\ge 0,$$

$$B_{NL}\equiv {\mathbf{W}}^T\left({\mathbf{F}}^{\ast }-\left(\text{accentset}\leftrightarrow \mathbf{F}\cdot \text{spacevec}n\right)\right)+\left(\text{spacevec}{F}^{\varsigma }\cdot \text{spacevec}n\right)\ge 0,$$

$$\text{accentset}\leftrightarrow \mathbf{F}\cdot \text{spacevec}n=\text{spacevec}\underline{A}\cdot \text{spacevec}n\mathbf{U}={\left[b{V}_n{n}_{1}Q{n}_{2}Q{n}_{3}Qa{V}_n\right]}^T,$$

$${\mathbf{F}}^{\ast }=\text{spacevec}\underline{A}\cdot \text{spacevec}n{\mathbf{U}}^{\ast }={\left[b{V}_n^{\ast }{n}_1{Q}^{\ast }{n}_2{Q}^{\ast }{n}_3{Q}^{\ast }a{V}_n^{\ast }\right]}^T.$$

Full text

1 Introduction

The ingredients for a reliable numerical method for the approximation of partial differential equations, e.g. one that will not blow up, include stable inter-element and physical boundary condition implementations. The recognition that the discontinuous Galerkin spectral element method (DGSEM) with Gauss-Lobatto quadratures satisfies a summation-by-parts (SBP) operators Gassner:2016ye ; gassner2010 has allowed for the analysis of these schemes and to connect them with penalty collocation and SBP finite difference schemes. For instance, in Gassner2018 , we showed that a split form approximation of the compressible Navier-Stokes equations was both linearly and entropy stable provided that the boundary conditions were properly imposed.

The importance of stable boundary condition procedures for hyperbolic equations has long been studied, especially in relation to finite difference methods, e.g. skew_sbp2 ; nordstroem:2006 ; Nordstrom:2016jk . Only recently have they been studied for discontinuous Galerkin approximations. In PARSANI201588 , the authors showed that the reflection approach is stable when using an entropy conserving flux and an additional entropy stable dissipation term (EC-ES). In CHEN2017427 , the authors show that the reflection condition is stable if the numerical flux is either the Godunov or HLL flux.

In this paper, we analyze both the linear and entropy stability of two types of commonly used wall boundary condition procedures used with the DGSEM applied to the compressible Euler equations. In both cases, wall boundary conditions are implemented through a numerical flux. The boundary condition might be implemented through a special wall numerical flux that includes the boundary condition, or a fictitious external state applied to a Riemann solver approximation. We show how to construct special wall numerical fluxes that are stable, and study the behavior of the approximations. In particular, we show that the use of Riemann solvers at the boundaries introduce numerical dissipation in an amount that depends on the size of the normal Mach number at the wall.

2 The compressible Euler Equations and the Wall Boundary Condition

We write the Euler equations as

[TABLE]

The state vector contains the conservative variables

[TABLE]

In standard form, the components of the advective fluxes are

[TABLE]

Here, $\varrho,\,\spacevec{v}=(v_{1},v_{2},v_{3})^{T},\,p,\,E$ are the mass density, fluid velocities, pressure and total energy. We close the system with the ideal gas assumption, which relates the total energy and pressure

[TABLE]

where $\gamma$ denotes the adiabatic coefficient. For a compact notation that simplifies the analysis, we define block vectors (with the double arrow)

[TABLE]

so that the system of equations can be written in the compact form

[TABLE]

The linear Euler equations are derived by linearizing about a constant mean state $(\bar{\varrho},\bar{v}_{1},\bar{v}_{2},\bar{v}_{3},\bar{p})$. We follow Nordstrom:2005qy for the symmetrization of the linearized equations, with the constants

[TABLE]

where $\bar{c}$ is the sound speed of the constant mean state. The state variables become

[TABLE]

where $\spacevec{v}$ is the velocity perturbation from the mean state, and we introduce

[TABLE]

which depend on the density and pressure perturbations $\tilde{\varrho},\tilde{p}$. The flux vectors are

[TABLE]

where Nordstrom:2005qy

[TABLE]

are constant symmetric matrices.

The linear equations have the property that the $L^{2}$ norm of the solution over a domain $\Omega$ is bounded by terms of the boundary data on $\partial\Omega$, only. Let

[TABLE]

represent the $L^{2}$ inner product of two state vectors $\mathbf{v}$ and $\mathbf{w}$ and two block vectors $\accentset{\,\leftrightarrow}{{\mathbf{f}}}$ and $\accentset{\,\leftrightarrow}{{\mathbf{g}}}$, respectively. Since the coefficient matrices are constant the product rule and symmetry of $\spacevec{\underline{A}}$ implies

[TABLE]

Then it follows from Gauss’ law (integration by parts) that

[TABLE]

where $\spacevec{n}$ is the outward normal to the surface of $\Omega$. The norm of the solution therefore satisfies

[TABLE]

Replacing the boundary terms by boundary conditions leads to a bound on the solution in terms of the boundary data. The argument of the boundary integral on the right of (15) is

[TABLE]

where $v_{n}$ is the wall normal velocity, $v_{n}=\spacevec{v}\cdot\spacevec{n}$. Note that here, the mean flow must be chosen such that the normal flow vanishes at the wall boundary $\spacevec{\bar{v}}\cdot\spacevec{n}=0$, so that the boundary condition makes physical sense.

Therefore, with the no penetration wall condition $v_{n}=0$ applied,

[TABLE]

and the (energy) norm of the solution is bounded for all time by its initial value.

The nonlinear equations, on the other hand, satisfy a bound on the entropy that depends only on the boundary data. For what follows, we assume that the solution is smooth so that we don’t have to consider entropy generated at shock waves. We introduce the entropy density (scaled with $(\gamma-1)$ for convenience) as

[TABLE]

where $\varsigma=\ln(p)-\gamma\ln(\varrho)$ is the physical entropy. (The minus sign is conventional in the theory of hyperbolic conservation laws to ensure a decreasing entropy function.) The entropy flux for the Euler equations is

[TABLE]

Finally the entropy variables are

[TABLE]

The entropy pair contracts the solution and fluxes, meaning that it satisfies the relations

[TABLE]

When we multiply (6) with the entropy variables and integrate over the domain,

[TABLE]

Next we use the properties of the entropy pair to contract 22 and use integration by parts to get

[TABLE]

showing that, in the continuous case, the total entropy in the domain can only change via the boundary conditions.

In the case of a zero-mass flux boundary condition, with $v_{n}=\spacevec{v}\cdot\spacevec{n}=0$, the entropy is not changed by the slip-wall boundary condition, since

[TABLE]

3 Stability bounds for the DGSEM

The DGSEM is described in detail in Gassner2018 and elsewhere Black1999 ; Kopriva:2009nx . We will only quickly summarize the approximation here. The domain, $\Omega$ is subdivided into non-overlapping, conforming, hexahedral elements. Each element is mapped to the reference element $E=[-1,1]^{3}$. Associated with the transformation from the reference element is a set of contravariant coordinate vectors, $\spacevec{a}^{i}$, and transformation Jacobian, $\mathcal{J}$. The equations (6) transform to another conservation law on the reference element as

[TABLE]

where $\accentset{\,\leftrightarrow}{\tilde{\mathbf{f}}}$ is the contravariant flux vector with components $\tilde{\mathbf{f}}^{i}=\mathcal{J}\spacevec{a}^{i}\cdot\accentset{\,\leftrightarrow}{{\mathbf{f}}}$.

The approximation of (25) proceeds as follows: A weak form is created by taking the inner product of the equation with a test function. The Gauss law is applied to the divergence term to separate the boundary from the interior contributions. The resulting weak form is then approximated: The solution vector is approximated by a polynomial of degree $N$ interpolated at the Legendre Gauss-Lobatto points. In the following, we will represent the true continuous solutions by lower case letter. Upper case letters will denote their polynomial approximations, except for the density, where the approximation is denoted by $\rho$. The volume fluxes are replaced by two-point numerical fluxes. In the linear case, the two point fluxes are immediately relatable to a split form of the equations. Integrals are replaced by Legendre-Gauss-Lobatto quadratures. Finally, the boundary fluxes are replaced by a numerical flux. See Gassner2018 and Gassner:2013uq for details.

The result is an approximation that is energy stable for the linearized equations if at every quadrature point along a physical boundary the numerical flux $\tilde{\mathbf{F}}^{*}$ satisfies the bound Gassner2018

[TABLE]

where $\accentset{\,\leftrightarrow}{\tilde{\mathbf{F}}}$ is the polynomial interpolation of the contravariant flux from the interior, $\hat{n}$ is the reference space outward normal direction, and $\mathbf{U}$ is the approximation of the state vector. Since the contravariant fluxes are proportional to the normal fluxes Kopriva:2009nx , we can change the condition (26) to

[TABLE]

For entropy stability of the nonlinear equations, the boundary stability condition shown in Gassner2018 is proportional to

[TABLE]

where $\spacevec{F}^{\,\varsigma}$ is the polynomial interpolation of the entropy flux, $\spacevec{f}^{{\,\varsigma}}$, and $\mathbf{W}$ is the interpolation of the entropy variables.

3.1 Linear Stability of Wall Boundary Condition Approximations

To find linearly stable implementations of the wall condition $v_{n}=0$, one needs only find a numerical flux that satisfies it and the condition (27). For the linear equations, the approximation of the state vector is $\mathbf{U}=[\rho^{\prime}\;\;\spacevec{V}\;\;P^{\prime}]^{T}$ and the normal contravariant flux is proportional to

[TABLE]

where $V_{n}$ is the approximation of the normal velocity at the wall computed from the interior, $Q=b\rho^{\prime}+aP^{\prime}$, and $(n_{1},n_{1},n_{3})$ are the three components of the physical space normal vector, $\spacevec{n}$. The numerical flux can be expressed as

[TABLE]

It then remains only to find $Q^{*}$ so that (27) is satisfied when the normal wall condition $V^{*}_{n}=0$ is applied. When we substitute the fluxes (29) and (30) into (27),

[TABLE]

Substituting the wall boundary condition $V^{*}_{n}=0$ yields the condition on $Q^{*}$ for stability

[TABLE]

Neutral stability is thus ensured if $\rho^{*}$ and $P^{*}$ are computed from the interior, i.e. $\rho^{*}=\rho^{\prime}$, $P^{*}=P^{\prime}$ so that $Q^{*}=Q$.

In practice, the boundary condition is also implemented through the use of a Riemann solver and external state designed to imply the physical boundary condition to construct the numerical boundary flux. The exact upwind ($\varepsilon=1$) normal Riemann flux and the central flux ($\varepsilon=0$) for the linear system of equations is

[TABLE]

where ${{\underline{A}}_{\,n}}\equiv{\spacevec{\underline{A}}}\cdot\spacevec{n}$ is the normal coefficient matrix. The external state is set by by using the interior values of the density and pressure and the negative of the value of the normal velocity,

[TABLE]

For $\varepsilon=0$, using the central (averaged) numerical flux, the interior flux contribution cancels and condition (27) reduces to

[TABLE]

which is neutrally stable, having no additional stabilizing dissipation. We note again, that the mean state for the linearization is chosen such that the normal mean velocity components are zero, resulting in the zeros on the diagonal of ${\underline{A}}_{\,n}$.

Substituting the exact upwind flux where $\varepsilon=1$ into (27) and rearranging,

[TABLE]

where ${\underline{A}}_{\,n}^{-}=\frac{1}{2}\left({\underline{A}}_{\,n}-\left|{\underline{A}}_{\,n}\right|\right)$ is negative semidefinite. The second term is non-negative, depends only on the interior state, and adds stabilizing dissipation. From the matrix absolute value, the dissipation term is

[TABLE]

where ${\operatorname{Ma}}_{n}=V_{n}/\bar{c}$ is the normal Mach number. Stability depends, then, on the value of the first term, which is where the boundary conditions are incorporated through the external state $\mathbf{U}^{\mathrm{ext}}$ written in (34). Then

[TABLE]

Therefore, using the upwind numerical flux, (36) becomes

[TABLE]

as required. The amount of dissipation depends on how far the interior computed normal velocity deviates from zero.

The combination of the reflective state and local Lax-Friedrichs flux is also linearly stable. In that case the exact matrix absolute value is replaced by a diagonal matrix, $\left|\underline{A}_{n}\right|\approx\left|\lambda\right|_{\mathrm{max}}\underline{I}$. The jump term is added to the central (averaged) flux so

[TABLE]

Finally, a dissipative version of the direct numerical flux (30) can be formed by looking at the reflective state approach. For instance, the equivalent to using the Lax-Friedrichs flux is to choose $\rho^{*}=\rho^{\prime}$ and

[TABLE]

Then $Q^{*}=Q+\bar{c}^{3}{\left|\lambda\right|_{\mathrm{max}}}{\operatorname{Ma}}_{n}$ and

[TABLE]

A similar, though more complicated, modified $P^{*}$ can be made to be equivalent to the exact upwind flux.

3.2 Entropy Stability of Wall Boundary Condition Approximations

As in the linear approximation, the wall boundary condition can be imposed for the nonlinear equations either by directly specifying the numerical flux or by computing it through a Riemann solver using a reflection external state that enforces the normal wall condition implicitly. Note that in this section, the discrete variables $(\rho,\spacevec{V},P)$ describe the full nonlinear state.

For the nonlinear equations, we construct the numerical flux for a slip-wall as

[TABLE]

where we imposed $V_{n}=0$ leading to a flux with no mass or energy transfer, and we introduce a wall pressure $P^{*}$, whose value will be chosen to ensure consistency and stability.

After some manipulations, the discrete entropy stability condition (28) becomes

[TABLE]

Therefore if we choose $P^{*}=P$, to be the internal pressure, the boundary flux does not contribute to the total entropy, independent of the inner normal velocity $V_{n}$. A value of $P^{*}$ that leads to a dissipative boundary condition can be found either through exact solution of the Riemann problem at the boundary, or through the use of an external state and an approximate Riemann solver.

Exact solution of the Riemann problem

In VanDerVegt2002_BC a symmetric 1D Riemann problem is exactly solved following Toro torobook , to get the wall pressure $P^{*}$, accounting for the fact that $V_{n}$ never vanishes discretely and therefore the wall pressure should be different from the interior pressure. The exact solution of the 1D Riemann problem reads as

[TABLE]

with the normal Mach number, ${\operatorname{Ma}}_{n}=\frac{V_{n}}{c}$, and the sound speed $c=\sqrt{\gamma\frac{P}{\rho}}$.

As shown by Toro torobook , the solution for the rarefaction has a limiting vacuum solution for ${\operatorname{Ma}}_{n}\leq-2(\gamma-1)^{-1}.$ We will restrict our analysis to normal Mach numbers yielding strictly positive pressure solutions only (${\operatorname{Ma}}_{n}>\,-5$ for $\gamma=\frac{7}{5}$).

It is easy to see that using $P^{*}$ from 45, the entropy inequality 44 is still satisfied for $|V_{n}|\neq 0$, and the added entropy scales with the discrete value of $V_{n}$ at the boundary. Hence, for $h\rightarrow 0$, the discrete boundary condition converges to its physical counterpart, since $V_{n}\rightarrow 0$. The choice of $P^{*}$ from 45 appears to stabilize under-resolved simulations, which can be now explained by the fact that the boundary flux always adds entropy for $|V_{n}|\neq 0$.

Using approximate Riemann solvers for the boundary flux

A well known strategy in finite volume methods is to mirror only the velocity of the internal state and solve an approximate Riemann problem to get the boundary flux, mostly just because of a simpler implementation, since an approximate Riemann solver is already available and used for the fluxes between the elements. For DG methods, see also, for example, CHEN2017427 and PARSANI201588 where reflection conditions are proved to be entropy stable.

The mirror state is set so that the mass and energy flux are zero. Let the inner state be labeled $L$ and the outer $R$. then the inner and outer states that satisfy the mirror condition are

[TABLE]

We show below under what conditions on the normal velocity $V_{n}$ that the reflection condition is entropy stable for the Lax-Friedrichs, HLL and HLLC, Roe and EC-ES fluxes.

Lax-Friedrichs Flux

We start with the simplest approximate Riemann solver, the Lax-Friedrichs or Rusanov flux, which reads as

[TABLE]

Inserting the states from 46, we get

[TABLE]

The maximum wave speed is normally approximated from the largest leftgoing and rightgoing wave speed,

[TABLE]

and thus gives a definition of $P^{*}$

[TABLE]

which shows that the Lax-Friedrichs flux satisfies the entropy inequality 44.

HLL and HLLC Flux

The HLL flux torobook is written as

[TABLE]

The leftgoing and rightgoing wave speeds are $S^{L}=V_{n}^{L}-c^{L}=-V_{n}^{R}-c^{R}=-S^{R}$ and the HLL flux reduces to

[TABLE]

If we would choose $S^{R}$ to be the maximum wave speed, the HLL flux would reduce to the Lax-Friedrichs flux. However, with $S^{R}=V_{n}^{R}+c^{R}=-V_{n}+c$, an even simpler relation for $P^{*}$ is found, which also satisfies the entropy inequality

[TABLE]

For the HLLC flux torobook , one can show that since the Riemann problem is symmetric, the approximate wave speed of the contact discontinuity is $\lambda^{*}=0$ and, choosing $S^{R}=-V_{n}+c$, HLLC reduces to the HLL flux.

[TABLE]

Roe Flux

For the original Roe method without entropy fix torobook , the mean values are

[TABLE]

After some manipulations,

[TABLE]

with $\tilde{\lambda_{1}}=\tilde{V}_{n}-\tilde{c}=-\tilde{c}$, $\alpha_{1}=\rho V_{n}/\tilde{c}$ and $\tilde{\mathbf{K}}^{1}$ from torobook . This leads again to a definition of $P^{*}$

[TABLE]

which fulfills the entropy inequality as long as

[TABLE]

Thus, the Roe flux is entropy stable for shocks, but not for supersonic rarefactions.

EC-ES Fluxes

We can also apply an entropy conservative (EC) flux that is used for interior element interfaces and add an entropy stable dissipation term (ES) to compute the boundary flux via the mirrored states 46. This is exactly the strategy proposed in Parsani et al. PARSANI201588 to get the boundary flux. Such an EC-ES flux is presented in Winters et al. winters2017_ESmatrix

[TABLE]

where $\underline{D}$ is a dissipation matrix and the matrix \underline{H}\left\llbracket\mathbf{w}\right\rrbracket\simeq\left\llbracket\mathbf{u}\right\rrbracket is carefully derived from the left and right states. Details are given in winters2017_ESmatrix , where two approaches for the dissipation are distinguished. One is a Lax-Friedrichs-type dissipation, scaling with the maximum eigenvalue $\lambda_{\mathrm{max}}=|V_{n}|+c$ (referred to as ’EC-LF’). The other is a Roe-type dissipation computed via the eigenstructure of the matrix $(\underline{D}\,\underline{H})$ (referred to as ’EC-Roe’).

If we carefully insert the two mirrored boundary states into 59, we again get an equation for the modified pressure

[TABLE]

for the Lax-Friedrichs-type dissipation and

[TABLE]

for the Roe-type dissipation. Both approaches lead to an entropy stable boundary flux when using a mirrored state. Note that the modified pressure of the EC-Roe flux (61) exactly matches the one of the HLL flux 53.

4 Discussion

In the previous section we have shown conditions under which a specified wall flux is stable. In the linear analysis, the central numerical flux adds no dissipation and is neutrally stable. In the nonlinear analysis, entropy is not generated if the numerical wall pressure is equal to the internal pressure, $P^{*}=P^{\prime}$. For upwinded approximations, the amount of energy or entropy dissipation depends on the normal Mach number. Since the boundary condition is only imposed weakly through the numerical flux, the normal Mach number will not be exactly zero except in the convergence limit. In fact, flow computations (especially steady state ones) are usually initiated with an impulsive start, where the initial state is a uniform flow, and the normal Mach number is not zero. This has proved over time to be very robust in practice. The analysis above gives an explanation why.

In the linear analysis the dissipation due to imposing the boundary condition is proportional to the square of the normal Mach number. With an impulsive start initialization, this dissipation will be large. As the flow develops and the boundary condition is better enforced, the dissipation reduces, going away only as the approximate solution converges.

A similar effect is observed for the use of the different approximate Riemann solvers in the nonlinear analysis. In Fig. 1, we compare the entropy contribution

[TABLE]

for the different wall boundary fluxes, over a range of normal Mach numbers for $(\rho c)=1$ and $\gamma=7/5$. When the boundary condition is exactly fulfilled (${\operatorname{Ma}}_{n}=0$), the entropy contribution is zero. For low normal Mach numbers, all fluxes have the same behavior. Compared to the exact Riemann problem (RP), the Lax-Friedrichs flux and the EC-LF flux always produce more entropy whereas the the HLL flux produces less entropy for impinging velocities ${\operatorname{Ma}}_{n}>0$. The results of HLLC and EC-Roe fluxes are not plotted, as they coincide with the HLL flux. As shown in the analysis, the Roe flux produces a negative entropy change for supersonic rarefactions, implying that it is not suitable for all flow configurations.

Acknowledgements.

This work was supported by a grant from the Simons Foundation (#426393, David Kopriva). G.G. has been supported by the European Research Council (ERC) under the European Union’s Eights Framework Program Horizon 2020 with the research project Extreme, ERC grant agreement no. 714487. Florian Hindenlang thanks Eric Sonnendrücker and the Max-Planck Institute for Plasma Physics in Garching for their constant support.