Stability of Wall Boundary Condition Procedures for Discontinuous Galerkin Spectral Element Approximations of the Compressible Euler Equations
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.
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
The state vector contains the conservative variables
In standard form, the components of the advective fluxes are
Here, 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
where denotes the adiabatic coefficient. For a compact notation that simplifies the analysis, we define block vectors (with the double arrow)
so that the system of equations can be written in the compact form
The linear Euler equations are derived by linearizing about a constant mean state . We follow Nordstrom:2005qy for the symmetrization of the linearized equations, with the constants
where is the sound speed of the constant mean state. The state variables become
where is the velocity perturbation from the mean state, and we introduce
which depend on the density and pressure perturbations . The flux vectors are
are constant symmetric matrices.
The linear equations have the property that the norm of the solution over a domain is bounded by terms of the boundary data on , only. Let
represent the inner product of two state vectors and and two block vectors and , respectively. Since the coefficient matrices are constant the product rule and symmetry of implies
Then it follows from Gauss’ law (integration by parts) that
where is the outward normal to the surface of . The norm of the solution therefore satisfies
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
where is the wall normal velocity, . Note that here, the mean flow must be chosen such that the normal flow vanishes at the wall boundary , so that the boundary condition makes physical sense.
Therefore, with the no penetration wall condition applied,
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 for convenience) as
where 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
Finally the entropy variables are
The entropy pair contracts the solution and fluxes, meaning that it satisfies the relations
When we multiply (6) with the entropy variables and integrate over the domain,
Next we use the properties of the entropy pair to contract creftype 22 and use integration by parts to get
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 , the entropy is not changed by the slip-wall boundary condition, since
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, is subdivided into non-overlapping, conforming, hexahedral elements. Each element is mapped to the reference element . Associated with the transformation from the reference element is a set of contravariant coordinate vectors, , and transformation Jacobian, . The equations (6) transform to another conservation law on the reference element as
where is the contravariant flux vector with components .
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 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 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 satisfies the bound Gassner2018
where is the polynomial interpolation of the contravariant flux from the interior, is the reference space outward normal direction, and 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
For entropy stability of the nonlinear equations, the boundary stability condition shown in Gassner2018 is proportional to
where is the polynomial interpolation of the entropy flux, , and 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 , 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 and the normal contravariant flux is proportional to
where is the approximation of the normal velocity at the wall computed from the interior, , and are the three components of the physical space normal vector, . The numerical flux can be expressed as
Substituting the wall boundary condition yields the condition on for stability
Neutral stability is thus ensured if and are computed from the interior, i.e. , so that .
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 () normal Riemann flux and the central flux () for the linear system of equations is
where 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,
For , using the central (averaged) numerical flux, the interior flux contribution cancels and condition (27) reduces to
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 .
Substituting the exact upwind flux where into (27) and rearranging,
where 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
where 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 written in (34). Then
Therefore, using the upwind numerical flux, (36) becomes
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, . The jump term is added to the central (averaged) flux so
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
A similar, though more complicated, modified 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 describe the full nonlinear state.
For the nonlinear equations, we construct the numerical flux for a slip-wall as
where we imposed leading to a flux with no mass or energy transfer, and we introduce a wall pressure , whose value will be chosen to ensure consistency and stability.
After some manipulations, the discrete entropy stability condition (28) becomes
Therefore if we choose , to be the internal pressure, the boundary flux does not contribute to the total entropy, independent of the inner normal velocity . A value of 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 , accounting for the fact that 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
with the normal Mach number, , and the sound speed .
As shown by Toro torobook , the solution for the rarefaction has a limiting vacuum solution for We will restrict our analysis to normal Mach numbers yielding strictly positive pressure solutions only ( for ).
It is easy to see that using from creftype 45, the entropy inequality creftype 44 is still satisfied for , and the added entropy scales with the discrete value of at the boundary. Hence, for , the discrete boundary condition converges to its physical counterpart, since . The choice of from creftype 45 appears to stabilize under-resolved simulations, which can be now explained by the fact that the boundary flux always adds entropy for .
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 and the outer . then the inner and outer states that satisfy the mirror condition are
We show below under what conditions on the normal velocity that the reflection condition is entropy stable for the Lax-Friedrichs, HLL and HLLC, Roe and EC-ES fluxes.
We start with the simplest approximate Riemann solver, the Lax-Friedrichs or Rusanov flux, which reads as
Inserting the states from creftype 46, we get
The maximum wave speed is normally approximated from the largest leftgoing and rightgoing wave speed,
and thus gives a definition of
which shows that the Lax-Friedrichs flux satisfies the entropy inequality creftype 44.
HLL and HLLC Flux
The HLL flux torobook is written as
The leftgoing and rightgoing wave speeds are and the HLL flux reduces to
If we would choose to be the maximum wave speed, the HLL flux would reduce to the Lax-Friedrichs flux. However, with , an even simpler relation for is found, which also satisfies the entropy inequality
For the HLLC flux torobook , one can show that since the Riemann problem is symmetric, the approximate wave speed of the contact discontinuity is and, choosing , HLLC reduces to the HLL flux.
For the original Roe method without entropy fix torobook , the mean values are
After some manipulations,
with , and from torobook . This leads again to a definition of
which fulfills the entropy inequality as long as
Thus, the Roe flux is entropy stable for shocks, but not for supersonic rarefactions.
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 creftype 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
where is a dissipation matrix and the matrix 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 (referred to as ’EC-LF’). The other is a Roe-type dissipation computed via the eigenstructure of the matrix (referred to as ’EC-Roe’).
If we carefully insert the two mirrored boundary states into creftype 59, we again get an equation for the modified pressure
for the Lax-Friedrichs-type dissipation and
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 creftype 53.
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, . 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
for the different wall boundary fluxes, over a range of normal Mach numbers for and . When the boundary condition is exactly fulfilled (), 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 . 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.
-  K. Black. A conservative spectral element method for the approximation of compressible fluid flow. KYBERNETIKA, 35(1):133–146, 1999.
-  Tianheng Chen and Chi-Wang Shu. Entropy stable high order discontinuous Galerkin methods with suitable quadrature rules for hyperbolic conservation laws. Journal of Computational Physics, 345:427 – 461, 2017.
-  T. Fisher, M.H. Carpenter, J. Nordström, N. K. Yamaleev, and C. Swanson. Discretely conservative finite-difference formulations for nonlinear conservation laws in split form: Theory and boundary conditions. J. Comput. Phys., 234:pp. 353–375, 2013.
-  Gregor J. Gassner, Andrew R. Winters, Florian J. Hindenlang, and David A. Kopriva. The BR1 scheme is stable for the compressible Navier–Stokes equations. Journal of Scientific Computing, Apr 2018.
-  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, 2016.
-  D. A. Kopriva and G. Gassner. On the quadrature and weak form choices in collocation type discontinuous Galerkin spectral element methods. Journal of Scientific Computing, 44(2):136–155, 2010.
-  D. A. Kopriva and G. J. Gassner. An energy stable discontinuous Galerkin spectral element discretization for variable coefficient advection problems. SIAM J. Sci. Comp., 36(4):A2076–A2099, 2014.
-  David A. Kopriva. Implementing Spectral Methods for Partial Differential Equations. Scientific Computation. Springer, May 2009.
-  Jan Nordström. Conservative finite difference formulations, variable coefficients, energy estimates and artificial dissipation. Journal of Scientific Computing, 29(3):375–404, 2006.
-  Jan Nordström. A roadmap to well posed and stable problems in computational physics. Journal Of Scientific Computing, DOI 10.1007/s10915-016-0303-9, 2016.
-  Jan Nordström and Magnus Svard. Well-posed boundary conditions for the Navier–Stokes equations. SIAM J. on Numerical Analysis, 43(3):1231–1255, 2005.
-  Matteo Parsani, Mark H. Carpenter, and Eric J. Nielsen. Entropy stable wall boundary conditions for the three-dimensional compressible Navier–Stokes equations. Journal of Computational Physics, 292:88 – 113, 2015.
-  E.F. Toro. Riemann Solvers and Numerical Methods for Fluid Dynamics. Springer, June 2009.
-  Jacobus J.W. van der Vegt and H. van der Ven. Slip flow boundary conditions in discontinuous Galerkin discretizations of the Euler equations of gas dynamics. In H.A. Mang and F.G. Rammenstorfer, editors, Proceedings of the 5th World Congress on Computational Mechanics (WCCM V), number NLR-TP in Technical Publications, pages 1–16. National Aerospace Laboratory, NLR, 7 2002.
-  Andrew R. Winters, Dominik Derigs, Gregor J. Gassner, and Stefanie Walch. A uniquely defined entropy stable matrix dissipation operator for high Mach number ideal MHD and compressible Euler simulations. Journal of Computational Physics, 332:274 – 289, 2017.