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 [4, 7] has allowed for the analysis of these schemes and to connect them with penalty collocation and SBP finite difference schemes. For instance, in [5], 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. [3, 9, 10]. Only recently have they been studied for discontinuous Galerkin approximations. In [12], 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 [2], 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

$$\displaystyle \begin{aligned} {\mathbf u_t} + \sum_{i = 1}^3 {\frac{{\partial {\mathbf f_{i}}}}{{\partial {x_i}}}} = 0. \end{aligned} $$
(1)

The state vector contains the conservative variables

$$\displaystyle \begin{aligned} \mathbf u = \left[ { \varrho \;\; {\varrho \vec v} \;\; {E} } \right]^{T} = \left[ { \varrho \;\; {\varrho v_1} \;\; {\varrho v_2} \;\; {\varrho v_3} \;\; {E} } \right]^{T}.\end{aligned} $$
(2)

In standard form, the components of the advective fluxes are

$$\displaystyle \begin{aligned} \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], \end{aligned} $$
(3)

Here, \(\varrho ,\,\vec {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

$$\displaystyle \begin{aligned} p = (\gamma-1)\left(E - \frac{1}{2}\varrho\left\|\vec{v}\right\|{}^2\right), {} \end{aligned} $$
(4)

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

$$\displaystyle \begin{aligned} \overset{\,\leftrightarrow}{{\mathbf f}} = \left[ { {{\mathbf f_1}} \;\; {{\mathbf f_2}} \;\; {{\mathbf f_3}} } \right]^{T}\,, \end{aligned} $$
(5)

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

$$\displaystyle \begin{aligned} \mathbf u_{t} + \vec \nabla_x\cdot\overset{\,\leftrightarrow}{{\mathbf f}}=0. {} \end{aligned} $$
(6)

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 [11] for the symmetrization of the linearized equations, with the constants

$$\displaystyle \begin{aligned} a = \sqrt {\frac{{\gamma - 1}}{\gamma }} \bar c,\quad b = \frac{\bar c}{{\sqrt \gamma }}\,,\bar c=\sqrt{\frac{\gamma\bar{p}}{\bar{\varrho}}}, \end{aligned} $$
(7)

where \(\bar c\) is the sound speed of the constant mean state. The state variables become

$$\displaystyle \begin{aligned} \mathbf u = \left[\varrho'\;\; v_{1}\;\; v_{2}\;\; v_{3}\;\; p'\right]^{T}, \end{aligned} $$
(8)

where \(\vec {v}\) is the velocity perturbation from the mean state, and we introduce

$$\displaystyle \begin{aligned} \varrho' = b\frac{\tilde{\varrho}}{\bar{\varrho}}\,,\quad p'=\frac{1}{\bar{\varrho}a}\tilde{p} -\frac{1}{\sqrt{\gamma-1}}\varrho' , \end{aligned} $$
(9)

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

$$\displaystyle \begin{aligned} \mathbf f_{i}=\underline{A}_{i}\mathbf u,\quad \overset{\,\leftrightarrow}{{\mathbf f}} = \vec {\underline{A}}\mathbf u = \left(\underline{A}_{1}\hat x + \underline{A}_{2}\hat y+\underline{A}_{3}\hat z\right)\mathbf u, \end{aligned} $$
(10)

where [11]

$$\displaystyle \begin{aligned}{\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]\end{aligned} $$
(11)

are constant symmetric matrices.

The linear equations have the property that the L 2 norm of the solution over a domain Ω is bounded by terms of the boundary data on ∂ Ω, only. Let

$$\displaystyle \begin{aligned} \left\langle \mathbf v,\mathbf w\right\rangle = \int\limits_{\Omega}{\mathbf v^{T}\mathbf w\,dxdydz },\quad \left\langle \overset{\,\leftrightarrow}{{\mathbf f}},\overset{\,\leftrightarrow}{{\mathbf g}} \right\rangle = \int\limits_{\Omega} {\sum_{i = 1}^3 {\mathbf f_i^T{\mathbf g_i}} \,dxdydz }. \end{aligned} $$
(12)

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

$$\displaystyle \begin{aligned} \left\langle \vec\nabla_x\cdot\overset{\,\leftrightarrow}{{\mathbf f}},\mathbf u\right\rangle = \left\langle \vec\nabla_x\cdot\left(\vec{\underline{A}}\mathbf u\right),\mathbf u\right\rangle=\left\langle \nabla_x\mathbf u,\overset{\,\leftrightarrow}{{\mathbf f}}\right\rangle. \end{aligned} $$
(13)

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

$$\displaystyle \begin{aligned} \left\langle \vec\nabla_x\cdot\overset{\,\leftrightarrow}{{\mathbf f}},\mathbf u\right\rangle=\frac{1}{2}\int_{\partial\Omega} \mathbf u^{T}\overset{\,\leftrightarrow}{{\mathbf f}}\cdot\vec ndS, \end{aligned} $$
(14)

where \(\vec n\) is the outward normal to the surface of Ω. The norm of the solution therefore satisfies

$$\displaystyle \begin{aligned} \frac{d}{dt}|| \mathbf u||{}^{2}=-\int_{\partial\Omega} \mathbf u^{T}\overset{\,\leftrightarrow}{{\mathbf f}}\cdot\vec ndS. {} \end{aligned} $$
(15)

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

$$\displaystyle \begin{aligned} \mathbf u^{T}\overset{\,\leftrightarrow}{{\mathbf f}}\cdot\vec n=\mathbf u^{T}\left(\vec{\underline{A}}\cdot\vec n\right)\mathbf u = 2\left(\varrho b + ap\right)v_{n}+(\vec{\bar{v}}\cdot\vec{n})(\varrho^2 + |\vec{v}|{}^2+p^2), \end{aligned} $$
(16)

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

Therefore, with the no penetration wall condition v n = 0 applied,

$$\displaystyle \begin{aligned} \frac{d}{dt}|| \mathbf u||{}^{2} = 0, \end{aligned} $$
(17)

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 (γ − 1) for convenience) as

$$\displaystyle \begin{aligned} s(\mathbf u) = -\frac{\varrho \varsigma}{(\gamma-1)}\,, \end{aligned} $$
(18)

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

$$\displaystyle \begin{aligned} \vec f^{\,\varsigma}(\mathbf u) = \vec v \, s = -\frac{\varrho \varsigma \vec v}{(\gamma-1)}\,. \end{aligned} $$
(19)

Finally the entropy variables are

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

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

$$\displaystyle \begin{aligned} \mathbf w^T\,\mathbf u_t = \left(\frac{\partial s}{\partial\mathbf u}\right)^T\mathbf u_t = s_t(\mathbf u),\quad \mathbf w^T\, \vec{\nabla}_x\cdot\overset{\,\leftrightarrow}{{\mathbf f}} = \vec\nabla_x\cdot\vec{f}^{\,\varsigma}. {} \end{aligned} $$
(21)

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

$$\displaystyle \begin{aligned} \left\langle \mathbf w(\mathbf u), {\mathbf u_t}\right\rangle + \left\langle \mathbf{w}(\mathbf u),\vec\nabla_x\cdot\overset{\,\leftrightarrow}{{\mathbf f}}\right\rangle = 0\,. {} \end{aligned} $$
(22)

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

$$\displaystyle \begin{aligned} \left\langle s_t(\mathbf u ),1\right\rangle =- \left\langle \vec\nabla_x\cdot\vec{f}^{\,\varsigma},1\right\rangle = -\int\limits_{\partial\Omega}\left(\vec{f}^{\,\varsigma}\cdot\vec{n}\right) \,\operatorname{dS} \, {} \end{aligned} $$
(23)

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=\vec v \cdot \vec n=0\), the entropy is not changed by the slip-wall boundary condition, since

$$\displaystyle \begin{aligned} -\vec f^{\,\varsigma}\cdot \vec n=\frac{\varrho \varsigma}{(\gamma-1)} v_n=0. \end{aligned} $$
(24)

3 Stability Bounds for the DGSEM

The DGSEM is described in detail in [5] and elsewhere [1, 6]. We will only quickly summarize the approximation here. The domain, Ω 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, \(\vec a^{i}\), and transformation Jacobian, \(\mathcal J\). Equation (6) transform to another conservation law on the reference element as

$$\displaystyle \begin{aligned} \mathcal J {{\mathbf u}_t} + {\vec\nabla _\xi } \cdot {\overset{\,\leftrightarrow}{\tilde{\mathbf f}}} = 0, {} \end{aligned} $$
(25)

where \(\overset {\,\leftrightarrow }{\tilde {\mathbf f}}\) is the contravariant flux vector with components \(\tilde { \mathbf f}^{i} = \mathcal J \vec a^{i}\cdot \overset {\,\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 ρ. 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 [5] and [8] 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 [5]

$$\displaystyle \begin{aligned} {\mathbf U^{T}{\left\{{\tilde{\mathbf F}}^{*}-\frac{1}{2}\overset{\,\leftrightarrow}{\tilde{\mathbf F}}\cdot\hat n\right\}}}\ge 0, {} \end{aligned} $$
(26)

where \(\overset {\,\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 U is the approximation of the state vector. Since the contravariant fluxes are proportional to the normal fluxes [6], we can change the condition (26) to

$$\displaystyle \begin{aligned} B_{L}\equiv{\mathbf U^{T}{\left\{{\mathbf F}^{*}-\frac{1}{2}\overset{\,\leftrightarrow}{{\mathbf F}}\cdot\vec n\right\}}}\ge 0, {} \end{aligned} $$
(27)

For entropy stability of the nonlinear equations, the boundary stability condition shown in [5] is proportional to

$$\displaystyle \begin{aligned} B_{NL} \equiv \mathbf W^T\left(\mathbf F^{*}-\left(\overset{\,\leftrightarrow}{{\mathbf F}}\cdot\vec{n}\right)\right) + \left(\vec{F}^{\,\varsigma}\cdot\vec{n}\right)\ge 0, {} \end{aligned} $$
(28)

where \(\vec {F}^{\,\varsigma }\) is the polynomial interpolation of the entropy flux, \(\vec f^{{\,\varsigma }}\), and 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 '\;\; \vec V \;\; P']^T\) and the normal contravariant flux is proportional to

$$\displaystyle \begin{aligned} \overset{\,\leftrightarrow}{{\mathbf F}}\cdot\vec n = \vec{ {\underline{A}}}\cdot\vec n \,\mathbf U = \left[ bV_{n}\;\; n_1 Q\;\; n_2 Q\;\;n_3 Q\;\; aV_{n} \right]^{T}, {} \end{aligned} $$
(29)

where V n is the approximation of the normal velocity at the wall computed from the interior, Q = bρ′ + aP′, and (n 1, n 1, n 3) are the three components of the physical space normal vector, \(\vec n\). The numerical flux can be expressed as

$$\displaystyle \begin{aligned} \mathbf F^{*} = \vec{ {\underline{A}}}\cdot\vec n \,\mathbf U^{*}= \left[ bV^{*}_{n}\;\; n_1 Q^{*}\;\; n_2 Q^{*}\;\;n_3 Q^{*}\;\; aV^{*}_{n} \right]^{T}. {} \end{aligned} $$
(30)

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),

$$\displaystyle \begin{aligned} B_{L}=\frac{1}{2}\left\{ Q\left( 2V^{*}_{n}-V_{n}\right) + V_{n}\left( 2Q^{*}-Q\right)\right\} = \frac{1}{2}\left\{ 2QV^{*}_{n} + 2V_{n}\left( Q^{*}-Q\right)\right\} {} \end{aligned} $$
(31)

Substituting the wall boundary condition \(V^{*}_{n}=0\) yields the condition on Q ∗ for stability

$$\displaystyle \begin{aligned} V_{n}\left( Q^{*}-Q\right)\ge 0. \end{aligned} $$
(32)

Neutral stability is thus ensured if ρ ∗ and P ∗ are computed from the interior, i.e. ρ ∗ = ρ′, P ∗ = P′ 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 (ε = 1) normal Riemann flux and the central flux (ε = 0) for the linear system of equations is

$$\displaystyle \begin{aligned} \mathbf F^{*}\left( \mathbf U, \mathbf U^{\mathrm{ext}}\right) = \frac{1}{2}\left\{\overset{\,\leftrightarrow}{{\mathbf F}}\left(\mathbf U\right)\cdot\vec n + \overset{\,\leftrightarrow}{{\mathbf F}}\left(\mathbf U^{\mathrm{ext}}\right)\cdot\vec n\right\} - \frac{\varepsilon}{2}\left|{{\underline{A}}_{\,n}}\right|\left(\mathbf U^{\mathrm{ext}} - \mathbf U\right), \end{aligned} $$
(33)

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

$$\displaystyle \begin{aligned} \mathbf U^{\mathrm{ext}}=\left[ \rho' \;\; {} \;\; \left(\vec V - 2V_{n}\vec n\right) {} \;\; P' \right]^{T}. {} \end{aligned} $$
(34)

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

$$\displaystyle \begin{aligned} B_{L,0} &=\frac{1}{2}\mathbf U^{T} {{\underline{A}}_{\,n}}\mathbf U^{\mathrm{ext}} = \left[ {\begin{array}{*{20}{c}} \rho' &{}&\vec V&{}&P' \end{array}} \right]\left[ {\begin{array}{*{20}{c}} 0&{n_1 b}&{n_2 b}&{n_3 b}&0 \\ {n_1 b}&0&0&0&{n_1 a} \\ {n_2 b}&0&0&0&{n_2 a} \\ {n_3 b}&0&0&0&{n_3 a} \\ 0&{n_1 a}&{n_2 a}&{n_3 a}&0 \end{array}} \right]\left[ {\begin{array}{*{20}{c}} \rho' \\ {} \\ {\vec V - 2V_n\vec n} \\ {} \\ P' \end{array}} \right] \\ &= Q\left( { - V \cdot \vec n} \right) + \left(V \cdot \vec n\right)Q = 0, \end{aligned} $$
(35)

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 A̲n.

Substituting the exact upwind flux where ε = 1 into (27) and rearranging,

$$\displaystyle \begin{aligned} B_{L,1}= - \mathbf U^{T}\left| {{\underline{A}}_{\,n}^{-}}\right| \mathbf U^{\mathrm{ext}}+\frac{1}{2}\mathbf U^{T} \left| {{\underline{A}}_{\,n}}\right|\mathbf U, {} \end{aligned} $$
(36)

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

$$\displaystyle \begin{aligned} \mathbf U^{T} \left| {{\underline{A}}_{\,n}} \right|\mathbf U = \frac{1}{\bar c}Q^{2}+\bar c^{3}{\operatorname{Ma}}^{2}_{n}, \end{aligned} $$
(37)

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 U ext written in (34). Then

$$\displaystyle \begin{aligned} \mathbf U^{T}\left| {{\underline{A}}_{\,n}^{-}}\right| \mathbf U^{\mathrm{ext}} = \frac{1}{2\bar c}Q^{2} -\frac{\bar c^{3}}{2}{\operatorname{Ma}}^{2}. \end{aligned} $$
(38)

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

$$\displaystyle \begin{aligned} B_{L,1}=\bar c^{3}\mathrm{Ma}^{2}_{n}\ge 0, \end{aligned} $$
(39)

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

$$\displaystyle \begin{aligned} B_{L,LF} =-\frac{\left|\lambda\right|{}_{\mathrm{max}}}{2}\mathbf U^{T}\left(\mathbf U^{\mathrm{ext}}-\mathbf U\right)=\bar c^{2}{\left|\lambda\right|{}_{\mathrm{max}}}{\operatorname{Ma}}_{n}^{2}\ge 0 \end{aligned} $$
(40)

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 ρ ∗ = ρ′ and

$$\displaystyle \begin{aligned} P^{*} = P' + \frac{\bar c^{3}}{a}{\left|\lambda\right|{}_{\mathrm{max}}}{\operatorname{Ma}}_{n}. \end{aligned} $$
(41)

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

$$\displaystyle \begin{aligned} V_{n}\left( Q^{*}-Q\right) = \bar c^{2}{\left|\lambda\right|{}_{\mathrm{max}}}{\operatorname{Ma}}_{n}^{2}\ge 0. \end{aligned} $$
(42)

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 ,\vec {V},P)\) describe the full nonlinear state.

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

$$\displaystyle \begin{aligned} \left(\overset{\,\leftrightarrow}{{\mathbf F}} \cdot \vec n\right)^* = \left[ {\begin{array}{*{20}c} 0 \\ P^*\,\vec n \\ 0 \\ \end{array} } \right] \end{aligned} $$
(43)

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

$$\displaystyle \begin{aligned} \begin{aligned} -\rho {V_n}\left(\frac{\gamma-\varsigma}{(\gamma-1)}-\beta|| {\vec V}||{}^2\right)&\\ + 2\beta V_n \left(P^*- {P}-\rho || {\vec V}||{}^2\right) + 2\beta V_n \left(\rho E+P\right) - \frac{\rho \varsigma V_n}{(\gamma-1)}=&\\ -\rho V_n\left(\frac{\gamma}{(\gamma-1)} -\beta|| {\vec V}||{}^2\right) + 2\beta V_n\left(P^*+\rho E-\rho|| {\vec V}||{}^2\right)=&\\ \frac{\rho V_n}{P}\left(-\frac{\gamma}{(\gamma-1)}P +\frac{1}{2}\rho|| {\vec V}||{}^2 + P^*+\frac{P}{(\gamma-1)}-\frac{1}{2}\rho|| {\vec V}||{}^2\right) =&\\ \rho V_n\left(\frac{P^*}{P}-1\right)&\geq 0 {} \end{aligned} \end{aligned} $$
(44)

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.

3.2.1 Exact Solution of the Riemann Problem

In [14] a symmetric 1D Riemann problem is exactly solved following Toro [13], 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

$$\displaystyle \begin{aligned} \left(\frac{P^*}{P}\right)_{\text{RP}}=\left\{\begin{array}{*{20}lll} 1+\gamma{\operatorname{Ma}}_n\left(\frac{(\gamma+1)}{4}{\operatorname{Ma}}_n+\sqrt{\left(\frac{(\gamma+1)}{4}{\operatorname{Ma}}_n\right)^2+1}\right)&> 1 &\quad \text{for}\quad V_n>0 \\ \left(1+\frac{1}{2}(\gamma-1){\operatorname{Ma}}_n\right)^{\frac{2\gamma}{(\gamma-1)}}&\leq 1 &\quad \text{for}\quad V_n\leq 0 \end{array}\right.{} \end{aligned} $$
(45)

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 [13], 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|≠ 0, and the added entropy scales with the discrete value of V n at the boundary. Hence, for h → 0, the discrete boundary condition converges to its physical counterpart, since V n → 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|≠ 0.

3.2.2 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, [2] and [12] 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

$$\displaystyle \begin{aligned} \mathbf U^L = \left[\rho\;\;\rho\vec V_n\;\; E\right]^T\,,\quad \mathbf U^R=\left[\rho\;\;\rho(\vec V - 2V_n\vec n )\;\; E\right]^T \end{aligned} $$
(46)

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.

3.2.2.1 Lax-Friedrichs Flux

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

$$\displaystyle \begin{aligned} \left(\overset{\,\leftrightarrow}{{\mathbf F}} \cdot \vec n\right)_{\text{LF}}^*=\frac{1}{2}\vec n\cdot\left(\overset{\,\leftrightarrow}{{\mathbf F}} (\mathbf U^L)+\overset{\,\leftrightarrow}{{\mathbf F}} (\mathbf U^R)\right)-\frac{\left|\lambda\right|{}_{\mathrm{max}}}{2} (\mathbf U^R-\mathbf U^L). \end{aligned} $$
(47)

Inserting the states from (46), we get

$$\displaystyle \begin{aligned} \left(\overset{\,\leftrightarrow}{{\mathbf F}} \cdot \vec n\right)_{\text{LF}}^* = \left[ {\begin{array}{*{20}c} {0} \\ {(\rho V_n^2 + P)\,\vec n} \\ {0} \\ \end{array} } \right]-\frac{\lambda_{\mathrm{max}}}{2} \left[ {\begin{array}{*{20}c} {0} \\ {-2\rho V_n\vec n} \\ {0} \\ \end{array} } \right] = \left[ {\begin{array}{*{20}c} {0} \\ {(\rho V_n^2+\rho V_n\lambda_{\mathrm{max}} + P)\,\vec n} \\ {0} \\ \end{array} } \right]. \end{aligned} $$
(48)

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

$$\displaystyle \begin{aligned} \lambda_{\mathrm{max}}&=\max(|V_n^L|+c^L,|V_n^R|+c^R)=|V_n|+c \,,\quad \text{since}\quad c^L=c^R=c\,, \\ &\quad V_n=V_n^L=-V_n^R \end{aligned} $$
(49)

and thus gives a definition of P ∗

$$\displaystyle \begin{aligned} \left(\frac{P^*}{P}\right)_{\text{LF}} &= 1+\gamma{\operatorname{Ma}}_n\left({\operatorname{Ma}}_n+|{\operatorname{Ma}}_n|+1\right) \\ &=\left\{\begin{array}{*{20}lll} 1+\gamma{\operatorname{Ma}}_n(2{\operatorname{Ma}}_n+1)&> 1 &\quad \text{for}\quad V_n>0 \\ 1+\gamma{\operatorname{Ma}}_n &\leq 1 &\quad \text{for}\quad V_n\leq 0 \\ \end{array}\right., {} \end{aligned} $$
(50)

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

3.2.2.2 HLL and HLLC Flux

The HLL flux [13] is written as

$$\displaystyle \begin{aligned} \left(\overset{\,\leftrightarrow}{{\mathbf F}} \cdot \vec n\right)_{\text{HLL}}^*=\frac{1}{S^R-S^L}\left(\vec n\cdot\left(S^R\overset{\,\leftrightarrow}{{\mathbf F}} (\mathbf U^L)-S^L\overset{\,\leftrightarrow}{{\mathbf F}} (\mathbf U^R)\right)+S^L S^R\left(\mathbf U^R-\mathbf U^L\right)\right). \end{aligned} $$
(51)

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

$$\displaystyle \begin{aligned} \left(\overset{\,\leftrightarrow}{{\mathbf F}} \cdot \vec n\right)_{\text{HLL}}^*=\frac{1}{2} \vec n\cdot\left(\overset{\,\leftrightarrow}{{\mathbf F}} (\mathbf U^L)+\overset{\,\leftrightarrow}{{\mathbf F}} (\mathbf U^R)\right) -\frac{S^R}{2}\left(\mathbf U^R-\mathbf U^L\right). \end{aligned} $$
(52)

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

$$\displaystyle \begin{aligned} \left(\frac{P^*}{P}\right)_{\text{HLL}} = 1+\gamma{\operatorname{Ma}}_n\left\{\begin{array}{*{20}ll} > 1 &\quad \text{for}\quad V_n>0 \\ \leq 1 &\quad \text{for}\quad V_n\leq 0 \\ \end{array}\right.{} \end{aligned} $$
(53)

For the HLLC flux [13], one can show that since the Riemann problem is symmetric, the approximate wave speed of the contact discontinuity is λ ∗ = 0 and, choosing S R = −V n + c, HLLC reduces to the HLL flux.

$$\displaystyle \begin{aligned} \left(\frac{P^*}{P}\right)_{\text{HLLC}} =1+\gamma{\operatorname{Ma}}_n=\left(\frac{P^*}{P}\right)_{\text{HLL}} {} \end{aligned} $$
(54)
3.2.2.3 Roe Flux

For the original Roe method without entropy fix [13], the mean values are

$$\displaystyle \begin{aligned} \tilde{V}_n &= \frac{\sqrt{\rho^L}V_n^L +\sqrt{\rho^R}V_n^R}{\sqrt{\rho^L} +\sqrt{\rho^R}} =0\,,\quad \tilde{V}_{t_1}=V_{t_1}\,,\quad \tilde{V}_{t_2}=V_{t_2}\,, \\ \tilde{c} &= c\sqrt{1+\frac{(\gamma-1)}{2}{\operatorname{Ma}}_n^2} . \end{aligned} $$
(55)

After some manipulations,

$$\displaystyle \begin{aligned} \left(\overset{\,\leftrightarrow}{{\mathbf F}} \cdot \vec n\right)_{\text{Roe}}^* &= \left(\overset{\,\leftrightarrow}{{\mathbf F}} \cdot \vec n\right) + \tilde{\lambda}_1\tilde{\alpha}_1\tilde{\mathbf K}^{1}= \left(\overset{\,\leftrightarrow}{{\mathbf F}}\cdot \vec n\right)+(-\tilde{c})\frac{\rho V_n}{\tilde{c}}\left[ {\begin{array}{*{20}c} {1} \\ {-\tilde{c}} \\ V_{t_1}\\V_{t_2} \\ { \frac{1}{\rho}(\rho E + P)} \\ \end{array} } \right] \\ &=\left[ {\begin{array}{*{20}c} {0} \\ {(\rho V_n^2 +\rho V_n \tilde{c} + P)\,\vec n} \\ {0} \\ \end{array} } \right]\,. \end{aligned} $$
(56)

with \(\tilde {\lambda _1}=\tilde {V}_n-\tilde {c}=-\tilde {c}\), \(\alpha _1=\rho V_n /\tilde {c}\) and \(\tilde {\mathbf K}^{1}\) from [13]. This leads again to a definition of P ∗

$$\displaystyle \begin{aligned} \left(\frac{P^*}{P}\right)_{\text{Roe}} =1+\gamma {\operatorname{Ma}}_n\left({\operatorname{Ma}}_n+\sqrt{1+\frac{(\gamma-1)}{2}{\operatorname{Ma}}_n^2}\right), {} \end{aligned} $$
(57)

which fulfills the entropy inequality as long as

$$\displaystyle \begin{aligned} {\operatorname{Ma}}_n \geq -\sqrt{\frac{2}{3-\gamma}} \,,\quad \text{for}\,\gamma=\frac{7}{5}\quad {\operatorname{Ma}}_n > -1.12 \,. \end{aligned} $$
(58)

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

3.2.2.4 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. [12] to get the boundary flux. Such an EC-ES flux is presented in Winters et al. [15]

$$\displaystyle \begin{aligned} \left(\overset{\,\leftrightarrow}{{\mathbf F}} \cdot \vec n\right)_{\text{ES}}^*=\overset{\,\leftrightarrow}{{\mathbf F}}_{\text{EC}}\left(\mathbf U_L,\mathbf U_R\right) \cdot \vec n - \frac{1}{2} \underline{D}\,\underline{H} \left(\mathbf W^R-\mathbf W^L\right) {} \end{aligned} $$
(59)

where D̲ is a dissipation matrix and the matrix \(\underline {H} \left [\![ \mathbf w\right ]\!] \simeq \left [\![ \mathbf {u}\right ]\!] \) is carefully derived from the left and right states. Details are given in [15], where two approaches for the dissipation are distinguished. One is a Lax-Friedrichs-type dissipation, scaling with the maximum eigenvalue λ max = |V n| + c (referred to as ‘EC-LF’). The other is a Roe-type dissipation computed via the eigenstructure of the matrix (D̲ 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

$$\displaystyle \begin{aligned} \left(\frac{P^*}{P}\right)_{\text{EC-LF}} =1+\gamma {\operatorname{Ma}}_n\left(|{\operatorname{Ma}}_n|+1\right) {} \end{aligned} $$
(60)

for the Lax-Friedrichs-type dissipation and

$$\displaystyle \begin{aligned} \left(\frac{P^*}{P}\right)_{\text{EC-Roe}} =1+\gamma {\operatorname{Ma}}_n {} \end{aligned} $$
(61)

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′. 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

$$\displaystyle \begin{aligned} \Delta s = (\rho c) {\operatorname{Ma}}_n\left(\frac{P^*}{P}-1\right) {} \end{aligned} $$
(62)

for the different wall boundary fluxes, over a range of normal Mach numbers for (ρc) = 1 and γ = 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 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.

Fig. 1
figure 1

Entropy contribution Δs (62) produced by the wall boundary flux. RP refers to the exact Riemann problem (45), LF to (50), EC-LF to (60), HLL to (53) and Roe to (57). Plotted over the normal Mach number ranges \(|{\operatorname {Ma}}_n|\leq 5\) on the top and restricted to \(|{\operatorname {Ma}}_n|\leq 1\) on the bottom

Full size image