Sub-Cell Shock Capturing for Discontinuous Galerkin Methods

TLDR: A shock capturing strategy for higher order Discontinuous Galerkin approximations of scalar conservation laws is presented and it is shown that the proposed approach is capable of capturing the shock as a sharp, but smooth profile, which is typically contained within one element.

Abstract: A shock capturing strategy for higher order Discontinuous Galerkin approximations of scalar conservation laws is presented. We show how the original explicit artificial viscosity methods proposed over fifty years ago for finite volume methods, can be used very eectively in the context of high order approximations. Rather than relying on the dissipation inherent in Discontinuous Galerkin approximations, we add an artificial viscosity term which is aimed at eliminating the high frequencies in the solution, thus eliminating Gibbs-type oscillations. We note that the amount of viscosity required for stability is determined by the resolution of the approximating space and therefore decreases with the order of the approximating polynomial. Unlike classical finite volume artificial viscosity methods, where the shock is spread over several computational cells, we show that the proposed approach is capable of capturing the shock as a sharp, but smooth profile, which is typically contained within one element. The method is complemented with a shock detection algorithm which is based on the rate of decay of the expansion coecients of the solution when this is expressed in a hierarchical orthonormal basis. For the Euler equations, we consider and discuss the performance of several forms of the artificial viscosity term.

Figures

Figure 4. Supersonic flow around a NACA0012 airfoil at Mach 1.5, with interpolation of degree p = 4. Again we see how the shocks are resolved within the elements, and the diffusion is only applied to elements around the shocks.

Figure 1. Solution of the inviscid Burgers’ equation. The interpolation order is p = 10, and the vertical lines indicate the element boundaries. Note how the shock is well-resolved within one element, and that the jumps between the elements are small.

Figure 3. Transonic flow around a NACA0012 airfoil at Mach 0.8, with interpolation of degree p = 4. Note how the shock is accurately resolved within the elements. The high order scheme provides low dissipation away from the shock, as can be seen in the entropy plot (middle). The bottom plot shows the applied diffusion, with gray color representing no diffusion added.

Figure 5. Mach 2 flow around a cylinder, with interpolation of degree p = 5. Note that the coarse grid shown is the actual grid used in the computation. Some of the oscillation observed are due to the coarseness of the grid. We observe that the shocks are resolved within the elements.

Figure 6. Mach 5 flow around a cylinder, with interpolation of degree p = 5. Again, we observe that the shocks are resolved within the elements.

Table 1. Normalized shock thickness for density, Mach number, entropy and enthalpy using the three different viscosity models as a function of the upstream Mach number.

Frequently asked questions

Q1. What have the authors contributed in "Sub-cell shock capturing for discontinuous galerkin methods" ?

A shock capturing strategy for higher order Discontinuous Galerkin approximations of scalar conservation laws is presented. The authors show how the original explicit artificial viscosity methods proposed over fifty years ago for finite volume methods, can be used very effectively in the context of high order approximations. Rather than relying on the dissipation inherent in Discontinuous Galerkin approximations, the authors add an artificial viscosity term which is aimed at eliminating the high frequencies in the solution, thus eliminating Gibbs-type oscillations. Unlike classical finite volume artificial viscosity methods, where the shock is spread over several computational cells, the authors show that the proposed approach is capable of capturing the shock as a sharp, but smooth profile, which is typically contained within one element. For the Euler equations, the authors consider and discuss the performance of several forms of the artificial viscosity term. 

Q2. What is the main drawback of the DG algorithm?

Reducing the order of the interpolating polynomial increases the inter-element jumps and hence the amount of dissipation naturally added by the DG algorithm. 

Q3. How can the authors obtain stable discretizations with optimal apriori error estimates?

By appropriately choosing the interface fluxes,13 it is possible to obtain stable discretizations with optimal apriori error estimates using equal order interpolations for both u and q. 

Q4. Why is the entropy boundary layer thicker than the other two models?

the authors note that because of the large overshoot in entropy, the entropy boundary layer is considerably thicker than with the other two models. 

Q5. What is the obvious solution to the problem?

An obvious solution, almost universally accepted, is that near the shocks, the solution can only be first order accurate and therefore, if one adaptively refines that region by locally reducing h, one may be able to alleviate the problem. 

Q6. What is the entropy coefficient for the three models?

For the physical models a single constant has been determined and then, the Mach number scaling described above has been employed to determine the value of the viscosity coefficient for each flow condition. 

Q7. What is the effect of the LDG approach on the viscosity coefficient?

The authors have found out however that, for high Mach numbers, the value of the viscosity coefficient needs to be increased to maintain stability and this results in wider shocks. 

Q8. What is the rationale behind introducing viscosity to the original equations?

Introducing viscosity to the original equations requires discretizing second order derivatives which is not straightforward when the approximating space is discontinuous. 

Q9. What is the simplest way to express the solution of order p within each element?

The authors express the solution of order p within each element in terms of an orthogonal basis asu = N(p)∑ i=1 uiψi , (5)where N(p) is the total number of terms in the expansion and ψi are the basis functions. 

Q10. How can the authors determine the value of u for a normal shock?

The value of ∆u for a normal shock can be easily determined as a function of the upstream Mach number, M∞, and the upstream velocity u∞, as ∆u = 2u∞(M2∞−1)/[(γ+1)M2∞]. 

Q11. What is the amount of viscosity required when the shock crosses between elements?

when the shock crosses between elements (as for T = 0.50, right plots) the amount of required viscosity is typically small. 

Q12. Why can't the discretization of the model term be carried out with the standard DG?

The authors note that the discretization of the model term in equation (2) can not be carried out with the standard DG method due to the presence of the higher derivatives. 

Q13. What is the value of u for a normal shock?

the authors have thatδs∆u ρ∗µ∗ ≈ 1 . (19)Here, ∆u is the change in normal velocity across the shock and ρ∗ and µ∗ are the values of the density and viscosity evaluated at sonic conditions. 

Q14. How do the authors determine a suitable sensor for discontinuities?

In order to determine a suitable sensor for discontinuities, the authors write the solution within each element in terms of a hierarchical family of orthogonal polynomials. 

Q15. What is the definition of a mesh adaptive strategy?

As a consequence, mesh adaptive strategies must incorporate some degree of directionality if they are to lead to an effective approach, specially in three dimensions.