0%" bgcolor=white>
Shock Capturing with Higher-Order, State Based Artificial Viscosity
Shock Capturing with Higher-Order, PDE-Based
Articial Viscosity
Garrett E. Barter and David L. Darmofal Massachusetts Institute of Technology, Cambridge, MA 02139
Articial viscosity can be combined with a higher-order discontinuous Galerkin (DG)
discretization to resolve a shock layer within a single cell. However, when a piecewise-
constant articial viscosity model is employed with an otherwise higher-order approxima-
tion, element-to-element variations in the articial viscosity arising at the shock induce
oscillations in state gradients and pollute the downstream ow. To alleviate these dif-
culties, this work proposes a new higher-order, state-based articial viscosity with an
associated governing PDE. In the governing PDE, the shock sensor acts as a forcing term,
driving the articial viscosity to a non-zero value where it is necessary. The decay rate
of the higher-order solution modes and edge-based jumps are both shown to be reliable
shock indicators. This new approach leads to a smooth, higher-order representation of the
articial viscosity, that evolves in time with the solution. Additionally, an articial dissipa-
tion operator that preserves total enthalpy is introduced. The combination of higher-order,
PDE-based artical viscosity and enthalpy-preserving dissipation operator is shown to over-
come the disadvantages of the piecewise-constant articial viscosity, while achieving greater
robustness on ows with strong shocks.
I.
Introduction
In the past decade, the discontinuous Galerkin (DG) nite element method (FEM) has matured to become
a viable alternative to nite volume schemes. Whereas higher-order accuracy in nite volume schemes is
typically achieved by polynomial reconstruction of neighboring cell averages, in the DG context, higher-
order approximations are realized by increasing the order of the approximating polynomial, p, within each
element. This serves to maintain a nearest neighbor numerical stencil for all solution orders at the cost of
additional degrees of freedom on a given mesh. This compactness makes DG well suited for unstructured
grids, h and p adaptivity, and parallelization. Much of the ground work for DG methods was laid down by
Karniadakis, Cockburn and Shu.
111
Bassi & Rebay and Bey & Oden demonstrated the capabilities of DG
for the Euler and Navier-Stokes equations (including RANS).
1216
Persson and Peraire extended this work
with sub-cell shock capturing and improved methods for discretizing turbulence models.
17, 18
Recent work
has also focused on improving DG solution methods.
1923
For more information on DG methods, consult
the review by Cockburn & Shu.
5
The focus of this paper is shock capturing for higher-order (p 1) DG approximations. For rst-
order, piecewise-constant solutions in DG, the inter-element jumps, combined with an appropriate upwind
ux function, introduce enough numerical dissipation to obviate the need for shock capturing. However,
for higher-order approximations, spurious numerical oscillations manifest themselves in the proximity of
discontinuities and additional shock capturing steps must be taken. Many researchers have borrowed from
the shock capturing experience of nite dierence and nite volume schemes to develop shock capturing
approaches for DG FEM.
One classical approach to shock capturing is the addition of articial viscosity, pioneered by von Neumann
and Richtmyer.
24
Articial viscosity expands the thickness of the shock layer so that it safely exceeds the
resolution length scales of the numerical method and eliminates the spurious oscillations. Persson and
Peraire
17
introduced a p-dependent articial viscosity and demonstrated that higher-order representations Research Assistant, Student Member AIAA. Email gbarter@mit.edu Associate Professor, Senior Member AIAA. Email darmofal@mit.edu
1 of 14
American Institute of Aeronautics and Astronautics
18th AIAA Computational Fluid Dynamics Conference
25 - 28 June 2007, Miami, FL
AIAA 2007-3823
Copyright � 2007 by the American Institute of Aeronautics and Astronautics, Inc. All rights reserved. and a piecewise-constant articial viscosity can be combined to produce sub-cell shock resolution. Specically,
by introducing an articial viscosity which scales with h/p, the shock width
s
= Ch/p, where C is usually
around 24. Thus, for suciently high p,
s
= Ch/p < h and the shock can be captured within a single
element. To locate the shocks in the oweld, Persson and Peraire developed a sensor based on the magnitude
of the highest-order coecients in an orthonormal representation of the solution.
This research builds upon the benets of articial viscosity for shock capturing in DG. As will be described
below, a piecewise-constant articial viscosity has some inherent shortcomings. Specically, element-to-
element variations can lead to oscillations in state gradients and disparate equilibrium shock-jump conditions
in neighboring elements. This can potentially corrupt the smoothness and accuracy of the downstream
oweld. We seek to develop a smoother representation of articial viscosity, without sacricing the compact
numerical stencil that makes DG an attractive scheme.
In Section II, we motivate a smooth representation of articial viscosity by highlighting the shortcomings
of the piecewise-constant approach using a simplied test problem. In Section III, we present a PDE-based
articial viscosity to obtain the necessary smoothness, and two dierent shock sensors. This is followed by
the discretization of our new articial viscosity model in Section IV and a discussion of the extra steps taken
to preserve total enthalpy through the shock. Numerical results using the PDE-based articial viscosity are
presented in Section V.
II.
Articial Viscosity Comparison
The one-dimensional viscous Burgers equation is employed to demonstrate the benets of a smooth vari-
ation in articial viscosity, compared to a piecewise-constant approach. The governing equation is modied
to support a steady-state shock solution with a forcing term, x
1
2 u
2
= u +
x
(x) u
x
+ f (x),
for x R
(1)
where u(x, t) is the state variable, (x) is the viscosity, is a constant and the forcing term, f (x), is set
such that the exact, steady-state solution has a shock at x = 0.
The viscosity, (x), is prescribed to be either a piecewise-constant or smooth Gaussian function, as
depicted in Figure 1. The piecewise-constant viscosity is applied to the cells immediately adjacent to the
shock location with adjustable amplitude. The Gaussian distribution of viscosity is specied to have a
standard deviation equal to the cell size and the same total area as the piecewise-constant rectangle between
x [h, h].
Figure 1. Distributions of piecewise-constant and Gaussian articial viscosity as
applied to the 1-D modied Burgers equation.
To perform the comparison, Equation (1) is discretized using sixth order Legendre polynomials in DG
FEM with the second method of Bassi and Rebay (BR2) for the elliptic terms. An H
1
norm of the error
2 of 14
American Institute of Aeronautics and Astronautics outside of the shock layer is measured between the discrete and exact solution for the two viscosity formu-
lations. The shock layer is dened to be the distance extending from x = 0 to where the discrete solution is
rst within 0.5% of the exact solution.
The results for low, moderate, and high values of viscosity are shown in Figure 2. At a low viscosity
amplitude, the numerical oscillations in u(x) are damped, but oscillations still remain in the derivative,
u
x
(x), for both solutions. At a higher viscosity amplitude, the Gaussian viscosity solution is smooth for both
u(x) and u
x
(x), but the piecewise-constant viscosity solution still has signicant oscillations in u
x
(x). These
oscillations are due to the conservation of the ux, (u
2
/2 + u
x
), across element boundaries. A jump in (x),
requires a similar jump in u
x
. For higher-order solutions, this jump in the derivative induces uctuations
throughout the element. The greater accuracy of the Gaussian viscosity solution is reected in the H
1
norm values as well. Finally, for much higher viscosity amplitudes, the Gaussian viscosity solution remains
well-behaved, but the piecewise-constant solution suers from oscillations in both u(x) and u
x
(x). At the
high viscosity amplitude, the H
1
norm of the error for the Gaussian viscosity solution is smaller than the
piecewise-constant solution by two or