A Comparison of Arti cial Viscosity Sensors for the Discontinuous Galerkin
Method
by
Joshua Favors
A thesis submitted to the Graduate Faculty of
Auburn University
in partial ful llment of the
requirements for the Degree of
Master of Science
Auburn, Alabama
December 14, 2013
Keywords: Discontinuous Galerkin Method, discontinuity sensor, entropy residual, modal
decay
Copyright 2013 by Joshua Favors
Approved by
Andrew Shelton, Assistant Professor of Aerospace Engineering
Roy J. Hart eld, Jr., Chair, Woltosz Professor of Aerospace Engineering
Brian Thurow, Reed Associate Professor of Aerospace Engineering
Abstract
The following work focuses on the development of an explicit Runge-Kutta, discontin-
uous Galerkin nite element method for the nonlinear, hyperbolic Euler equations. The
discontinuous Galerkin method is a control-volume formulation, such as in the nite-volume
method, but achieves high-order accuracy through the use of high-order polynomial approx-
imations to the computed solution in each discretized element.
Contrary to the Navier-Stokes equations, in which viscosity is present, the Euler equa-
tions do not contain any natural dissipation terms. Numerical methods of second-order or
higher are susceptible to oscillations when discontinuities, such as shocks or large gradient
features form in the solution. The lack of a natural dissipation mechanism requires the use
of a stabilizing method to ensure physical solutions can be obtained.
In order to preserve real physical solutions, the addition of an arti cial viscosity was
applied in the elements where discontinuities or large/sharp gradients of the ow variables
occurred. The addition of arti cial viscosity was accomplished through the use of a sen-
sor that detected where oscillations in the approximated solution occur due to the Gibbs-
phenomenon. The present work incorporates two di erent arti cial viscosity sensors. The
rst sensor is based on the work by Klockner, Warburton, and Hesthaven [1] which looks at
the modal decay rate of the higher-order solution modes. The second sensor was based on
the work presented by Zingan, Guermond, and Popov [2] which uses the entropy production
residual to determine where arti cial viscosity should be applied.
In comparing the two arti cial viscosity sensors, several benchmark problems were mod-
eled: the 2-D isentropic vortex, the 2-D Kelvin-Helmholtz shear layer instability, the 1-D
Sod?s shock tube, the 1-D Planar Shu-Osher problem, and the 2-D shock-vortex interaction.
The results include error analysis, computational cost, robustness, and user-parameters.
ii
Acknowledgments
I have always heard that behind every great man is a strong woman. It must stand to
reason that behind every crazy man is an even stronger and more patient woman. Kelly,
thank you for taking my last name and enduring what I am sure has felt like one never
ending semester for the past few years. You stood beside me as I chased a dream and for
that I will always appreciate what we have accomplished together. I also would like to thank
my daughter, Hannah, who has become the greatest inspiration in my life. You have brought
a joy into our lives that we have never experienced before.
I also owe my deepest thanks to my parents who always supported me throughout my
life. Without your help and con dence in me, the road I have traveled would have been
impassible. Every accomplishment I have made in life has been due to the qualities instilled
through you both. I also must thank my extended family for the life lessons and motivational
support you have provided throughout the years.
The idea of graduate school was always a distant hope, but I must thank you Dr. Shelton
for making it a reality. I appreciate you showing faith in my studies and approaching me
about the possibility of graduate school as I was completing my undergraduate degree. It
has been an honor to study and learn under your guidance.
timshel
iii
Table of Contents
Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ii
Acknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iii
List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vi
List of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ix
List of Abbreviations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . x
1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.2 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
2 Discontinuous Galerkin Method . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
2.1 DG Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
2.1.1 Choice of Basis Function . . . . . . . . . . . . . . . . . . . . . . . . . 11
2.1.2 Numerical Quadrature . . . . . . . . . . . . . . . . . . . . . . . . . . 15
2.1.3 Operator Form . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
2.1.4 Element-wise Operations . . . . . . . . . . . . . . . . . . . . . . . . . 21
3 Numerical Implementation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
3.1 Conservation Laws . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
3.2 Shock-Capturing via Arti cial Viscosity . . . . . . . . . . . . . . . . . . . . . 27
3.3 Normalization of Flow Variables . . . . . . . . . . . . . . . . . . . . . . . . . 29
3.4 Time Evolution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30
3.5 Numerical Flux Function : AUSM+ . . . . . . . . . . . . . . . . . . . . . . . 31
3.6 Arti cial Viscosity Sensors . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
3.6.1 Modal Decay Sensor . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
3.6.2 Entropy Residual Sensor . . . . . . . . . . . . . . . . . . . . . . . . . 37
iv
3.7 Computational E ciency . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
3.7.1 Sparse-Matrix Multiplication . . . . . . . . . . . . . . . . . . . . . . . 42
3.7.2 Open-MP . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43
4 Test Cases . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
4.1 Isentropic Vortex and Kelvin Helmholtz Shear Layer . . . . . . . . . . . . . 45
4.1.1 Isentropic Vortex . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
4.1.2 Kelvin-Helmholtz Shear Layer . . . . . . . . . . . . . . . . . . . . . . 46
4.2 Sod Shock Tube . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
4.3 Shu-Osher Planar Problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48
4.4 Shock-Vortex Interaction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49
5 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51
5.1 Isentropic Vortex and Kelvin Helmholtz Shear Layer . . . . . . . . . . . . . 51
5.1.1 Isentropic Vortex . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51
5.1.2 Kelvin-Helmholtz Shear Layer . . . . . . . . . . . . . . . . . . . . . . 53
5.2 Sod Shock Tube . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58
5.2.1 Case for N = 11, NI = 100 . . . . . . . . . . . . . . . . . . . . . . . . 60
5.3 Shu-Osher Planar Problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65
5.4 Shock-Vortex Interaction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70
6 Conclusions and Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79
6.1 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79
6.2 Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81
Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84
v
List of Figures
1.1 Dissipation and Dispersion in the Transport Equation . . . . . . . . . . . . . . 3
2.1 One-Dimensional Computational Domain . . . . . . . . . . . . . . . . . . . . . 8
2.2 Sketch of Discontinuous Solutions at Element Interfaces . . . . . . . . . . . . . 9
2.3 1-D Transformation from Physical to Standard Element . . . . . . . . . . . . . 11
2.4 Representation of Modal and Nodal Basis for N = 5 . . . . . . . . . . . . . . . 14
3.1 Sparse Sti ness Matrices for N = 9 . . . . . . . . . . . . . . . . . . . . . . . . . 43
3.2 Sparse Lifting Matrices for N = 9 . . . . . . . . . . . . . . . . . . . . . . . . . . 43
3.3 Decrease in Wallclock Time Due to Parallelization . . . . . . . . . . . . . . . . 44
5.1 L2 Density Errors for Approximating Polynomial Order . . . . . . . . . . . . . . 53
5.2 Kelvin-Helmholtz Energy Concentration at t = 61.0043 . . . . . . . . . . . . . . 56
5.3 Kelvin-Helmholtz Arti cial Viscosity at t = 61.0043 . . . . . . . . . . . . . . . . 56
5.4 Kelvin-Helmholtz Energy Concentration at t = 300 . . . . . . . . . . . . . . . . 57
5.5 Kelvin-Helmholtz Arti cial Viscosity at t = 300 . . . . . . . . . . . . . . . . . . 57
5.6 Kelvin-Helmholtz Stats for Arti cial Viscosity Sensors . . . . . . . . . . . . . . 57
5.7 Percentage Distribution of Applied Arti cial Viscosity for Sod Shock Tube . . . 58
vi
5.8 Arti cial Viscosity Statistics per Order of the Approximating Polynomial for Sod
Shock Tube . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59
5.9 Density Pro le of Shock tube (Rarefaction Wave) at Final Time 0.4 . . . . . . . 61
5.10 Density Pro le of Shock tube (Contact Discontinuity) at Final Time 0.4 . . . . 61
5.11 Density Pro le of Shock tube (Shock Wave) at Final Time 0.4 . . . . . . . . . . 61
5.12 Stats for Arti cial Viscosity Sensors for Sod Shock tube . . . . . . . . . . . . . 62
5.13 Time History of Arti cial Viscosity for Sod Shock tube Case . . . . . . . . . . . 63
5.14 Time History of Activated Elements for Sod Shock tube Case . . . . . . . . . . 64
5.15 Percentage Distribution of Applied Arti cial Viscosity for Kelvin-Helmholtz Case 65
5.16 Density Pro le of Shu-Osher Case at Final Time 1.8 . . . . . . . . . . . . . . . 66
5.17 Stats for Arti cial Viscosity Sensors for Shu-Osher Test Case . . . . . . . . . . . 67
5.18 Time History of Arti cial Viscosity for Shu-Osher Case . . . . . . . . . . . . . . 68
5.19 Time History of Elements Activated with Arti cial Viscosity for the Shu-Osher
Case . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69
5.20 Percentage Distribution of Applied Arti cial Viscosity for Shu-Osher Case . . . 70
5.21 Stats for Arti cial Viscosity Sensors for Shock-Vortex Test Case, cmax = 1 . . . 71
5.22 Stats for Arti cial Viscosity Sensors for Shock-Vortex Test Case, cmax = 4 . . . 72
5.23 Percentage Distribution of Applied Arti cial Viscosity for Shock-Vortex Case . . 73
5.24 Density Field for Modal Decay Sensor at Final Time = 16.2 . . . . . . . . . . . 74
vii
5.25 Density Field for Entropy Residual Sensor at Final Time = 16.2 . . . . . . . . 74
5.26 Arti cial Viscosity for Modal Decay Sensor at Final Time = 16.2 . . . . . . . . 75
5.27 Arti cial Viscosity for Entropy Residual Sensor at Final Time = 16.2 . . . . . 76
5.28 Elements Activated with Arti cial Viscosity for Modal Decay Sensor at Final
Time = 16.2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76
5.29 Elements Activated with Arti cial Viscosity for Entropy Residual Sensor at Final
Time = 16.2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77
5.30 Density Field for Entropy Residual Sensor at Final Time = 16.2 . . . . . . . . 78
viii
List of Tables
1.1 RMS Error for Transport Equation Cases . . . . . . . . . . . . . . . . . . . . . 3
5.1 Degrees of Freedom for Grid Discretization . . . . . . . . . . . . . . . . . . . . . 52
5.2 Density L2 Error for Number of Elements, K, and Polynomial Order, N . . . . . 52
5.3 Wall-Clock for Number of Elements, K, and Polynomial Order, N . . . . . . . . 53
5.4 Sod Shock Tube L2-Error for Density Values . . . . . . . . . . . . . . . . . . . . 59
5.5 Wallclock Percent Decrease of Entropy Sensor from Modal Sensor for Sod Shock
Tube . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60
5.6 Percent Di erence of Entropy Sensor Compared to the Modal Sensor . . . . . . 75
ix
List of Abbreviations
bkj basis function
CFD Computational Fluid Dynamics
CFL Courant-Fredrichs-Lewy
FD Finite Di erence
FOU First- Order Upwind
FV Finite Volume
^i, ^j, ^k Cartesian coordinate directions
LGL Legendre-Gauss-Lobatto
Mk transformation map
N approximating polynomial order
NI, NJ, NK number of elements that span each spatial direction ^i, ^j, and ^k
P function space
~um modal coe cient
^un nodal coe cient
RK Runge-Kutta
s entropy/smoothness coe cient
SOU Second-Order Upwind
x
TVD total variation diminishing
V Vandermonde matrix
physical domain
h computational domain
xi
Chapter 1
Introduction
1.1 Motivation
The vast majority of commercial computational uid dynamics (CFD) software used in
the engineering eld today are second-order numerical methods. This is due to the great
amount of research performed over the last few decades to making the nite di erence (FD)
and nite volume (FV) methods highly robust and well understood. As more intensive ow
elds are being analyzed, the need for higher-order methods has gained popularity. The
interest in higher-order methods results from the superior accuracy obtained in a shorter
simulation run-time when compared to a lower-order method achieving the same accuracy
through grid and temporal re nements. The higher-order methods provide higher accuracy
per degree of freedom when compared to classical lower-order methods [3].
The extension to higher-order methods comes at a loss of numerical dissipation that
is present in lower-order numerical approximations. The inherent dissipation in lower-order
methods is due to the even derivative terms in the approximating schemes truncation error.
Since the truncation error in higher-order methods contains higher-order even derivatives, the
dissipation inherent in the numerical scheme decreases as the order of the method increases.
Flow features that were once smoothed by the inherent numerical dissipation in low-order
schemes become more susceptible to producing spurious oscillations in higher-order methods.
These spurious oscillations that occur near discontinuities and sharp gradients are known as
the Gibb?s phenomenon.
To help grasp the idea of the Gibb?s phenomenon, Figure 1.1 presents the numerical
solution of the hyperbolic linear convection equation for a sine wave packet (Figure 1.1a)
and a jump discontinuity (Figure 1.1b). The numerical solution is obtained using two nite
1
di erence approximations known as the rst-order and the second-order upwind methods.
The one-dimensional transport equation is given by
@u
@t +a
@u
@x = 0 (1.1)
where a is the convective velocity which transports a ow variable, u. Two forms of errors
present in the numerical approximation are the dissipation and dispersion errors. If the lead-
ing term in the truncation error contains an even derivative, then the error in the numerical
method will be dominated by dissipative errors which tend to smooth high gradients in the
ow variable. Conversely, if the leading derivative of the truncation error is odd, the approx-
imation will feature dispersive errors which result in spurious oscillations around the large
gradient features. The rst-order upwind (FOU) approximation, with a leading second-order
di erential in the truncation error, is dominated by dissipation errors. Figure 1.1a shows the
e ect of the dissipative error as the amplitude of the wave has decreased compared to the
exact solution. Figure 1.1b shows how the once sharp shock is now smeared, or widened,
due to the dissipative error. The second-order upwind(SOU) approximation contains less
inherent dissipation when compared to the FOU method, but instead is dominated by the
dispersive error. The SOU method given in Figure 1.1a presents a good approximation in
the amplitude of the wave, but minor oscillations occur where the wave meets the constant
solution. Figure 1.1b truly captivates the e ects of the dispersive error. The SOU method,
while maintaining a closer approximation prior to the jump discontinuity, is riddled with
high frequency oscillations immediately following the discontinuity.
Due to the linear nature of the transport equation, as long as the Courant-Fredrichs-
Lewy (CFL) condition is met the approximation remains bounded and will not cause an
unstable approximation to develop. In contrast, oscillations occurring in the highly nonlin-
ear Euler equations of uid motion can quickly become unstable. The resulting instabilities
2
(a) Wave Packet for Transport Equation (b) Jump Discontinuity for Transport Equation
Figure 1.1: Dissipation and Dispersion in the Transport Equation
require a method to keep the unphysical oscillations from destroying the numerical approx-
imation. Another reason of concern occurs in Figure 1.1b where the oscillations cause the
transport quantity to take on negative values. This becomes an issue if, for example, the
Euler equations are being modeled and density takes on a negative value. Negative values
of density or pressure changes the nature of the partial di erential equations locally which
results in the hyperbolic nature of the equation being lost [4]. The approximation cannot be
continued when such an unphysical solution is encountered.
Despite these challenges, the RMS error for the above two cases yields an important
result. Table 1.1 shows that although the SOU method may not be visually appealing, the
method is still more accurate than the FOU method. For the smooth wave packet test case,
the solution is an order of magnitude more accurate. The FOU method is overly dissipative
due to the implicit arti cial viscosity that is inherent in the numerical method, while the
SOU method develops spurious oscillations due to the lack of implicit arti cial viscosity.
Case Wave Packet Jump Discontinuity
FOU 5.074E-2 3.943E-2
SOU 4.710E-3 3.078E-2
Table 1.1: RMS Error for Transport Equation Cases
3
Using the previous results, one approach in dealing with the oscillations is to add explicit
arti cial viscosity in order to stabilize the approximation while still maintaining higher-order
accuracy. The explicit arti cial viscosity, represented by a dissipative term, is added to the
governing equation and proves to be a robust method for higher-order methods. The arti cial
viscosity will dampen the oscillations such that the discontinuity or sharp gradient can be
properly resolved in the numerical method. The addition of the arti cial viscosity controls the
spurious oscillations and will prevent the numerical approximation from becoming unstable
or taking on unphysical values. Unlike the implicit arti cial viscosity that is applied without
bias throughout the entire discretized domain, the explicit arti cial viscosity should be added
only in the regions where discontinuous or sharp gradients threaten the stability of the
numerical method in order to ensure optimal accuracy. In order to apply arti cial viscosity
to the elements containing the spurious oscillations, a sensor determines where in the ow
domain dissipation is required. Robustness, accuracy, and computational expense become
important characteristics when determining which sensor to use.
1.2 Background
The Euler equations of uid motion, when coupled with an appropriate equation of
state, describe the conservation laws for an inviscid compressible ow. The equations are a
simpli cation of the Navier-Stokes equations when viscosity and thermal conduction terms
are neglected. The nonlinear hyperbolic Euler equations are useful for approximating ows
where the e ects of viscosity are negligible, as in high-speed ows and regions outside of
the boundary layer. When approximating the conservative form of the Euler equations, ow
discontinuities such as shock waves and contact discontinuities may develop in the solution.
Unlike the full Navier-Stokes equations, the Euler equations do not contain any natural
di usion operators. Arti cial viscosity is the inclusion of dissipating terms in a numerical
method, whether caused by the numerical scheme itself (implicit arti cial viscosity) or by a
purposely added term for numerical stability (explicit arti cial viscosity). Arti cial viscosity
4
exhibits a di usive behavior and tends to reduce all gradients in the ow, whether physically
correct or as a result of the numerical scheme [5]. Since the previously discussed Gibb?s
phenomenon lead to dispersive errors resulting in spurious oscillations, arti cial viscosity
can be used to selectively dampen the nonphysical approximations.
The DG method is a higher-order method that is closely related to the FV method.
Whereas the FV method uses a piecewise linear approximation across each discretized cell,
the DG method instead approximates the solution using a polynomial representation inside
each discretized element. The DG method is capable of hp-adaptivity: increased accuracy
of the numerical method can be accomplished by increasing the number of approximat-
ing elements, h-adaptivity, or by increasing the order of the approximating polynomial,
p-adaptivity.
The addition of arti cial viscosity is a robust method to contain the spurious oscillation
that occur in higher-order schemes such as the DG method. Arti cial viscosity di uses a
discontinuity or sharp gradient to the extent that the numerical method may approximate
the solution without the resulting Gibb?s phenomenon causing instabilities. While numerical
methods approximate shocks as discontinuities, shocks occur in nature over a length scale on
the order of the mean-free-path of the medium in which it occurs. Discretizing a ow eld to
accurately capture a shock wave would require computing times that become unrealistic. The
role of arti cial viscosity is to smear the sharp pro le over such a length that the numerical
method can properly approximate the discretized ow eld condition.
Other methods that have been used to stabilize numerical schemes include reducing the
order of the approximation down to rst-order in the region of sharp gradients via the use
of limiters or non-oscillatory polynomial reconstructions (ENO,WENO). Reducing the order
of the approximating scheme down to rst-order introduces the large amount of implicit
arti cial viscosity into the solution. This approach defeats the idea of using a higher-order
method as the solution is reduced to rst-order in the vicinity of Gibb?s phenomenon. The
order of the applied arti cial viscosity scales like O(h) for rst-order methods. In contrast,
5
the needed explicit arti cial viscosity added to a high-order method to properly resolve a
shock pro le is ofO(h=p) [6]. Unless mesh re nement is coupled with the rst-order reduction
in the vicinity of troubled cells, the accuracy of the solution may become severely degraded.
Limiters, while the preferred method in the FV setting, have been found to only work in the
DG scheme for structured, one-dimensional discretizations. The limiter methods developed
for DG schemes are usually dependent upon grid structure, polynomial order, and the choice
of quadrature sets [7]. This dependency reduces the overall robustness of the method and
does not allow a general approach to solving a wide range of ow situations. Limiters also
reduce the approximation to linear or constant in a wide region around the trouble elements
[8]. If high-order approximations are sought then mesh re nement must then be carried out
in the regions where limiting is occurring. Weighted essentially non-oscillatory (WENO)
approaches have garnered interest due to their ability to retain higher-order solutions near
oscillatory regions. This is accomplished by using additional degrees of freedom to resolve
the sharp pro les which results in preserving the nonlinear stability [8]. Unfortunately the
computational expense becomes increasingly excessive for higher-order methods. In addition,
the methods have only proven compatible with one-dimensional structured grids.
In this work, the implementation of the arti cial viscosity was constructed using a sen-
sor to detect where sharp gradients or discontinuities were occurring in the approximated
solution. Two separate sensors were analyzed to determine accuracy, robustness, and com-
putational expense. The rst sensor analyzed uses the modal expansion of the approximate
solution inside each discretized element. The modal decay of the expansion is de ned as the
decreasing magnitude of the modal coe cients as the mode number increases. This method is
based on the resulting numerical behavior of the modal coe cients of the solution. The sen-
sor approximates the modal decay of the solution using an exponential equation,j^unj n s.
The resulting estimated decay exponent, s, is then used to determine if arti cial viscosity
should be applied to the approximated solutions in each element. The second method uses
the residual of an associated entropy equation to determine where arti cial viscosity should
6
be applied. The residual of the entropy equation takes a value less than zero when under-
resolved features, such as shocks or steep gradients in the ow variables, are encountered.
For smooth regions, the residual is nonnegative. The arti cial viscosity sensor is then built
on the maximum magnitude of the residual in each element. Contrary to the modal sensor,
the entropy sensor is based on the physics of the ow eld by applying arti cial viscosity
when the Second Law of Thermodynamics is being violated.
Several test cases will be presented in one and two-dimensions. The 2-D isentropic vortex
case presents a smooth solution where, given adequate resolution, no activation of arti cial
viscosity should occur. The 2-D Kelvin-Helmholtz test case presents a developing shear
layer originating from a horizontal jet. An instability in the ow is created by small vertical
velocity perturbations that are added onto the horizontal jet. The shock-free ow eld is still
prone to numerical instabilities due to the sharp ow gradient present in the shear layer.
The ability of each sensor to detect and di use these sharp regions is examined. The classic
1-D Sod shock tube is presented due to the ability to compare the numerical results with a
known exact solution. The resulting ow eld is shown to develop a shock wave, a contact
discontinuity, and a rarefaction wave. The shock-vortex interaction and Shu-Osher 2-D test
cases present ow elds that incorporate a shock wave interacting with small-scale features.
These test cases give insight into the ability of each scheme to di use the shock wave while
not overly dampening the small-scale structures that occur in the same proximity.
7
Chapter 2
Discontinuous Galerkin Method
2.1 DG Method
The following analysis of the DG method follows the work of Hesthaven and Warburton
in [4]. Consider a physical domain, , that is discretized using k = 1;:::;K non-overlapping,
boundary conforming elements, Dk. For the structured meshes used in this work, the total
number of elements is given by K = (NI NJ NK), where NI, NJ, and NK are the number
of elements that span each spatial direction, ^i, ^j, and ^k. A one-dimensional schematic is
shown in Figure 2.1. The physical domain is approximated by the discretized computational
domain, h.
Figure 2.1: One-Dimensional Computational Domain
To assist in the explanation of the DG method, a general one-dimensional scalar con-
servation law of the following form is speci ed.
@u
@t +
@f(u)
@x = 0; x2 (2.1)
The global solution, u(x;t), is approximated by the direct sum of each element?s piecewise
continuous, N-th order polynomial approximations, uh(x;t), and is represented as
u(x;t)?uh(x;t) =
KM
k=1
ukh(x;t) (2.2)
8
Let bkj represent a high-order local basis function and ~ukj represents the modal coe cients,
then the approximate local solution can be given by
ukh(x;t) =
NX
j=0
~ukj(t)bkj(x) (2.3)
The local residual of the approximate solution can be formed on each element by
Rkh(x;t) = @u
k
h
@t +
@f(ukh)
@x ; x2D
k (2.4)
The residual is then required to vanish locally in a Galerkin sense by performing a weighted
average. The Galerkin method speci es that the weighting function is chosen from the same
basis as the state. The foundation of determining the local solution is then given by
Z
Dk
Rkhbki dx = 0 i = 0;:::;N (2.5)
By requiring the residual to vanish on each element, the elemental local solutions are discon-
nected from the local solutions on adjacent elements. This results in double-valued solutions
at each element interface as shown in Figure 2.2.
Figure 2.2: Sketch of Discontinuous Solutions at Element Interfaces
In order to obtain a global solution, Equation 2.5 is rst integrated by parts and then the
divergence theorem is applied resulting in
Z
Dk
Rkhbki dx =
Z
Dk
bki @f(u
k
h)
@t dx+
Z
Dk
bki @u
k
h
@x dx = 0 (2.6)
9
Z
Dk
bki @u
k
h
@t dx+
Z
Dk
@
@x
bkif(ukh)
dx
Z
Dk
@bki
@xf(u
k
h) dx = 0 (2.7)
Z
Dk
bki @u
k
h
@t dx+
I
@Dk
bkif(ukh)nkx ds
Z
Dk
@bki
@xf(u
k
h) dx = 0 (2.8)
where @Dk denotes the element edge.
The term f(ukh)nkx presented in Equation 2.8 represents the non-unique ux that oc-
curs at the element interfaces. This term can be replaced with a physics-based numerical
ux function, f , that approximates the physical ux through the solution of the Riemann
problem. Many numerical ux functions have been developed by the FV community for use
on uid ow problems. The numerical ux function is de ned in terms of interior (-) and
exterior (+) interface states and is given by
x2@Dk; fkhnkx!f (uk; h ;uk;+h ;nk; x ) (2.9)
The minimum requirements imposed for consistency are then given by
f (uk; h ;uk;+h ;nk; x ) = f (uk;+h ;uk; h ;nk;+x )
f (u;u;nx) f(u)nx
(2.10)
The rst equation states that any ux leaving the right interface of an element, is the negative
of the ux leaving the left interface of the adjacent element. This ensures a conservation of
the ux between two adjacent cells. The second condition ensures that the numerical ux
of two identical states is simply equal to the ux. One convenient aspect of the DG method
is that the boundary conditions are simply applied through the uk;+h term in the numerical
ux function. The stability of the DG method is enforced through the local ux choice [4].
10
With the numerical ux function inserted, the approximating semi-discrete DG scheme
takes the form
Z
Dk
bki @u
k
h
@t dx+
I
@Dk
bkif (uk; h ;uk;+h ;nk; x )ds
Z
Dk
@bki
@xf
k
hdx = 0
i = 0;:::;N; k = 1;:::;K
(2.11)
The solution of Equation 2.11 are piecewise smooth inside each element, but at the element
interfaces the solution is discontinuous. The addition of the numerical ux function allows
the discontinuous elements to communicate and the global solution can be obtained from
each individual element solution.
2.1.1 Choice of Basis Function
Consider a standard element such that x! 2( 1;1), and expand the state using a
degree N polynomial function space such that, uh( )2PN. A transformation map, Mk, is
used to get from the standard element to each physical element, Dk, by
x2Dk : x =Mk( ) (2.12)
and is shown in Figure 2.3.
Figure 2.3: 1-D Transformation from Physical to Standard Element
Using a linear interpolation, the transformation map is given by
Mk( ) = xka + 1 + 2 (xkb xka) (2.13)
11
The choice of the local basis is important because it is the local basis which gives the DG
method it?s accuracy. At each element, the local solution can be represented as a Nth-order
polynomial by either a modal (2.14) or nodal basis (2.15).
uh( ) =
NX
m=0
~um m( ) (2.14)
uh( ) =
NX
n=0
^un?n( ) (2.15)
where m is a local polynomial basis, and ?n is a Lagrange interpolating polynomial. The
weighting functions are also selected from either the modal class, bi! i, or the nodal class,
bi!?i, where m;?n2PN.
Modal Basis
The modal coe cients, ~um, can be found using an L2 projection,
Z 1
1
i( )u( )d =
Z 1
1
i( )uh( )d
=
Z 1
1
i( )
NX
m=0
~um m(
!
d
=
NX
m=0
Z 1
1
i( ) m( )d
~um
=
NX
m=0
Mim~um
(2.16)
where Mim is the mass matrix. To recover the modal coe cients, the inverse of the mass
matrix must be computed. An orthonormal basis is used to insure that the mass matrix is
well-conditioned in order to preserve accuracy. The three term recurrence relation is used
12
to nd an orthonormal basis and is given by [9]
p
m+1 m+1( ) = ( m) m( )
p
m 1 m 1( )
1 := 0; 0 := 1=
p
0
(2.17)
for m = 0;:::;N 1. By choosing the coe cients ( m; m), di erent types of polynomials
can be obtained from Equation 2.17 such as the Legendre, Chebychev, or Hermite polyno-
mials. Following the work of Hesthaven [4], the Legendre polynomials were chosen to form
a hierarchical modal basis. The coe cients that result in the Legendre polynomials are
m = 0; m =
8
><
>:
2 if m = 0
1=(4 m 2) if m 1
(2.18)
An important property, resulting from the use of Legendre polynomials, is given by
Z 1
1
i( ) j( )d = ij :=
8
><
>:
1 if i = j
0 if i6= j
(2.19)
Nodal Basis
The nodal basis, given in Equation 2.15, uses the Lagrange interpolating polynomials.
The resulting polynomial takes the form
?n( ) =
NY
i=0i6=n
i
n i (2.20)
for the nodal locations given by o;:::; N. The important property resulting from using the
Lagrange polynomials is given by
?j( i) = ij :=
8
><
>:
1 if i = j
0 if i6= j
(2.21)
13
Modal/Nodal Relationship
While the modal and nodal representations are both mathematically equivalent, they are
computationally di erent. Figure 2.4 shows the representations of the modal and nodal basis
for a fth-order polynomial. The modal basis, using Legendre polynomials, is a hierarchical
basis. This implies that an Nth-order approximation will use N + 1 curves of increasing
polynomial order. Figure 2.4a shows a 5th-order modal approximation that is composed
of a constant line and then the curves of polynomial orders 1;2;:::;5. When computing
a solution in an element, the modal representation will use the combination of the N + 1
curves to obtain a solution on the element. The nodal representation will nd a solution
by generating N + 1 curves of each polynomial order N. The nodal representation then
computes an approximation for only the node whose solution is being required and zeros at
the other node locations.
?1 0 1
?2
0
2pi
0
?1 0 1
?2
0
2pi
1
?1 0 1
?2
0
2pi
2
?1 0 1
?2
0
2pi
3
?1 0 1
?2
0
2pi
4
?1 0 1
?2
0
2pi
5
(a) A modal basis with m( )2P5
?1 0 1?1
0
1l
0
?1 0 1?1
0
1l
1
?1 0 1?1
0
1l
2
?1 0 1?1
0
1l
3
?1 0 1?1
0
1l
4
?1 0 1?1
0
1l
5
(b) A nodal basis with ?n( )2P5
Figure 2.4: Representation of Modal and Nodal Basis for N = 5
A connection between the modes, ~u, and the nodal values, ^u, can be made by de ning the
Vandermonde matrix (V). The Vandermonde matrix is formed by evaluating the Legendre
polynomials at the same location that de ne the Lagrange polynomials. The Vandermonde
14
matrix is then evaluated as
Vij := j( i) (2.22)
Substituting the Vandermonde matrix into the modal and nodal representations yields
uh( i) =
NX
m=0
~um m( i) =
NX
n=0
^un?n( i)()
NX
m=0
Vim~um =
NX
n=0
in^un (2.23)
Rewriting these in matrix-vector form results in
~uh := uh(~ ) = ^u = V ~u () ~u = V 1^u (2.24)
In order to successfully transition between the modal and nodal values the Vandermonde
matrix must be well-conditioned. The conditioning of the Vandermonde matrix is dependent
on the choice of the nodal interpolation points. The Lebesque constant, given by
= maxr
NpX
i=1
j?i( )j (2.25)
indicates how far away from the interpolation may from the best possible polynomial representation[4].
Minimizing the Lebesque constant will result in the best possible interpolation points so that
the Vandermonde matrix is well-conditioned. The resulting points are the Legendre-Gauss-
Lobatto (LGL) quadrature points [10].
2.1.2 Numerical Quadrature
Following the work by Gautschi [9], the integration of a function, f, may be approxi-
mated via summation by Z
1
1
f( )d =
NpX
j=1
!jf( j) (2.26)
where !j are the quadrature weights and j are the quadrature nodes. A key property of
polynomials is that Gaussian quadrature is exact.
15
The Gauss rule, which does not account for the interface nodes since f2P2Np 1 (with
1 < 1; Np < 1), is given by
Z 1
1
f( )d ?
NpX
i=1
!if( i) (2.27)
The nodes and weights are found by using the polynomial recursion formula. Begin by letting
~ ( ) := [ 0( ); 1( );:::; p 1( )]T and de ne a Jacobi matrix
Jp :=
2
66
66
66
64
0 p 1 0
p
1 1
p
2
...
0 p p 1 p 1
3
77
77
77
75
(2.28)
The rst p terms of the recursion in matrix-vector form is
~ ( ) = Jp~ ( ) +p p p( )~ep (2.29)
Finding the eigenvalues, i, and the corresponding eigenvectors, ~v(i), of the Jacobi matrix
Jp, then the nodes and weights of the p-point Gauss formula are
i = i; !i = 0(v(i)1)2; i = 1;:::;p (2.30)
with 1 < 1 and Np < 1.
The Gauss-Lobatto rule, which includes the interface nodes since f 2 P2Np 3 (with
0 = 1; Np+1 = 1), is given by
Z 1
1
f( )d = !0f( 0) +
NpX
i=1
!if( i) +!Np+1f( Np+1) (2.31)
16
The polynomial recursion formula is also used to nd the nodes and weights for the Gauss-
Lobatto rule. A Jacobi-Lobatto matrix is given by
JLp+2 :=
2
64 Jp+1
p
p+1~ep+1
p
p+1~e
T
p+1
p+1
3
75 (2.32)
where for Legendre polynomials the coe cients are given as
p+1 = 0; p+1 = 1=2 + 1=(4p+ 2) (2.33)
Finding the eigenvalues, Li , and the corresponding eigenvectors, ~vL(i), of the Jacobi-
Lobatto matrix JLp+2, then the nodes and weights of the p+ 2-point Gauss-Lobatto formula
are
i = Li ; !i = 0(vL(i)1)2; i = 0;:::;p+ 1 (2.34)
with 0 = 1 and Np+1 = 1.
2.1.3 Operator Form
Restating the model problem and DG discretization:
@u
@t +
@f(u)
@x = 0 (2.35)
Z
Dk
bki @u
k
h
@t dx
Z
Dk
@bki
@xf
k
hdx+
I
@Dk
bkif (uk; h ;uk;+h ;nk; x )ds = 0 (2.36)
In order to evaluate the integrals on the physical elements, a transformation for the di er-
entials is given by
x2Dk : dx =Jkd ; Jk = x
k
b x
k
a
2 (2.37)
17
where Jk is the Jacobian of the transformation. The state and ux decomposition is found
by using the coordinate transformation so that
ukh(x) = uh(Mk( )) =
NX
j=0
~ukj(t)bj( ) (2.38)
fkh(x) =
NX
j=0
~fkj (t)bkj( ) (2.39)
where M is the transformation map de ned in Equation 2.12.
The time derivative present in Equation 2.36 are now evaluated using the modal expan-
sion given in Equation 2.3:
Z
Dk
bki @u
k
h
@t dx =
Z
Dk
bki (x) @@t
NX
j=0
~ukj(t)bkj(x)
!
dx
=
NX
j=0
Z
Dk
bki (x)bkj(x)dx
d~uk
j
dt
=
NX
j=0
Z 1
1
bi( )bj( )Jkd
d~uk
j
dt
=Jk
NX
j=0
Mijd~u
k
j
dt
(2.40)
where Mij is the mass matrix de ned by
Mij :=
Z 1
1
bki ( )bkj( )d (2.41)
18
The spatial derivative term goes through the same process as the time derivative term, and
is given by
Z
Dk
dbki
dxf
k
hdx =
Z
Dk
bki (x)db
k
i (x)
dx
NX
j=0
~fkjbkj(x)
!
dx
=
NX
j=0
Z
Dk
dbki (x)
dx b
k
j(x)dx
~fkj
=
NX
j=0
Z 1
1
dbi( )
d bj( )d
~fkj
=
NX
j=0
(ST)ij ~fkj
(2.42)
where Sij is the sti ness matrix and is de ned by
STij :=
Z 1
1
bi( )dbj( )d d (2.43)
Noticing the closed interval, the numerical ux, f , does not exist over the whole element,
Dk, but only on the interface @Dk. The numerical ux is expanded out regardless as
f (s) = f (uk; h ;uk;+h ;nk; x )
s
=
NX
j=0
g(f s)k
jbj(s) (2.44)
19
where the variable, s, simply denotes the restriction x2@Dk. The remaining formulation is
given by
I
@Dk
bkif (uk; h ;uk;+h ;nk; x )ds
= [bkif (uk; h ;uk;+h ;nk; x )]xka + [bkif (uk; h ;uk;+h ;nk; x )]xk
b
= bki (xka)
NX
m=0
g(f a)k
mb
k
m(x
k
a) +b
k
i (x
k
b)
NX
n=0
g(f b )k
nb
k
n(x
k
b)
=
NX
m=0
bk
i (x
k
a)b
k
m(x
k
a)
g(f
a)
k
m +
NX
n=0
bk
i (x
k
b)b
k
n(x
k
b)
g(f
b )
k
n
=
NX
m=0
(bi( 1)bm( 1)) g(f a)km +
NX
n=0
(bi(1)bn(1)) g(f b )kn
=
NX
m=0
(La)img(f a)km +
NX
n=0
(Lb)ing(f b )kn
(2.45)
The new variables (La)ij and (Lb)ij are the lifting matrices. Since the solution is approx-
imated by a Nth-order polynomial, the data obtained at the boundaries must be extended
throughout the rest of the domain. These matrices \lift" the lower dimensional boundary
data to the higher dimensional element data. The lifting matrices are given by
(La)ij := bi( 1)bm( 1); (Lb)ij := bi(1)bn(1) (2.46)
where (La)ij corresponds to the left interface of the standard element and (La)ij corresponds
to the right interface of the standard element With the basis expansion complete, Equation
2.36 can now be expressed as
Jk
NX
j=0
Mijd~u
k
j
dt
NX
p=0
(ST)ip ~fkp +
NX
m=0
(La)img(f a)km +
NX
n=0
(Lb)ing(f b )kn = 0 (2.47)
20
The DG discretization is further simpli ed into the matrix vector form by using the mass,
sti ness, and lifting matrices:
JkMd~u
k
dt S
T ~fk +Lag(f
a)
k +L
bg(f b )
k = 0 (2.48)
Equation 2.48 represents a ODE in time after the spatial discretization has been carried out.
This nal form is the method implemented in the current computational work.
2.1.4 Element-wise Operations
The current work employed the nodal basis for the computational implementation of the
DG method. At element interfaces, the nodal approach is more e cient due to the presence
of nodal values at the interface locations. The modal approach would require extrapolation
in order to obtain the interface values needed for the numerical ux function. The nodal
approach, which outputs point-wise physical values, is more convenient to deal with when
the ux is a nonlinear function. The advantage of the modal approach lies in it having the
property of hierarchical construction. This would be a bene t if resolution adaptivity where
to be implemented.
Using the same methods that allow the representation of the function uh by,
uh( i) =
NX
j=0
~uj j( ) =
NX
i=0
^ui?i( ) (2.49)
then the same method is applicable to the basis functions resulting in the expression
j( ) =
NX
i=0
d( j)
i?i( ) (2.50)
The nodal coe cient is then simply the pointwise value of the function:
d( j)
i = j( i)
:= V
ij (2.51)
21
where Vij is again the Vandermonde matrix. The relationship between the nodal and the
modal basis functions are given by
j( ) =
NX
i=0
Vij?i( ) =
NX
i=0
(VT)ji?i( )()?i( ) =
NX
j=0
(V T)ij j( ) (2.52)
To formulate the DG method using the nodal expression relies on the relationship
?i( ) =
NX
j=0
(V T)ij j( ) (2.53)
Mij :=
Z 1
1
?i( )?j( )d
=
Z 1
1
NX
m=0
(V T)im m( )
NX
n=0
(V T)jn n( )d
=
NX
m=0
NX
n=0
(V T)im
Z 1
1
m( ) n( )d
(V 1)nj
=
NX
m=0
NX
n=0
(V T)im mn(V 1)nj
=
NX
m=0
NX
n=0
(V T)im(V 1)mj () M = (VVT) 1
(2.54)
In order to express the sti ness matrix in a nodal form, rst a di erentiation matrix that
transforms a point value to the derivative at that point is de ned by
Dij := ?0i( j)
=
NX
m=0
(V T)jm 0m( i)
=
NX
m=0
(V T)jm(V )im
=
NX
m=0
(V )im(V 1)mj () D = (V V 1)
(2.55)
22
The sti ness matrix in a nodal basis is given by
Sij :=
Z 1
1
?i( )?0j( )d () (ST)ij :=
Z 1
1
?0i( )?j( )d (2.56)
Using the mass matrix with the di erentiation matrix de ned in Equation 2.55 the sti ness
matrix can be formulated by
NX
n=0
MijDnj =
NX
n=0
Z 1
1
?i( )?n( )d
?0j( n)
=
Z 1
1
?i( )
NX
n=0
?0j( n)?n( )
!
d
=
Z 1
1
?i( )
NX
n=0
d ?0
j
n?n( )
!
d
=
Z 1
1
?i( )?0j( )d
= Sij () S = (MD)
(2.57)
Since the transpose of the sti ness matrix is used in the DG implementation, the formula
for the transpose of the sti ness matrix is given by
S = (MD) () ST = (DTMT)
= (V V 1)T(V TV 1)T
= (V V 1)T(V TV 1)
= (V V 1)T(VVT) 1
(2.58)
In order to incorporate the solution of the numerical ux function at the interface throughout
the element, the lifting matrices are formulated for the two element interfaces given by \a"
and \b". The \a" interface corresponds to the standard element left interface location of 1
and the \b" interface corresponds to the standard element right interface location of 1. The
23
formulation of the \a" and \b" interfaces follow the same procedure by
(La)ij := ?i( 1)?j( 1)
= ?i( 0)?j( 0)
= 0i 0j
(Lb)ij := ?i(1)?j(1)
= ?i( N)?j( N)
= Ni Nj
(2.59)
De ning the numerical ux function using the interface locations of \a" and \b", their
formulation is given by
(f a)k =
NX
j=0
d(f a)k
j?j( 1)
=
NX
j=0
d(f a)k
j?j( 0)
=
NX
j=0
d(f a)k
j 0j
= d(f a)k0
(f b )k =
NX
j=0
d(f b )k
j?j(1)
=
NX
j=0
d(f b )k
j?j( N)
=
NX
j=0
d(f b )k
j Nj
= d(f b )kN
(2.60)
Combining the results for the lifting matrix and the ux function at the \a" interface results
in
NX
i=0
(La)ijd(f a)kj =
NX
i=0
0i 0j 0jd(f a)kj =
NX
i=0
i0d(f a)k0 (2.61)
Lad(f a)k =
2
66
66
66
64
1 0 ::: 0
0 0 ::: 0
... ... ...
0 ::: 0
3
77
77
77
75
8
>>>>
>>>
<
>>>
>>>
>:
d(f a)k
0
0
...
0
9
>>>>
>>>
=
>>>
>>>
>;
= d(f a)k0
8
>>>>
>>>
<
>>>
>>>
>:
1
0
...
0
9
>>>>
>>>
=
>>>
>>>
>;
= d(f a)k0~e0
Following the same method as above, the \b" interface becomes
NX
i=0
(Lb)ijd(f b )kj =
NX
i=0
Ni Nj Njd(f b )kj =
NX
i=0
iNd(f b )kN (2.62)
24
Lbd(f b )k =
2
66
66
66
64
0 0 ::: 0
0 0 ::: 0
... ... ...
0 ::: 1
3
77
77
77
75
8
>>>>
>>>
<
>>>
>>>
>:
0
...
0
d(f b )k
N
9
>>>>
>>>
=
>>>
>>>
>;
= d(f b )kN
8>
>>>
>>>
<
>>>
>>>
>:
0
...
0
1
9>
>>>
>>>
=
>>>
>>>
>;
= d(f b )kN~eN
With the results from the above formulations, the local form of the 1-D DG statement
in using a nodal basis becomes
JkMd^u
k
dt S
T ^fk +La(f
a)
k +Lb(f
b )
k = 0 (2.63)
and resummarizing the components as
Mass Matrix M = (VVT) 1
Sti ness Martrix ST = (V V 1)T(VVT) 1
Lifting Martrices - (\a") La = ~e0
Lifting Martrices - (\b") Lb = ~eN
Vandermonde Martrix Vij := j( i)
Derivative Vandermonde Martrix (V )ij := 0j( i)
(2.64)
25
Chapter 3
Numerical Implementation
3.1 Conservation Laws
The DG code developed in this work was used to approximate the Euler equations. The
Euler equations are the non-linear, hyperbolic equations that describe the conservation of
mass, momentum, and energy using the continuum hypothesis for uid ows. The Euler
equations are a simpli cation of the Navier-Stokes equation in which viscosity and thermal
conduction terms are neglected. The assumption of inviscid ow is a standard simpli ca-
tion for high-speed ows in which the e ect of viscosity becomes increasingly negligible with
increasing Reynolds number. The challenge of numerically approximating the Euler equa-
tions is the occurrence of localized discontinuities, such as shocks and contact discontinuities,
which may occur in the approximated solution.
The Euler equations, which express the conservation of mass (Eq. 3.1), momentum (Eq.
3.2 & 3.3), and energy (Eq. 3.4), are given in two dimensions as
@
@t +
@( u)
@x +
@( v)
@y = 0 (3.1)
@( u)
@t +
@( u2 +p)
@x +
@( uv)
@y = 0 (3.2)
@( v)
@t +
@( uv)
@x +
@( v2 +p)
@y = 0 (3.3)
@( e0)
@t +
@u( e0 +p)
@x +
@v( e0 +p)
@y = 0 (3.4)
26
where the conserved variables are density, , momentum, u, and the total energy of the gas,
e0. Rearranging the equations in ux vector (strong conservation) form results in
@U
@t +
@
@xF(U) +
@
@yG(U) = 0 (3.5)
The state vector, U, and the inviscid ux vectors, F and G, are expressed by
U =
8>
>>>
>>>
<
>>>
>>>
>:
u
v
e0
9>
>>>
>>>=
>>>
>>>
>;
F =
8>
>>>
>>><
>>>
>>>>
:
u
u2 +p
vu
( e0 +p)u
9>
>>>
>>>=
>>>
>>>>
;
G =
8>
>>>
>>><
>>>
>>>>
:
v
uv
v2 +p
( e0 +p)v
9
>>>
>>>
>=
>>>
>>>>
;
Noting that the ratio of speci c heats is de ned by = cp=cv, the system of equations can
be closed with the ideal gas equation of state
p = ( 1)
e0 u
2 +v2
2
(3.6)
3.2 Shock-Capturing via Arti cial Viscosity
In the approach of shock-capturing methods, shock waves are computed as part of the
solution using the conservative form of the Euler equations. The conservative form of the
Euler equations are required so that the solution will satisfy the Rankine-Hugoniot equations
[5].
When shock-capturing methods are used, high frequency solutions will form near shock
waves and other high gradient features. The high frequency solutions are caused by the
Gibbs phenomena, and a method must be used to keep the solution from taking on unreal
values. In this work, arti cial viscosity is used to smooth out the high frequency areas in
the solution. The resulting piecewise discontinuous arti cial viscosity that is applied to the
needed elements is accomplished with the addition of an unphysical di usion term to the
27
Euler equations such that
@U
@t +
@
@xF(U) +
@
@yG(U) =r ( rU) (3.7)
The resulting di usive scheme is numerically consistent with damping methods provided by
some schemes. Expanding out the viscosity term on the right hand side yields
@U
@t +
@
@xF(U) +
@
@yG(U) =
@
@x
@U@x
+ @@y
@U@y
(3.8)
where the arti cial viscosity is given by .
Due to the unbiased nature that arti cial viscosity possesses, an appropriate method is
required to ensure that the di usion is applied to the unphysical high-frequency solutions
instead of the physically correct features present in the ow eld. With the possibility of
arti cial viscosity dampening relevant information, the method in which di usion is applied
becomes vital to ensuring that the arti cial viscosity is applied only in areas where it is
required. Although the Gibb?s phenomenon tends to produce oscillations around sharp
gradients, the accuracy of the approximation is still preserved. The key to ensuring a robust
method while still maintaining accuracy is to apply the smallest amount of arti cial viscosity
such that the solution remains stable.
With the addition of the arti cial viscosity components given in Equation 3.8, the
numerical ux function used for the inviscid ux will be avoided for these terms. When
the numerical ux function encounters the arti cial viscosity terms in Equation 3.8, simple
averages of the adjacent states are used in place of the numerical ux function and is shown
by
U = U
L +UR
2
Fvis = Fartvis(U
L) +Fartvis(UR)
2
(3.9)
28
The present work uses a simple nonlinear arti cial ux, r ( rU), applied to each
state. It is worth noting that other methods can be used to apply the arti cial viscosity
to the conservation equations. Zingan [11] uses the following viscous ux to augment the
inviscid ux, and is given by
g(c) =
8
>>>
><
>>>>
:
r
rsu
rsu u rT
9
>>>
>=
>>>>
;
(3.10)
where rsu is a symmetric tensor of rank 2 given by
(rsu)ij = 12
@u
i
@xj +
@uj
@xi
(3.11)
The new arti cial viscosity terms are a di usion viscosity, a dynamic viscosity, and a
thermal viscosity. The temperature, T, is found using the standard Equation of State for a
perfect gas.
3.3 Normalization of Flow Variables
The ow variables present in the Euler equations are normalized in this work for con-
venience. The normalization assumes a uniform gas constant, R = R1, and speci c heat
ratio, = 1. The choice of normalization used in this work is
xy = xL; ty = ta1L ; uy = ua
1
y =
1
; py = p
1a21
= p p
1
; Ty = TR1a2
1
= T T
1
ey = ea2
1
; cyp = cpR
1
= 1; cyv = cvR
1
= 1 1 (3.12)
29
By normalizing the ow variables in this way, the ideal gas equation and the speed of sound
in an ideal gas remain the same and are given by
py = ( 1) yey; ay =
s
py
y ; e
y = cy
vT
y (3.13)
Any reference to these ow variables will assume the normalized form and the normalization
notation (:)y will be abandoned throughout the following.
3.4 Time Evolution
Determining the timestep when solving the Euler equations, the Courant-Friedrichs-
Lewy (CFL) criteria is usually employed. Due to the e ects of arti cial viscosity, a more
restrictive formulation is employed to ensure stability. Using the same method as [4], the
time-step is determined by
t 1
maxN
2
h +k kL1
N4
h2
(3.14)
where max refers to the largest global characteristic velocity, is the arti cial viscosity, and
h is the smallest local cell size. The determination of the timestep using Equation 3.14 was
found to be overly restrictive for some of the test cases due to the sensor never applying the
maximum arti cial viscosity throughout the solution. In order to ensure continuity in the
sensor comparison, Equation 3.14 was used for every method regardless.
The simple forward Euler time stepping is only rst-order accurate in time and is given
by
un+1 = un + tR(un) (3.15)
where R is given by
du
dt =R(u) (3.16)
Since the DG method focuses on obtaining higher-order spatial accuracy, the time integration
should also employ a method that obtains higher-order accuracy. Because discontinuities
30
are being modeled with the Euler equations in this work, it is important that the temporal
scheme does not also admit spurious oscillations into the solution. The current work employs
a three-stage, third-order Runge-Kutta (RK) method that maintains the TVD property [12].
The TVD property ensures monotonic solutions, thus eliminating the spurious oscillations
caused by the time evolution. The 3rd-order TVD-RK method is given by
u(1) = un + tR(u(n))
u(2) =
3
4
un +
1
4
u(1) +
1
4
tR(u(1))
un+1 =
1
3
un +
2
3
u(2) +
2
3
tR(u(2))
(3.17)
3.5 Numerical Flux Function : AUSM+
The upwind numerical inviscid ux function, AUSM+ (Advection Upstream Splitting
Method +), developed by Liou in [13] was used in the present work. The method is based
on the idea of the inviscid ux consisting of two physically distinct parts: the convective
uxes and the pressure uxes. The convective uxes are associated with the advection speed
( ow speed) while the pressure uxes are associated with the acoustic speed. The convective
and pressure uxes are formulated using the eigenvalues of the ux Jacobian matrices. The
resulting AUSM+ scheme has proven to be highly robust and accurate for varying types
of uid dynamics problems. The decreased computational cost, along with the approved
accuracy, when compared to other ux-vector or ux-di erence splitting methods was the
driving force to implement the method in the present work.
The Euler equations given by
@U
@t +
@
@xF(U) +
@
@yG(U) = 0 (3.18)
consist of the state vector of conservative variables, U, along with the components of the
inviscid ux terms, F(U) and G(U). The numerical uxes through each element edge can
31
be written as combinations of convective and pressure terms
^f = nxF(U) +nyG(U) = Vn
8>
>>>
>>><
>>>
>>>>
:
u
v
h0
9>
>>>
>>>=
>>>
>>>>
;
+p
8>
>>>
>>><
>>>
>>>>
:
0
nx
ny
0
9
>>>
>>>
>=
>>>
>>>>
;
= ^fc + ^p
8
>>>
>>>
><
>>>
>>>>
:
0
nx
ny
0
9
>>>
>>>
>=
>>>
>>>>
;
(3.19)
where Vn = unx +vny is the convective velocity normal to the corresponding element inter-
faces and h0 is the total enthalpy given by h0 = e0 + p= . Examining Equation 3.19, the
numerical ux, ^f, has been decomposed to the sum of the numerical convective ux, ^fc, and
the numerical pressure ux, ^p.
The formulation begins with the computation of the average element interface speed of
sound, with the left and right side of the interface denoted by L and R respectively.
aL;R = aL +aR2 (3.20)
The interface speed of sound is found by
aface = min
a2
L
c L;
a2R
c R
(3.21)
where the variable c is computed using the corresponding left and right states, c =
max (a;j(nxu+nyu)j). The speed of sound, a, for the left and right faces is also found
by using the corresponding sides by
a =
s
2( 1)
( + 1)h0 (3.22)
The Mach number in the left and right elements are then found by
ML = Vn;La
L;R
; MR = Vn;Ra
L;R
(3.23)
32
where the value aL;R is the average of the speed of sound. The Mach number and pressure
splittings can then be found by
M (M) =
8
><
>:
1
2(M jMj; if jMj 1;
14(M 1)2 (M2 1)2; otherwise
P (M) =
8
><
>:
1
2(1 sign(M)); if jMj 1;
1
4(M 1)
2(2 M) M(M2 1)2; otherwise
(3.24)
where and are two constants with bounds such that ( 3=4 3=16) and ( 1=16
1=2). All of the test cases presented use the values ( ; ) = (3=16;1=8). Liou [13] stated that
these values show an improvement over the original AUSM method and are comparable to
the Roe ux splitting scheme. The interface values are then found by
Mface = M+ +M
pface = (P+pL) + (P pR)
(3.25)
The nal numerical ux is then computed by
^f = aface
0
BB
BB
BB
B@
1
2(Mface +jMfacej)
8
>>>>
>>>
<
>>>
>>>
>:
u
v
h0
9
>>>>
>>>
=
>>>
>>>
>;
L
+ 12(Mface +jMfacej)
8
>>>>
>>>
<
>>>
>>>
>:
u
v
h0
9
>>>>
>>>
=
>>>
>>>
>;
R
1
CC
CC
CC
CA
:::+pface
8
>>>
>>>>
<
>>>
>>>
>:
0
nx
ny
0
9
>>>
>>>>
=
>>>
>>>
>;
(3.26)
33
3.6 Arti cial Viscosity Sensors
3.6.1 Modal Decay Sensor
The following formulation follows the work presented by Klockner, Warburton, and
Hesthaven in [1]. The idea behind their work builds upon a numerical analysis of the modal
coe cients of the approximate solution in each element. Since this work focuses on the
approximate solution of the Euler equations, the nodal values of density are converted to
their modal coe cient form, f^qngNp 1n=0 , through the Vandermonde matrix relationship given
in Equation 2.24. This requires a matrix multiplication of order [Np;Np] to convert the nodal
values into the respective modal coe cients. The conversion of the coe cients is carried out
at each Runge-Kutta iteration. It is worth noting that other ow variables, such as Mach
number, entropy, or the residual of the conservation equations may be used as inputs into
the modal decay sensor. Given the di erent test cases analyzed, the density was found to
be the most reliable ow variable to base the sensor on.
The method focuses on the modal decay of the approximate solution in each element.
Modal decay is the shrinking of modal coe cient magnitudes,j^qnj, as the modal number, n,
increases and is given by
j^qnj cn s (3.27)
and taking the logarithm of the relationship yields
logj^qnj log(c) slog(n) (3.28)
34
The coe cients, s and log(c), can then be found by solving the overdetermined system of
equations of the form Ax = b where
A =
8
>>>
>>>
><
>>>
>>>>
:
1 log(1)
... ...
1 log(N)
1 log(N + 1)
9
>>>
>>>
>=
>>>
>>>>
;
; x =
8
><
>:
log(c)
s
9
>=
>;; (3.29)
The coe cients then satisfy
min
N
p 1X
n=1
jlogj^qnj (log(c) slog(n))j2
!
(3.30)
Two methods are applied to the initial modal coe cients prior to the solution of Equa-
tion 3.27. The methods, skyline pessimization and baseline decay, counteract certain oc-
currences that may be encountered by using the raw modal coe cients acquired by the
numerical approximation. The baseline decay procedure is rst applied to the initial modal
coe cients, ^qn, resulting in the new coe cients, ~qn. The skyline pessimization procedure is
then applied to ~qn with the nal resulting modal coe cients, qn, becoming the inputs into
Equation 3.27.
In the instance that the modal coe cients contain noise that is small compared to the
magnitude of the modal coe cient, the method of baseline decay is applied to the modal
coe cients. Baseline decay adjusts for this by adding a \sense of scale" by distributing
energy according to a \perfect modal decay", which is given by
j^bnj 1qP
Np 1
i=1
1
i2N
1
nN (3.31)
35
where N is the polynomial degree of the modal expansion. The new input to the skyline
pessimization simply controls coe cients that are relatively small compared to the element-
wise norm of the approximation and takes the form
j~qnj2 :=j^qnj2+kqnk2L2(Dk)j^bnj2 for n2f1;:::;Np 1g (3.32)
The approximation of the modal decay, represented by Equation 3.27, is only capable of
generating monotone modal decay approximations. Skyline pessimization is used to ensure
that the modal decay is indeed exhibiting a monotone behavior. The method is formulated
by considering two modes, n and m, where m>n. If j~qmj j~qnj, then the small coe cient
j~qnjwas likely spurious and can be replaced. The method ensures monotone decay by raising
each modal coe cient up to the largest higher-numbered modal coe cient. The resulting
procedure for skyline pessimization is then given by
qn := max
i2fmin(n;Np 2);:::;Np 1g
j~qij for n2f1;2;:::;Np 1g (3.33)
In the case where only the rst few modes are used a nal correction must be made to
counter odd-even e ects of the modal data. This is accomplished by ensuring that the last
modal coe cient is larger than the second-to-last modal coe cient.
Now that the modal decay can be approximated by Equation 3.27, the arti cial viscosity
can be computed using the estimated decay exponent, s, in an element-wise manner. The
coe cient, s, was scaled above such that s = 1 corresponds to a discontinuous solution,
s = 2 corresponds to a C0 solution, s = 3 corresponds to a C1 solution, and so forth. The
arti cial viscosity is then determined using the formulation
(s) = 0
8
>>>
><
>>>
>:
1 s2( 1;1);
1
2(1 + sin( (s 2)
2)) s2[1;3];
0 s2(3;1):
(3.34)
36
The work by Klockner[1] references Barter and Darmofal [14] as also de ning the term, 0
as
0 =
2
2t =
h
N (3.35)
where h is the local element size and N the approximating polynomial. A reference time scale
is given by t = (N=2) t, and the nal two variables to determine the maximum arti cial
viscosity are de ned as = h=N and t h=( N2). The present work deviates slightly
in the computation of the maximum viscosity. Previous work had shown a tendency for
the maximum arti cial viscosity computed by Equation 3.35 to be overly dissipative. The
maximum arti cial viscosity will be determined using the entropy residual form of the 0
presented in the next section. Using the same maximum arti cial viscosity enables a clear
comparison to be drawn between the two methods.
3.6.2 Entropy Residual Sensor
The entropy residual sensor is taken from the work of Guermond, Zinghan, and Pasquetti
[2, 15, 16]. The premise behind the entropy production residual is based on the Second
Law of Thermodynamics which states that the entropy of a system must be nondecreasing
with time. The method monitors R for negative values which would indicate a violation
of physics. In the computational approximation, this condition would be violated when a
feature of the ow is under-resolved. In the case of a shock wave, the properties of the ow
are undergoing signi cant changes across an area on the order of mean-free paths. This is
orders of magnitude smaller than the ow eld domain discretization. The resulting under-
resolution results in a negative entropy production across the shock wave. The result is a
discontinuity sensor that detects discontinuities and sharp gradients based on the physical
characteristics of the ow.
The entropy residual sensor is incorporated as a state variable along with the Euler
equation state vector for the system of conservation laws. The addition of the entropy
residual to the state vector, U, and the inviscid ux vectors, F(U) and G(U), is expressed
37
as
U =
8
>>>
>>>
>>>>
<
>>>
>>>
>>>
>:
u
v
e0
s
9
>>>
>>>
>>>>
=
>>>
>>>
>>>
>;
F(U) =
8
>>>
>>>>
>>>
<
>>>
>>>
>>>
>:
u
u2 +p
vu
( e0 +p)u
u s
9
>>>
>>>>
>>>
=
>>>
>>>
>>>>
;
G(U) =
8
>>>
>>>>
>>>
<
>>>
>>>
>>>>
:
v
uv
v2 +p
( e0 +p)v
v s
9
>>>>
>>>
>>>
=
>>>
>>>
>>>>
;
(3.36)
where the entropy, s, is given by
s = 1 1 ln
p
(3.37)
The residual of the entropy production equation is then given in two-dimensions by
R= @@t( s) + @@x(u s) + @@y(v s) 0 (3.38)
Using a nonlinear transformation of the conservative variables, the time derivative in Equa-
tion 3.38 can be expressed in terms of the spatial derivatives. The transformation begins by
using the chain rule on the time derviative such that
@( s)
@t =
@( s)
@U
@U
@t
= @( s)@V
@U
@V
1 @U
@t
= LM 1@U@t
= W@U@t
= W @@xF(U) W @@yG(U)
(3.39)
38
The variables present in Equation 3.39 are the vector of conservative variables, U, along with
the vector of primitive variables, V, given by
U =
8
>>>
>>>
><
>>>
>>>>
:
u
v
e0
9
>>>
>>>
>=
>>>
>>>>
;
; V =
8
>>>
>>>
><
>>>>
>>>
:
u
v
p
9
>>>
>>>
>=
>>>>
>>>
;
(3.40)
the transformation matrix, M, given by
M = @U@V =
2
66
66
66
64
1 0 0 0
u 0 0
v 0 0
k u v 1=( 1);
3
77
77
77
75
; (3.41)
along with the inverse
M 1 = @V@U =
2
66
66
66
64
1 0 0 0
u= 1= 0 0
v= 0 1= 0
k ( 1) u ( 1) v ( 1) ( 1);
3
77
77
77
75
(3.42)
and the resulting entropy variables [17], W, given by
W =
s+ kp 1; up ; vp ; p
(3.43)
where k is the kinetic energy given by k = 1=2(u2 + v2) . The entropy residual is then
formulated using only the spatial derivative terms by substituting the transformation from
39
Equation 3.39 into the initial Equation 3.38. This results in
R= @@x(u s) + @@y(v s) W @@xF(U) W @@yG(U) 0 (3.44)
where F(U) and G(U) are the original Euler inviscid uxes. After the entropy residual is
determined, the state vector value s is then recomputed from the updated state and the
next Runge-Kutta step is performed.
Once the residual, R, of the transformed entropy production equation is determined
using Equation 3.44, the arti cial viscosity can then be determined. The maximum arti cial
viscosity is determined by setting a maximum arti cial viscosity limit, 0, found by
0 = cmaxh N2 (3.45)
where cmax is a constant user-de ned parameter dependent on the test case. The value of 0
computed in Equation 3.45 was also used as the maximum arti cial viscosity for the modal
decay sensor. The additional factor of N in the denominator, when compared to Equation
3.35, not only seemed to bene t the modal sensor, but allowed an even comparison between
the two arti cial viscosity methods.
The entropy sensor proceeds by de ning a viscosity term, s, in each element that is
determined by using the residual values,R, that have been developed at each nodal location.
The formulation of s for each element is given by
s = cs h
2
N2 max
0; Rj sj
ref
(3.46)
where j sjref is a normalizing coe cient ensuring that s has the proper viscosity units
and cs is a user-de ned, problem speci c constant. The nal piecewise constant value of
40
arti cial viscosity in each element is then found through
reg = min ( s; 0) (3.47)
The entropy sensor, while based in the physics of the ow problem, has the unfortunate
property of containing the two user-de ned parameters, cmax and cs, in order to scale the
arti cial viscosity. This is opposed to the modal sensor that is a parameter free method
since 0 is determined through known ow variables. The parameters are only dependent
on the temporal and spatial integration methods and are independent of the timestep and
mesh size [16]. Guermond has suggested a simple method of nding the parameters by
using a coarse grid that should enable a quick determination of the values. Through his
research, Guermond suggests the parameters typically fall in the range of cmax2[0:15;0:5]
and cs2[0:1;1] while also neglecting the normalizing term j sjref.
In the present work a slightly di erent approach was taken in regards to the user-
de ned parameters. Through experiments the constant, cmax, was found to vary between 1
and 4 depending on the type of test case. The present work used cmax = 1 as a basis for
all test cases, but results were included for cmax = 4 for the shock-vortex interaction as a
comparison. It was found that the cs term could be set to one such that the scaling factor
could be incorporated into the normalizing coe cient j sjref.
3.7 Computational E ciency
The interest in higher-order methods stems from the ability to achieve higher accuracy
at the cost of increased work for each degree of freedom. If a given accuracy is required for
a solution, a higher-order method can achieve the accuracy more e ciently while using less
degrees of freedom than a lower-order method. This is due to the increased work required per
degree of freedom when using a higher-order method. A lower-order method must increase
41
the number of approximating discretized elements to achieve the same order of accuracy.
Increasing the number of elements, results in an increase in degrees of freedom.
3.7.1 Sparse-Matrix Multiplication
Sparse-matrix multiplication was used in the present work to improve computational
e ciency. The sti ness matrices, Sx and Sy, are of the order [(N + 1)2;(N + 1)2] with
(N + 1)(N + 1)2 nonzero entries. The sparse layout of the sti ness matrices are shown in
Figure 3.1. The lifting matrices, Lx and Ly, are of the order [(N + 1)2;(N + 1);2] with
(N + 1)(N + 1) nonzero entries for each of the third dimensions. The sparsity of the lifting
matrices can be seen in Figure 3.2.
The sparse matrix multiplication routine was used in the main DG scheme when com-
puting the local solution in each cell with the sti ness matrix. The routine was also employed
when the results of the numerical ux function at the element interfaces is \lifted" through
the local solution in each element. The wallclock time savings that the sparse matrix mul-
tiplication returned scaled directly with the choice of the approximating polynomial order
used in the DG method. To gauge the e ectiveness of the sparse matrix multiplication rou-
tine, the standard isentropic vortex test case, presented later, showed a 25% reduction in
wall-clock time for a 5th-order polynomial, a 34% reduction for a 7th-order polynomial, and
a 45% reduction for a 9th-order polynomial.
42
(a) Sparse X-Sti ness Matrix (b) Sparse Y-Sti ness Matrix
Figure 3.1: Sparse Sti ness Matrices for N = 9
(a) Sparse X-Lifting Matrix (b) Sparse Y-Lifting Matrix
Figure 3.2: Sparse Lifting Matrices for N = 9
3.7.2 Open-MP
The discontinuous nature of the DG method promotes a computational environment
that easily incorporates a parallelization of the computational scheme. The local solutions
computed in each element are independent of adjacent elements allowing the solutions to
43
be computed independently and thus each element may be handled by a single processor.
After the local solutions are computed, the global solution is computed with each element?s
edge data by use of a numerical ux function. The numerical ux function can then be
computed in a parallel setting because each element interface is independently solved using
the local element solutions that were previously solved. The present numerical code, being
developed in Fortran 90, utilized the Open-MP software to parallelize the code. The Open-
MP functionality was simply adapted by adding command lines into the base code.
To highlight the time-savings, the run-times are presented in Figure 3.3. The values
are given as the percentage of time normalized by the serial case (1 processor used). Using
the results, it was found that a better e ciency of running two cases using three or four
processors each became more e cient than using all eight processors on a single case run.
0 2 4 6 80.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
Number of Processors
Normalized Run Time
Figure 3.3: Decrease in Wallclock Time Due to Parallelization
44
Chapter 4
Test Cases
4.1 Isentropic Vortex and Kelvin Helmholtz Shear Layer
4.1.1 Isentropic Vortex
The isentropic vortex is a two-dimensional test case where velocity perturbations are
added into a uniform ow eld. The initial vortex is convected through the ow eld by
the mean velocity eld without any change to the vortex?s size or strength. Due to the
absence of any di usive mechanisms present in the test case, the isentropic vortex presents a
good measurement of the amount of inherent arti cial viscosity that occurs in the numerical
discretization method being used. The freestream values are determined by
u1 = M1cos( ); v1 = M1sin( ); p1 = 1= ; T1 = 1= (4.1)
where is the ow angle relative to the x-axis. The test case was performed using M1 = 0:5
and = 0:0. The above relations result in a density of 1 = 1 and an acoustic speed a1 = 1.
The vortex is formed by the addition of perturbations de ned by
u = yf; v = xf; T = 12 f2 (4.2)
The vortex potential, f is de ned by
f = B2 expA(1 r2); r2 = x2 + y2 (4.3)
45
The parameters A and B are referred to as the vortex strength and gradient, respectively,
and the vortex Mach number can be determined from those parameters by Mv = B=2 . The
parameters used here are A = 0:5 and B = 5:0. The pressure may then be computed using
the isentropic relation
p = p1
T
T1
1
= 1 ( T) 1 (4.4)
The domain is given by x;y2[ 5;5] x [ 5;5] with doubly-periodic boundary conditions in
the x and y directions. The nal time was computed such that the vortex would return to
the initial starting location after traversing the domain. This corresponded to a tfinal = 20:0.
4.1.2 Kelvin-Helmholtz Shear Layer
The Kelvin-Helmholtz shear layer is a two-dimensional test case with periodic boundary
conditions in both the spatial directions. The domain has dimensions in x,y of [0;2 ] x [0;2 ].
Initial ambient density and pressure values are = 1 and p = 1= . An instability appears in
the ow eld due to a free shear layer that is caused by a velocity gradient between two ows.
The ow domain consists of an horizontal inner jet moving through an ambient region. The
initial horizontal and vertical velocity conditions are given by
s =
8
><
>:
M tanh 3 2y2 for y>
M tanh 2y 2 otherwise
u = 12(s+ 1)M
v = M
sin
x
3
sin
2x
3
+ sin
3x
3
(4.5)
The ow parameters used in this work are = =15, = 0:1=4:9, and a jet Mach number of
M = 0:25. In order to initialize the instabilities, vertical velocity perturbations are seeded
into the initial ow eld. The perturbations are given by
v = M sin(x) (4.6)
46
The Kelvin-Helmholtz test case presents an example of a shock free ow that will still
develop numerical instabilities due to the sharp gradients that occur as a result of the shear
layer. As the solution progresses, the shear layer becomes thinner and thinner against a xed
resolution. This enables a comparison of the two sensors as to the onset and application of the
arti cial viscosity. In order to examine the e ects of each sensor over long time integrations,
a nal time of t = 300 was used. The arti cial viscosity parameter, j sjref, was found
experimentally to be 0.01.
4.2 Sod Shock Tube
The Sod shock tube is a one-dimensional test case in which a gas at two di erent states
is initially held separate from each other. At t = 0, the initially stagnant gases are allowed to
interact resulting in an unsteady ow. The ow eld exhibits a shock wave that moves into
the low-pressure state, a rarefaction wave that expands into the high-pressure state, and a
contact discontinuity that separates the two regions. The exact solution for the test case is
found iteratively using the method explained by Toro in [18]. The initial conditions for the
Sod shock tube problem, where x2[ 1;1], are
8
><
>:
= 1:0; u = 0:0; p = 1:0; if x< 0
= 0:125; u = 0:0; p = 0:1; if x> 0
(4.7)
The computation is carried out to a nal time of 0.4. The entropy sensor parameter,j sjref,
is taken as 0.095 which corresponded to the middle value in between the minimum and
maximum initial values of s.
The test case is well suited for examining unsteady compressible computational uid
dynamics solvers. The shock tube case enables a comparison of each sensor?s ability to
separate between the contact discontinuity and shock wave. The presence of the expansion
region also allows a comparison of each arti cial viscosity method and their ability to not
unduly di use the important ow features. With the resulting ow eld exhibiting many ow
47
features, each sensor is tested for the ability to separate and properly apply di usion where
needed.
4.3 Shu-Osher Planar Problem
The one-dimensional Shu-Osher planar wave test case follows the work by Shu and Osher
in [19]. The present test case exhibits a Mach 3 shock wave propagating across the domain
from left to right. The upstream and downstream boundary conditions are both xed by the
initial conditions. Small scale features originate as the shock encounters a sinusoidal density
eld. As the shock passes through the density eld, high frequency oscillations occur in the
immediate vicinity behind the shock wave until oscillations with the same frequency as the
initial sinusoidal pro le occur. The small scale features that occur behind the shock provide
insight of each arti cial viscosity sensor ability to not unduly di use these structures. The
Shu-Osher test case presents the ability of the arti cial viscosity sensors to apply di usion in
the region of the shock wave while not over di using the resulting small-scale features that
result around the shock wave.
The domain is given as x2[0;10]. The initial condition are given as
8
><
>:
= 3:857143; u = 2:629367; p = 10:333333; if x< 1
= 1 + 0:2 sin(5x); u = 0:0; p = 1:0; if x> 1
(4.8)
The simulation is carried out to a simulation time of t = 1.8. The entropy sensor parameter,
j sjref, is taken as 0.5 and corresponds to the middle value of the minimum and maximum
s initial values.
48
4.4 Shock-Vortex Interaction
The shock-vortex interaction test case follows Case G presented in the work by Inoue
and Hattori in [20]. The test case is a crude model for sound generation in compress-
ible turbulent ows. The interaction shows an acoustic wave that is composed of four
alternating compression and rarefaction regions [20]. The test case occupies the domain
x;y2 [ 10;20]x[ 15;15] with a stationary vortex located at (x;y) = (5;0) and a moving
shock starting at (x;y) = ( 5;0) and moving downstream to the right. The Mach numbers
of the shock and vortex are Ms = 1:29 and Mv = 0:39. The ow parameters used for the
vortex in the present case, which were de ned in the isentropic vortex case, are A = 0:5 and
B = 2:45044. The perturbations that de ne the vortex are added to the initial conditions
downstream of the shock wave. The ow parameters downstream of the normal shock wave
are given by
pL = pR
1 + 2 + 1(M2s 1)
L = R
1 + + 1 1 pLp
R
+ 1
1 +
pL
pR
1
uL = Vs
1 R
L
(4.9)
with the upstream conditions given by
pR = 1; R = 1 ; uR = 0 (4.10)
The shock is moving into the stagnant region with a speed
Vs = Ms
r p
R
R (4.11)
The upstream and downstream boundary conditions at x = 10;20 are xed while the
upper and lower boundaries at y = 15 are periodic. The simulation is run out until a nal
49
time of tfinal = 16:20. The entropy sensor parameter, j sjref, was determined through
experimentation and is taken as 0.5. While the choice used here was less dissipative and
allowed oscillations to form on the element interface locations it did not overly di use a
resulting double shock inside the vortex and was run at this value to show this.
50
Chapter 5
Results
5.1 Isentropic Vortex and Kelvin Helmholtz Shear Layer
5.1.1 Isentropic Vortex
The test case for the isentropic vortex was performed for polynomial approximation
orders of 3 through 9 with 100, 225, 400, 625, and 900 total number of elements. The degrees
of freedom for each set of discretized parameters are given, for quadrilateral elements, by
DoF = (N+1)2(NI NJ) and is presented in Table 5.1. Because the isentropic vortex case is
simply a convected vortex with no dissipation present, using a nal time of t = 20 convects
the vortex back to the initial position when periodic boundary conditions are employed.
Using the initial conditions, the L2-errors are calculated for the resulting ow eld density
values and are presented in Table 5.2. The bold face values are used to highlight the grid
parameters that exhibit the same magnitude errors that are presented in Table 5.2. The
computational wall-clock times are given in Table 5.3 for each run case. All of the test cases
were run using 7 processors with Open-MP.
Analyzing the bold faced value at N = 6 and N = 9, the 9th-order approximation
maintains an L2-error with the same order of magnitude as the 6th-order approximation.
The 9th-order approximation does this with 67:3% fewer degrees of freedom and 72:5% less
time than the 6th-order approximation. Looking at the bold face cases of N = 7 and N = 8,
the order of the error is still of the same magnitude, but the DG method does enjoy a
slightly lower value. Comparing these two cases, the 8th-order approximation results in a
28:8% decrease in degrees of freedom and a 16:8% decrease in runtime over the 7th-order
approximation. As can be determined from the above two analysis, using a higher-order
51
polynomial approximation will result in a decreasing amount of runtime to determine the
solution. This is achieved due to the decreasing number of needed degrees of freedom to
compute the same order of magnitude for the error values.
The exponential convergence of the DG method is seen in Figure 5.1 where the L2-
errors are plotted for the computed density eld. Figure 5.1 presents the convergence rate
for the L2-errors of the computed density eld. As the polynomial order is increased, the
DG method enjoys an exponential convergence rate for the smoothly varying solution eld.
Degrees of Freedom for Quadrilateral Elements
K N=3 N=4 N=5 N=6 N=7 N=8 N=9
100 1600 2500 3600 4900 6400 8100 10000
225 3600 5625 8100 11025 14400 18225 22500
400 6400 10000 14400 19600 25600 32400 40000
625 10000 15625 22500 30625 40000 50625 62500
900 14400 22500 32400 44100 57600 72900 90000
Table 5.1: Degrees of Freedom for Grid Discretization
L2 Error for Density Values
K N=3 N=4 N=5 N=6 N=7 N=8 N=9
100 2.32E-01 4.84E-02 1.21E-02 3.07E-03 8.64E-04 2.21E-04 4.69E-05
225 5.91E-02 1.08E-02 2.16E-03 3.38E-04 8.16E-05 1.06E-05 2.01E-06
400 2.40E-02 3.99E-03 5.59E-04 8.2E-05 1.37E-05 1.30E-06 3.18E-07
625 1.22E-02 1.77E-03 1.83E-04 2.39E-05 2.96E-06 3.65E-07 1.83E-07
900 7.60E-03 9.25E-04 7.82E-05 8.91E-06 8.44E-07 8.83E-08 8.99E-09
Table 5.2: Density L2 Error for Number of Elements, K, and Polynomial Order, N
52
Wall-Clock Times (seconds)
K N=3 N=4 N=5 N=6 N=7 N=8 N=9
100 1.007 2.470 4.114 8.419 11.867 24.064 39.592
225 2.767 8.527 13.438 28.005 42.946 82.471 128.859
400 7.781 20.520 37.131 70.809 99.169 214.139 334.289
625 14.386 34.029 70.325 143.711 207.170 419.264 678.581
900 24.90 73.87 127.36 249.07 460.86 594.30 1138.67
Table 5.3: Wall-Clock for Number of Elements, K, and Polynomial Order, N
2 4 6 8 10 1210
?10
10?8
10?6
10?4
10?2
100
Polynomial Order
L2 Error
NI,NJ = 10
NI,NJ = 15
NI,NJ = 20
NI,NJ = 25
NI,NJ = 30
Figure 5.1: L2 Density Errors for Approximating Polynomial Order
5.1.2 Kelvin-Helmholtz Shear Layer
The results presented for the Kelvin-Helmholtz test case were run using 27 elements in
the i-direction (NI=27) and 18 elements in the j-direction (NJ = 18). The modal sensor took
a total wall-clock time of 23938.839 seconds compared to 23346.527 seconds for the entropy
53
sensor. This equated to the modal sensor having a 2:5% longer wall-clock computation
time compared to the entropy sensor, or a 9.87 minute savings from a 6.65 hour run. The
modal sensor applied arti cial viscosity to 53:2% fewer total elements over the simulation run
than the entropy sensor, but did apply 32:5% more total arti cial viscosity compared to the
entropy sensor. The reasons for this can be seen through select timesteps of the simulation
procedure.
Figure 5.2 shows the energy concentration at t = 61.004. A qualitative inspection shows
that the solution utilizing the modal sensor exhibits more di usive e ects compared to the
entropy sensor. The modal sensor shows a higher e ect of di usion at the edges of the
shear layer and in the developing vortices. Contrary to what may be expected, Figure 5.3
shows the magnitude and element locations where arti cial viscosity was being applied in
the ow domain. The modal sensor applied a total amount of arti cial viscosity equaled to
0:1106358 across the ow eld while the entropy sensor applied a total value of 0:1510923.
This corresponds to a 36:57% increase of the total sum of arti cial viscosity across the
ow eld for the entropy sensor compared to the modal sensor at the current timestep.
Figure 5.4 presents the energy concentration at the nal simulation time of t = 300. The
energy concentration for the modal sensor shows a great deal of di usive errors in comparison
to the entropy sensor. Using the vortex in Figure 5.4 that resides in the region (x;y) (9;3)
as a comparison, the modal sensor, Figure 5.4a, has lost the sharp gradients surrounding the
vortex and the features have been lost into one blurred or smeared structure. The solution
computed with the entropy sensor, given in Figure 5.4b, still possesses sharp features of the
shear layer. The total arti cial viscosity applied at the nal time step for the entropy sensor
is 63:3% less compared to the modal sensor.
Figure 5.6 presents the total number of elements activated with arti cial viscosity (Fig-
ure 5.6a) and the total amount of arti cial viscosity applied throughout the ow domain
(Figure 5.6b) for selected timesteps throughout the simulation. The modal sensor activated
fewer elements per timestep with arti cial viscosity for almost the entirety of the simulation.
54
Conversely, with the exception of the beginning of the simulation corresponding to the ini-
tial vortex formulation, the modal sensor applied more total arti cial viscosity per timestep
compared to the entropy sensor. Over the course of the entire simulation, the entropy sensor
applied arti cial viscosity to 53:18% more elements, while also applying 32:46% less total
arti cial viscosity. The entropy sensor applied arti cial viscosity with a much smaller value,
orders of magnitude, in majority of the elements than the modal sensor. To emphasize this
point, returning to Figure 5.5 (t=300), the modal sensor is applying arti cial viscosity to
331 total elements as opposed to the entropy sensor which is activating all 486 elements
in the ow domain. Throughout the entire simulation the modal sensor never applied the
maximum value of arti cial viscosity to an element, whereas the entropy sensor applied the
maximum amount of arti cial viscosity, 0, to 158 elements throughout the simulation. Out
of the 158 elements activated with 0, 150 of the elements were activated between t = 20
and t = 70 which corresponds to the initial vortex formations.
Figure 5.7 presents the percentage of the maximum arti cial viscosity that was applied
throughout the simulation. Due to the presence of sharp gradients as opposed to disconti-
nuities such as shock waves, both sensors apply smaller fractions of the maximum arti cial
viscosity than the other test cases to be presented. The results further support the that the
entropy sensor is applying a smaller amount of arti cial viscosity throughout the simulation.
The entropy sensor is sensitive to the resolution of the discretized domain, with small values
of the entropy residual leading to small values of arti cial viscosity being introduced into
the approximation. It is also evident that the entropy sensor has signi cantly less elements
that are not being activated compared to the modal sensor.
The modal sensor applied arti cial viscosity throughout the ow simulation by apply-
ing a high, but not maximum, amount of arti cial viscosity in the select regions once the
spurious oscillations become large enough to warrant di usion. The entropy sensor applied
the maximum amount of arti cial viscosity in the beginning of the ow simulation as the
vortices started forming. As the solution progressed in time, the entropy sensor applied a
55
relatively small amount of di usion in the elements once any sign of oscillations occurred.
The threshold for activating arti cial viscosity proved to be much lower for the entropy
sensor when compared to the modal sensor.
(a) Modal Decay Sensor (b) Entropy Residual Sensor
Figure 5.2: Kelvin-Helmholtz Energy Concentration at t = 61.0043
5 10 15 20 25
2
4
6
8
10
12
14
16
18
i elements
j elements
2
4
6
8
10
x 10?3
(a) Modal Decay Sensor
5 10 15 20 25
2
4
6
8
10
12
14
16
18
i elements
j elements
2
4
6
8
10
x 10?3
(b) Entropy Residual Sensor
Figure 5.3: Kelvin-Helmholtz Arti cial Viscosity at t = 61.0043
56
(a) Modal Decay Sensor (b) Entropy Residual Sensor
Figure 5.4: Kelvin-Helmholtz Energy Concentration at t = 300
5 10 15 20 25
2
4
6
8
10
12
14
16
18
i elements
j elements
2
4
6
8
10
x 10?3
(a) Modal Decay Sensor
5 10 15 20 25
2
4
6
8
10
12
14
16
18
i elements
j elements
2
4
6
8
10
x 10?3
(b) Entropy Residual Sensor
Figure 5.5: Kelvin-Helmholtz Arti cial Viscosity at t = 300
0 50 100 150 200 250 3000
100
200
300
400
500
600
Timestep Sample
Number of Elements with Artificial Viscosity
Modal
Entropy
(a) Number of Elements with Arti cial Viscosity
0 50 100 150 200 250 3000
0.05
0.1
0.15
0.2
0.25
Timestep Sample
Sum of Artificial Viscosity
ModalEntropy
(b) Total Sum of Arti cial Viscosity
Figure 5.6: Kelvin-Helmholtz Stats for Arti cial Viscosity Sensors
57
0
20
40
60
80
100
Percentage of nu0
Percentage of Total Sum of Artificial Viscosity
0
0?1010?2020?3030?4040?5050?6060?7070?8080?9090?100
Modal Sensor
Entropy Sensor
Figure 5.7: Percentage Distribution of Applied Arti cial Viscosity
for Sod Shock Tube
5.2 Sod Shock Tube
The Sod shock tube test case was run using four polynomial approximation orders
of 5, 7, 9 and 11 for a ow domain consisting of 100 elements. The resulting L2-errors
for the density values are presented in Table 5.4. The modal sensor resulted in a higher
accuracy solution for the given test cases. The results indicate that as the polynomial order
increased, the percent di erence of the L2-errors between the two sensors decreased. The
explanation behind this correlation can be deduced from the results given in Figure 5.8. As
the order of the polynomial approximation increased, the number of elements activated with
arti cial viscosity and the total sum of the applied arti cial viscosity for the modal sensor
increased at a much higher rate than the entropy sensor. The modal sensor appears to be
better suited for lower polynomial orders, whereas when the polynomial order increases, the
entropy sensor appears to be the better choice. This sensor arrangement is also supported by
58
Table 5.5 which presents the percent di erence of the entropy sensor when compared to the
modal sensor. As the polynomial approximation increases, the entropy sensor increasingly
becomes the more e cient scheme to stabilize the DG approximation. The resulting time
savings of the entropy sensor are expected due to the required matrix multiplication that
the modal sensor must perform. The order of the matrix, [(N +1)2;(N +1)2], increases with
an increasing order of the polynomial approximation, N. Therefore, in actuality, it is the
modal sensor becoming less e cient as the resulting larger matrix multiplication becomes
increasingly more computationally expensive to perform.
Sensor N=5 N=7 N=9 N=11
Modal 9.658E-04 7.954E-04 7.268E-04 5.899E-04
Entropy 1.168E-03 8.699E-04 7.696E-04 6.280E-04
Table 5.4: Sod Shock Tube L2-Error for Density Values
5 7 9 110
2
4
6
8
10
12
14x 105
Polynomial Order
Number of Elements with Artificial Viscosity Applied
Modal SensorEntropy Sensor
(a) Total Number of Elements Activated
5 7 9 110
20
40
60
80
100
120
Polynomial Order
Total Amount of Artificial Viscosity Applied
Modal SensorEntropy Sensor
(b) Total Amount of Arti cial Viscosity Applied
Figure 5.8: Arti cial Viscosity Statistics per Order of the
Approximating Polynomial for Sod Shock Tube
59
N=5 N=7 N=9 N=11
-1.27% -1.93% -5.75% -9.62%
Table 5.5: Wallclock Percent Decrease of Entropy Sensor from
Modal Sensor for Sod Shock Tube
5.2.1 Case for N = 11, NI = 100
Focusing on a comparison of the test case corresponding to N = 11 and NI = 100, the
entropy sensor had a 9.62% decrease in wallclock time, a 66.07% decrease in the number
of activated elements, and a 77.76% decrease in total amount of arti cial viscosity applied
throughout the simulation. The expansion wave, contact discontinuity, and shock wave fea-
tures are magni ed in Figures 5.9, 5.10, and 5.11, respectively, so that any di usion e ects
may be highlighted. The exact solution in each gure is given by the solid black line, whereas
the multicolored lines are each elemental approximated local solution. The modal and en-
tropy sensors both exhibit di usive e ects near the onset of the expansion region (Fig.5.9).
The solution resulting from the modal sensor has lost the initial sharpness and is occurring
further upstream compared to the exact solution. The solution resulting from the entropy
sensor shows a sharper pro le at the start of the expansion and a closer approximation,
compared to the modal sensor, of the exact solution. The modal sensor applied less di usive
e ects in the region of the contact discontinuity compared to the entropy sensor (Fig. 5.10).
This implies that the entropy sensor applied a larger total amount of viscosity at the contact
discontinuity compared to the modal sensor. The results for both sensors appear almost
exact for the approximated shock wave (Fig. 5.11). The only noticeable di erence is a single
spurious nodal coe cient for the entropy sensor at the beginning of the element upstream of
the shock wave. This spurious point present in the entropy based approximation has been
smeared out by the modal sensor due to the increased arti cial viscosity.
60
?0.5?0.4?0.3?0.2?0.1 0 0.10.4
0.5
0.6
0.7
0.8
0.9
1
x?location
Density
(a) Modal Decay Sensor
?0.5?0.4?0.3?0.2?0.1 0 0.10.4
0.5
0.6
0.7
0.8
0.9
1
x?location
Density
(b) Entropy Residual Sensor
Figure 5.9: Density Pro le of Shock tube (Rarefaction Wave) at Final Time 0.4
0.25 0.3 0.35 0.4 0.45 0.50.25
0.3
0.35
0.4
0.45
x?location
Density
(a) Modal Decay Sensor
0.25 0.3 0.35 0.4 0.45 0.50.25
0.3
0.35
0.4
0.45
x?location
Density
(b) Entropy Residual Sensor
Figure 5.10: Density Pro le of Shock tube (Contact Discontinuity) at Final Time 0.4
0.65 0.7 0.750.1
0.15
0.2
0.25
x?location
Density
(a) Modal Decay Sensor
0.65 0.7 0.750.1
0.15
0.2
0.25
x?location
Density
(b) Entropy Residual Sensor
Figure 5.11: Density Pro le of Shock tube (Shock Wave) at Final Time 0.4
61
The arti cial viscosity statistics are presented in Figure 5.12 for the modal and entropy
sensors. The results show that the modal sensor activated more elements with arti cial
viscosity and also applied a substantially larger total sum of arti cial viscosity at each
sampled timestep when compared to the entropy sensor. To aid in explaining these results,
Figure 5.13 presents a graphical plot of the total amount of arti cial viscosity applied per
element at each timestep. The modal sensor applied the largest amounts of arti cial viscosity
to the expansion region and shock wave continuously, while the contact discontinuity was
di used with varying amounts of viscosity corresponding to the crossing of the element
interfaces. The entropy sensor activated elements with the highest amounts of arti cial
viscosity in a region con ned to the shock wave as it moved downstream.
102030405060708090100 120 1401540
20
40
60
80
100
120
140
Timestep Sample
Number of Elements with Artificial Viscosity
ModalEntropy
(a) Total Number of Elements Activated
102030405060708090100 120 1401540
0.002
0.004
0.006
0.008
0.01
0.012
Timestep Sample
Sum of Artificial Viscosity
ModalEntropy
(b) Total Amount of Arti cial Viscosity Applied
Figure 5.12: Stats for Arti cial Viscosity Sensors for Sod Shock tube
62
(a) Modal Decay Sensor (b) Entropy Residual Sensor
Figure 5.13: Time History of Arti cial Viscosity for Sod Shock tube Case
One selling point of the development of the entropy based sensor is in the ability to
distinguish between a shock wave and a contact discontinuity, with the entropy sensor being
the only method known to apply arti cial viscosity only to the shock wave [2, 11]. Unfor-
tunately, the results of the present work did not fully support this. Figure 5.14 presents
the elements activated with arti cial viscosity regardless of magnitude. Any element acti-
vated with arti cial viscosity is shown in red and the elements with no arti cial viscosity
are shown in blue. The modal sensor activated elements in the expansion and shock wave
regions consistently throughout the simulation while also consistently activating elements
before and after the contact discontinuity. In the region of the expansion, the modal sensor
decreased the amount arti cial viscosity as the simulation was marching through time. As
the simulation progressed the modal sensor applied a larger amount of arti cial viscosity
over a longer simulation time compared to the entropy sensor for the region encompassing
the expansion. The entropy sensor shows intermittent activation in the expansion region
at the initial timesteps. Once the expansion was di used to the extent that the numerical
method could properly resolve the ow feature, the entropy sensor stopped activating ele-
ments. The region of elements starting at the contact discontinuity and extending through
the shock wave were activated consistently up to a timestep of approximately 40. After
the timestep of 40, the regions of activated elements began to separate for the two separate
63
ow features. The activation of elements around the contact discontinuity exhibits a fairly
consistent pattern while the elements surrounding the shock wave were activated constantly.
Element Number
time?step
20 40 60 80 100
20
40
60
80
100
120
140
(a) Modal Decay Sensor
Element Number
time?step
20 40 60 80 100
20
40
60
80
100
120
140
(b) Entropy Residual Sensor
Figure 5.14: Time History of Activated Elements for Sod Shock tube Case
Figure 5.15 presents the percentage occurrences of the total amount of arti cial viscosity
applied throughout the sampled timesteps. The modal sensor activates more elements overall
than the entropy sensor. The modal sensor shows that 8.5% of the total number of activated
elements where withing 0 to 10% of 0, whereas the entropy sensor activated 7.4% of the total
number of elements throughout the simulation where within the same range. The modal
sensor also exhibits a 8.8% occurence of elements within 90% to 100% of the maximum
arti cial viscosity compared to 2.8% for the entropy sensor.
64
0
20
40
60
80
100
Percentage of nu0
Percentage of Total Sum of Artificial Viscosity
0
0?1010?2020?3030?4040?5050?6060?7070?8080?9090?100
Modal Sensor
Entropy Sensor
Figure 5.15: Percentage Distribution of Applied Arti cial Viscosity
for Kelvin-Helmholtz Case
5.3 Shu-Osher Planar Problem
The results of the Shu-Osher test case presented here utilized a 9th-order polynomial ap-
proximation and the ow domain was discretized using 200 elements in the i-direction. Over
the entire simulation, the entropy sensor took 2.4% less wallclock time, activated 8.68% more
elements with arti cial viscosity, and applied 67.7% less total arti cial viscosity compared
to the modal sensor.
The approximated density ow eld values for both sensors at the nal simulation time
of t = 1.8 are given in Figure 5.16. Moving from upstream to downstream, the nal solution
consists of a constant density post-shock region, a post-shock low-frequency region, a post-
shock high-frequency region, the normal shock, and the pre-shock sinusoidal region. Upon
a visual inspection, the resulting ow elds do not exhibit any noticeable di erences. Both
sensors successfully di use all pre- and post-shock spurious activity. The only resulting
65
di erence requires zooming into the post-shock low-frequency/high-frequency interface. The
interface point approximated using the modal sensor exhibits a slightly more di usive e ect
than the entropy based approximation. The resulting percent di erence between the two
density values at the interface is 0.9285%, which shows a very similar approximation.
?5 0 50
1
2
3
4
5
x?Location
Density
(a) Modal Decay Sensor
?5 0 50
1
2
3
4
5
x?Location
Density
(b) Entropy Residual Sensor
Figure 5.16: Density Pro le of Shu-Osher Case at Final Time 1.8
The total number of elements activated and the total amount of applied arti cial viscos-
ity of the modal and entropy sensors at selected timesteps can be seen in Figure 5.17. The
modal sensor stays reasonably consistent in the number of elements being activated with
arti cial viscosity, whereas the entropy sensor increases the number of activated elements
throughout the simulation. Despite the increasing number of activated elements, the en-
tropy sensor applies at least half the amount of arti cial viscosity throughout the simulation
compared to the modal sensor.
66
204060801001201401601802002202402680
20
40
60
80
100
120
140
160
Timestep Sample
Number of Elements with Artificial Viscosity
ModalEntropy
(a) Total Number of Elements Activated
204060801001201401601802002202402680
0.02
0.04
0.06
0.08
0.1
0.12
0.14
0.16
0.18
Timestep
Sum of Artificial Viscosity
Modal
Entropy
(b) Total Amount of Arti cial Viscosity Applied
Figure 5.17: Stats for Arti cial Viscosity Sensors for Shu-Osher Test Case
The timestep history of the total amount of arti cial viscosity applied in each element
is presented in Figure 5.18. The modal sensor activated elements at the shock wave and
applied the maximum arti cial viscosity in a 10-15 element region downstream of the shock
wave. Conversely, the entropy sensor only activated the maximum arti cial viscosity in a
3-8 element region downstream of the shock wave. The variation in these two number sets
corresponds to the shock wave crossing the element boundaries, with the higher number
occurring as the shock crosses the element interfaces. The modal sensor also activated at the
interface of the resulting ow where the post-shock high-frequency solution meets the post-
shock low-frequency solution. The entropy sensor did not activate with a signi cant amount
of viscosity for this same region. At the location where the smooth upstream density eld
was meeting the post-shock result, both sensors began activating around an approximate
timestep of 150.
67
(a) Modal Decay Sensor (b) Entropy Residual Sensor
Figure 5.18: Time History of Arti cial Viscosity for Shu-Osher Case
To highlight the fact that the entropy sensor is activating in 8.68% more elements,
Figure 5.19 shows the location of activated elements regardless of the magnitude of arti cial
viscosity applied. It becomes apparent that the entropy sensor is applying arti cial viscosity
in regions that were not initially evident, whereas the modal sensor still exhibits the same
pro les. The modal sensor applied arti cial viscosity at and downstream of the shock wave
location, whereas the entropy sensor activates elements at the shock wave and continues
activating upstream across the high-frequency post-shock solution as the timestep increases.
Both sensors exhibit arti cial viscosity activation at the constant density/initial-frequency
post-shock interface as well as the initial-frequency post-shock/high-frequency post-shock
interface.
68
Element Number
time?step
50 100 150 200
50
100
150
200
250
(a) Modal Decay Sensor
Element Number
time?step
50 100 150 200
50
100
150
200
250
(b) Entropy Residual Sensor
Figure 5.19: Time History of Elements Activated with
Arti cial Viscosity for the Shu-Osher Case
Figure 5.20 presents the percentage distribution of the range of arti cial viscosity mag-
nitudes used throughout the entire simulation. Both sensors indicate a similar approximately
85% of the total number of elements where not activated with arti cial viscosity for the sam-
pled timesteps. The entropy sensor does show that 11.15% of the elements were activated
within 0 to 10% of 0 compared to 3.2% for the modal sensor. It is again apparent that
the modal sensor activates more elements within 90% to 100% of 0 as was discovered in
previous results.
69
0
20
40
60
80
100
Percentage of nu0
Percentage of Total Sum of Artificial Viscosity
0
0?1010?2020?3030?4040?5050?6060?7070?8080?9090?100
Modal Sensor
Entropy Sensor
Figure 5.20: Percentage Distribution of Applied Arti cial Viscosity
for Shu-Osher Case
The results have shown that although a considerable di erence existed in the way the
two sensors applied arti cial viscosity to the test case, the resulting ow elds exhibited the
same solution. The entropy sensor did show a time-savings over the modal sensor while also
applying a smaller sum of the total amount of arti cial viscosity applied.
5.4 Shock-Vortex Interaction
The results of the shock-vortex interaction presented below utilized a 9th-order polyno-
mial approximation inside each element. The ow domain was discretized using 60 elements
in the i-direction and 45 elements in the j-direction. The arti cial viscosity sensors were
tested using the maximum arti cial viscosity parameter of cmax = 1 and 4. The increased
value of 0 allowed comparisons to be made for the resulting behaviors of modal and entropy
70
sensors. Results were also included for the entropy sensor parameter cs = 0.1 in order to
show the e ects of the choice in the parameter.
The arti cial viscosity statistics are presented for cmax = 1 in Figure 5.21 and cmax =
4 in Figure 5.22. The results show the entropy sensor is activating elements increasingly
through the simulation for both cmax values, while the modal sensor uses roughly the same
amount of elements at each timestep for each cmax case. Focusing on Figure 5.21a and
Figure 5.22a, a signi cant di erence is noticed using the two cmax values for the modal
sensor while the entropy sensor stays more consistent. The modal sensor activated 66.11%
more elements while the entropy sensor activated 3.68% more elements for using cmax = 4.
Comparing Figures 5.21b and 5.22b for the total arti cial applied at each timestep results in a
310.74% increase for the modal sensor and only a 17.72% increase for the entropy sensor when
using cmax = 4. The percent increases are for the average values at each timestep because
an evaluation over the entire simulation run would be skewed by the increased number of
timesteps due to a larger 0 value. The modal sensor showed a substantial increase in both
the total number of elements activated and the total amount of arti cial viscosity compared
to entropy sensor when cmax is increased from 1 to 4.
2025 31364145505560657075 810
500
1000
1500
2000
2500
Timestep Sample
Number of Elements with Artificial Viscosity
ModalEntropy
(a) Total Number of Elements Activated
2025 31364145505560657075 810
0.5
1
1.5
2
2.5
3
Timestep
Sum of Artificial Viscosity
Modal
Entropy
(b) Total Amount of Arti cial Viscosity Applied
Figure 5.21: Stats for Arti cial Viscosity Sensors for
Shock-Vortex Test Case, cmax = 1
71
2025 31364145505560657075 810
500
1000
1500
2000
2500
Timestep Sample
Number of Elements with Artificial Viscosity
ModalEntropy
(a) Total Number of Elements Activated
2025 31364145505560657075 810
2
4
6
8
10
12
Timestep Sample
Sum of Artificial Viscosity
ModalEntropy
(b) Total Amount of Arti cial Viscosity Applied
Figure 5.22: Stats for Arti cial Viscosity Sensors for
Shock-Vortex Test Case, cmax = 4
Expanding upon the results of the previous section, Figure 5.23 presents the percentage
distribution of the range of arti cial viscosity magnitudes that each sensor used throughout
the selected timesteps. The modal sensor exhibits an increase in the number of elements
not activated with arti cial viscosity compared to the entropy sensor as was demonstrated
in previous the previous results. The entropy sensor shows a signi cant increase in the
percentage of elements (49.1%) that are activated with 0 to 10% of 0 compared to the
modal sensor (5.75%). The percentage of elements activated without viscosity or with the 0
to 10% range constitutes 96.45% of the total number of elements sampled, with the remainder
of the ranges being relatively close. The modal sensor activates 5.4% of the elements with the
upper range of arti cial viscosity compared to 0.56% for the entropy sensor in the same range.
The entropy sensor tends to activate a large percentage of elements with small amounts of
viscosity whenever sharp gradients are present. This trend was also evident for the Kelvin-
Helmholtz test case. For test cases with sharp gradients, the entropy sensor presents a better
capability of handling the resulting ow features.
72
0
20
40
60
80
100
Percentage of nu0
Percentage of Total Sum of Artificial Viscosity
0
0?1010?2020?3030?4040?5050?6060?7070?8080?9090?100
Modal Sensor
Entropy Sensor
Figure 5.23: Percentage Distribution of Applied Arti cial Viscosity
for Shock-Vortex Case
The density pro le at the end of the simulation are presented in Figure 5.24 for the
modal sensor and Figure 5.25 for the entropy sensor. One ow feature worth pointing out is
the wiggle present in each density pro le. This feature is the rst vertical line present in each
of the gures. The resulting wiggle is a numerical remnant of starting each simulation run
using a discontinuous solution for the shock wave. This may be avoided by rst solving for
a steady state solution of the shock wave and then introducing the vortex and allowing the
shock wave to convect downstream. Figure 5.24 shows the density eld computed using the
modal sensor for both cmax values. Severe dissipation is present for the case of cmax = 4 as
can be seen by the increase of the normal shock width and the dissipation of the double shock
inside the vortex. Using cmax of 4 results in a 161.87% increase in wallclock time compared
to using cmax = 1. The density eld is given Figure 5.25 and shows the comparison using
the two cmax values for the entropy sensor. The cmax = 4 results show a reduction in the
number of spurious oscillations when compared to the cmax = 1 test case without noticeable
73
added e ects of dissipation. Comparing the cmax = 4 to the cmax = 1 case for the entropy
sensor resulted in a wallclock time increase of 175.19%. Compared to the modal sensor, the
entropy sensor showed only a slightly higher increase in wallclock time, but a signi cantly
lower total number of activated elements and total arti cial viscosity applied when the cmax
value was increased.
(a) cmax = 1 (b) cmax = 4
Figure 5.24: Density Field for Modal Decay Sensor at Final Time = 16.2
(a) cmax = 1 (b) cmax = 4
Figure 5.25: Density Field for Entropy Residual Sensor at Final Time = 16.2
Table 5.6 gives the values of the percent di erence of the entropy sensor compared to
the modal sensor for each value of cmax. The percent di erences given are the di erence in
wallclock computation time, the di erence in the total number of activated elements over
74
the course of the simulation, and the total sum of arti cial viscosity applied throughout the
simulation.
cmax Wallclock Activated Elements Sum of A.V.
1 -3.44% 159.58% -79.33%
4 1.47% 55.38% -94.77%
Table 5.6: Percent Di erence of Entropy Sensor Compared to the Modal Sensor
The arti cial viscosity values at the nal simulation time are given in Figure 5.26 for the
modal sensor and Figure 5.27 for the entropy sensor. When cmax was increased to 4, both
methods show a substantial decrease in activating the elements with the maximum amount of
arti cial viscosity. One important aspect of the results is the activation of arti cial viscosity
that occurs in the vortex for the modal sensor. The vortex, while undergoing a disturbance
caused by the shock interaction, still maintains a smooth pro le. The modal sensor applies
a noticeable amount of arti cial viscosity to the vortex whereas the entropy sensor does not
apply a noticeable amount of arti cial viscosity observed in the gure.
10 20 30 40 50 60
5
10
15
20
25
30
35
40
45
i elements
j elements
0
2
4
6
8
x 10?3
(a) cmax = 1
10 20 30 40 50 60
5
10
15
20
25
30
35
40
45
i elements
j elements
0
0.005
0.01
0.015
0.02
0.025
0.03
0.035
(b) cmax = 4
Figure 5.26: Arti cial Viscosity for Modal Decay Sensor at Final Time = 16.2
75
10 20 30 40 50 60
5
10
15
20
25
30
35
40
45
i elements
j elements
0
2
4
6
8
x 10?3
(a) cmax = 1
10 20 30 40 50 60
5
10
15
20
25
30
35
40
45
i elements
j elements
0
0.005
0.01
0.015
0.02
0.025
0.03
0.035
(b) cmax = 4
Figure 5.27: Arti cial Viscosity for Entropy Residual Sensor at Final Time = 16.2
Figures 5.28 and 5.29 give the plots indicating which elements are activated with ar-
ti cial viscosity, regardless of magnitude, for the modal and entropy sensors, respectively.
Comparing these gures to Figure 5.26 for the modal sensor and Figure 5.27, the modal
sensor is applying arti cial viscosity mainly in large magnitudes and is applying viscosity
downstream of the normal shock. The entropy sensor is activating a small magnitude of ar-
ti cial viscosity focusing in the region upstream of the normal shock with large magnitudes
of arti cial viscosity activating in elements containing large gradients of the sensing variable.
10 20 30 40 50 60
5
10
15
20
25
30
35
40
45
i elements
j elements
(a) cmax = 1
10 20 30 40 50 60
5
10
15
20
25
30
35
40
45
i elements
j elements
(b) cmax = 4
Figure 5.28: Elements Activated with Arti cial Viscosity for
Modal Decay Sensor at Final Time = 16.2
76
10 20 30 40 50 60
5
10
15
20
25
30
35
40
45
i elements
j elements
(a) cmax = 1
10 20 30 40 50 60
5
10
15
20
25
30
35
40
45
i elements
j elements
(b) cmax = 4
Figure 5.29: Elements Activated with Arti cial Viscosity for
Entropy Residual Sensor at Final Time = 16.2
The entropy sensor has the drawback of being a method that requires a user de ned
parameter in order to scale the arti cial viscosity for a given ow simulation. One advantage
that does come from the selection of the parameter is the ability to ne tune the amount
of arti cial viscosity that is applied to the ow eld. The entropy sensor was found to
continuously activate more elements with arti cial viscosity when compared to the modal
sensor. Whereas both sensors typically applied the maximum arti cial viscosity when a
strong shock was encountered, the ability to alter the j sjref value will adjust the small
and intermediate values of arti cial viscosity that occur with the entropy sensor. Figure
5.30 presents the nal density eld values for the entropy sensor with j sjref = 0.5 and
0.1. The entropy sensor is activating approximately the same amount of elements, -0.35%
di erence from the cs = 0.5 case, but is applying 69.18% more total arti cial viscosity. The
increase in viscosity removes the spurious oscillations at the cost of di using the ow eld.
Using cs = 0.1 results in the nal density eld being subjected to more di usive errors while
reducing the dispersive errors. The increased di usion is especially noticeable in both the
normal shock and the double shock wave emanating from the vortex. The resulting pro le
of both ow features have been widened, or smeared, over the course of the simulation.
77
The ability to adjust the scaling of the arti cial viscosity allows the user to experimen-
tally determine the amount of arti cial viscosity that is applied throughout the approxima-
tion. While having to determine the cs parameter reduces the robustness of the entropy
sensor, it does allow a way to ensure that the minimum amount of arti cial viscosity is being
applied throughout the ow eld and physically correct features are not being overly di used.
(a) cs = 0.5 (b) cs = 0.1
Figure 5.30: Density Field for Entropy Residual Sensor at Final Time = 16.2
78
Chapter 6
Conclusions and Future Work
6.1 Conclusions
The present work focused on the application of arti cial viscosity as a way to stabilize
the high-order approximation of the Euler equations by the discontinuous Galerkin method.
In order to apply arti cial viscosity in the elements where the Gibb?s phenomenon are oc-
curring, as opposed to a constant arti cial viscosity across the entire domain, two methods
of detecting where the spurious oscillations where occurring where analyzed.
The modal sensor, based on the modal decay of the approximate solution, resulted in
stable solutions for all the test cases presented by analyzing the numerical behaviour of the
modal coe cients at each timestep. Although the modal sensor proved to be overly di usive
for certain test cases, the sensor is considered robust since no user parameters are needed
in order to scale the arti cial viscosity. The modal sensor was more e cient for lower-order
polynomial approximations inside each element due to the required computational expense
involved in acquiring the modal coe cients. The modal sensor also was more e cient in
the arti cial viscosity applied per element thus resulting in a large decrease in the amount
of elements that are activated for each timestep. This translates into saved computational
expense due to the calculation of the second-order derivatives in the elements where arti cial
viscosity is being applied.
The entropy sensor, which was based on the entropy production residual, resulted in a
method that is based on the physics of the resulting ow eld. The entropy sensor successfully
di used each test case so that the resulting DG scheme remained stable. The entropy
sensor become computationally more e cient as the approximating polynomial increased.
Although the entropy sensor requires the user input of the parameter,j sjref, the amount
79
of applied arti cial viscosity may be tailored to more accurately apply the minimum amount
of di usion in order that the scheme remains stable. The entropy sensor applied less total
arti cial viscosity during the entire simulation as compared to the modal sensor. This came
at the cost of applying a smaller magnitude of arti cial viscosity to more elements than
the modal sensor. The larger number of activated elements results in more computational
expense being required at each timestep.
Both arti cial viscosity sensors tested had bene ts as well as drawbacks when comparing
the two methods. The modal sensor is ideally suited as a plug-and-chug sensor that requires
no user input but still results in a stable numerical approximation. Unfortunately, this
comes at the cost of an over di usive behavior for certain ow conditions. The modal sensor
also exhibited a higher threshold before activating arti cial viscosity in the approximation
as encountered in the isentropic vortex test case. This is bene cial if a quick turnaround
is required and a coarse mesh is to be used. If a low-order polynomial is also used, the
computational cost required by the modal sensor may also be o set. Because the goal of
the DG method is to use higher-order polynomial approximations in each element in order
to achieve higher accuracy, the entropy sensor is favored due to the reduced computational
cost for higher-order polynomial approximations. This comes at the cost of the user-de ned
parameter,j sjref, that must be determined for each resulting ow eld. The choice of the
parameter seems to be non-trivial, as a starting value may be estimated by inspection of
the initial ow eld values, s, and then re ned using a course mesh. The entropy sensor
also exhibited a less di usive nature when the approximated solution resulted in small-scale
features as was encountered in the shock-vortex test case. Another ow feature that proved
advantageous for the entropy sensor was sharp gradients opposed to discontinuities. This was
evident in the Kelvin-Helmholtz test case where the entropy sensor signi cantly exhibited
less di usive e ects compared to the modal sensor. The entropy sensor proved to be more
accurate for the isentropic vortex case and also for the expansion and shock wave approxi-
mations of the shock tube case, whereas the modal sensor more accurately approximated the
80
contact discontinuity for the shock tube test case. Considering the di erence in the accuracy
at the contact discontinuity, and the fact that the di erence decreases with grid re nement,
the entropy sensor is viewed as the advantageous method.
The two sensors presented have both demonstrated the ability to successfully stabilize
the numerical approximation of the Euler equations when the DG method is employed.
Although the requirement of a user-de ned parameter of the entropy sensor leads to an
additional unknown, the advantages presented outweigh the complication of determining
the unknown parameter. In addition, the determination of j sjref enables the minimum
amount of arti cial viscosity to be exploited such that physically correct ow features are
not unduly di used.
6.2 Future Work
The two sensors presented have been compared as close to the original authors methods
as allowable. It has been observed that other steps may be taken to improve the results
obtained from each method. Because only the elements containing arti cial viscosity undergo
the determination of a second-order di erential to apply the arti cial viscosity, it is bene cial
to activate only the necessary elements in order to maintain stability. The entropy sensor
applied a small amount of arti cial viscosity, compared to the maximum value in the entire
domain, to a large number of elements. Throughout the research, the reason behind these
small magnitudes seemed to be corresponding with the grid resolution of the test case. As
the number of elements increased, the number of elements with small amounts of arti cial
viscosity would decrease due to the higher resolution. Instead of having to apply a ne mesh
in order to achieve this, it may be bene cial to lter out any arti cial viscosity occurring
below a certain tolerance. It has been observed that ltering o any value of arti cial
viscosity approximately below 10% of 0 would signi cantly reduce the number of elements
activated. An investigation of how to determine what the given tolerance should be and the
e ects that this would have on the stability of the approximated solution should be carried
81
out. The reduction of the activated elements would decrease the associated computational
expense resulting in a decrease in wallclock time.
Extending the sensor comparison to 3-D would be bene cial in comparing the wallclock
times for each method. The modal sensor, having to perform an [(N + 1)2;(N + 1)2] matrix
multiplication for a quadrilateral 2D case, would have to perform an [(N + 1)3;(N + 1)3]
matrix multiplication for a cubic 3D case. Considering that each sensor is being computed
inside each Runge-Kutta step, the increased cost of computing the modal sensor should
prove to be increasingly more expensive than the entropy sensor. It is also pointed out that
these problems were run for relatively simple test cases. The addition of the third dimension
would result in a signi cant number of elements being added to the ow eld and allow a
better comparison to the bene ts of the entropy sensor. Both sensors can be extended to
3-D with non-trivial changes to the methods already employed.
The present work computed the needed arti cial viscosity at each Runge-Kutta step.
It has been observed through the research that the methods were stable if the sensor was,
for instance, utilized at only two of the three Runge-Kutta steps. The modal sensor may
also be computed at the beginning of each Runge-Kutta step with the same value used for
the entire time integration step. An analysis of the overall e ect of these procedures on the
time savings, resulting accuracy, and needed stabilization would provide additional insight.
It appears the modal sensor would bene t from this type of analysis due to the reduction in
the number of times the matrix multiplication would need to be carried out.
The arti cial viscosity applied in this work was piecewise discontinuous across the dis-
cretized domain. Barter and Darmofal [14] have shown that the application of a non-smooth
arti cial viscosity can itself cause oscillations in the solution approximation downstream.
Upgrading the arti cial viscosity such that a smooth transition occurs throughout the entire
ow domain would bene t the overall goal of reducing unwanted oscillations in the ow eld.
Whether this is accomplished through the use of a simple linear interpolation between each
element or through a pde-based model [14], should be determined by the improvement in
82
the stabilization and accuracy versus the incurred computational expense required by each
method. In addition the application of the arti cial viscosity may be applied using the form
presented by Zinghan and Guermond in Equation 3.10. This form may be bene cial at re-
ducing oscillations that occur upstream due to the parabolic regularization form of viscosity
terms that were employed here.
The modal sensor, being a parameter free method, does not allow the user to adjust
the amount of applied arti cial viscosity for di erent ow eld conditions. Both methods
exhibited di erent approaches to stabilizing the numerical approximation: the modal sensor
using a large amount of arti cial viscosity in a small number elements and the entropy sensor
applying a smaller amount of arti cial viscosity in a larger number of elements. It has been
observed that a simple scaling factor to the determination of the maximum arti cial viscosity,
0, could successfully limit the over di usive nature of the method. Scaling 0 by multiplying
the original value by a factor of 1=100 for the Kelvin-Helmholtz test case produced a result
that rivaled the entropy sensor. With the introduction of a scaling parameter, the increased
computational expense may prove to be o set by the reduction in the number of activated
elements for 2-D problems.
83
Bibliography
[1] Kl ockner, A., Warburton, T., and Hesthaven, J., \Viscous Shock Capturing in a Time-
Explicit Discontinuous Galerkin Method," Mathematical Modelling of Natural Phenom-
ena, 2010, pp. 1{42.
[2] Zingan, V., Guermond, J. L., Morel, J., and Popov, B., \Implementation of the Entropy
Viscosity Method with the Discontinuous Galerkin Method," Computational Methods
Applied Mechanical Engineering, Vol. 253, 2013, pp. 479{490.
[3] Shelton, A. B., A Multi-Resolution Discontinuous Galerkin Method for Unsteady Com-
pressible Flows, Ph.D. thesis, Georgia Institute of Technology, 2012.
[4] Hesthaven, J. and Warburton, T., Nodal Discontinuous Galerkin Methods - Algorithms,
Analysis, and Applications, Springer, New York, NY, 1st ed., 2008.
[5] Tannehill, J. C., Anderson, D. A., and Pletcher, R. H., Computational Fluid Dynamics
and Heat Transfer, Taylor and Francis, Philadelphia, PA, 2nd ed., 1997.
[6] Persson, P. O. and Peraire, J., \Sub-Cell Shock Capturing for Discontinuous Galerkin
Method," In Proc. of the 44th AIAA Aerospace Sciences Meeting and Exhibit, Vol. 112,
2006.
[7] Atkins, H. L. and Pampell, A., \Robust and Accurate Shock Capturing Method for
High-Order Discontinuous Galerkin Methods," AIAA, 2011.
[8] Huerta, A., Peraire, J., and Casoni, E., \A Simple Shock-Capturing Technique for High-
Order Discontinuous Galerkin Methods," International Journal of Numerical Methods
in Fluids, 2011, pp. 1614{1632.
[9] Gautschi, W., Orthogonal Polynomials Computation and Approximation, Oxford Uni-
versity Press, Oxford, NY, 1st ed., 2004.
[10] Hesthaven, J. and Warburton, T., \From Electrostatics to Almost Optimal Nodal
Sets for polynomial Interpolation in a Simplex," SIAM Journal of Numerical Analy-
sis, Vol. 35, No. 2, 1998, pp. 655{676.
[11] Zingan, V., Discontinuous Galerkin Finite Element Method for the Nonlinear Hyper-
bolic Problems with Entropy-Based Arti cial Viscosity Stabilization, Ph.D. thesis, Texas
A&M University, 2012.
84
[12] Gottlieb, S. and Shu, C.-W., \Total Variation Diminishing Runge-Kutta Schemes,"
Mathematics of Computation of the American Mathematical Society, Vol. 67, No. 221,
1989, pp. 73{84.
[13] Liou, M.-S., \A sequel to AUSM, Part II: AUSM+-up for all speeds," Journal of Com-
putational Physics, Vol. 214, No. 1, 2006, pp. 137{170.
[14] Barter, G. E. and Darmofal, D. L., \Shock Capturing with PDE-based Arti cial Vis-
cosity for DGFEM: Part I. Formulation," Journal of Computational Physics, Vol. 229,
No. 5, 2010, pp. 1810{1827.
[15] Guermond, J. L., Pasquetti, R., and Popov, B., \Entropy Viscosity Method for Nonlin-
ear Conservation Laws," Journal of Computational Physics, 2011.
[16] Guermond, J.-L., Pasquetti, R., and Popov, B., \From Suitable Weak Solutions to
Entropy Viscosity," Journal of Scienti c Computing, Vol. 49, No. 1, 2011, pp. 35{50.
[17] Fidkowski, K. and Roe, P., \An Entropy Adjoint Approach to Mesh Re nement," SIAM
Journal of Scienti c Computing, 2010.
[18] Toro, E. F., Riemann Solvers and Numerical Methods for Fluid Dynamics: A Practical
Introduction, Springer-Verlag Berlin Heidelberg, 3rd ed., 2009.
[19] Shu, C.-W. and Osher, S., \E cient Implementation of Essentially Non-Oscillatory
Shock-Capturing Schemes, II," Journal of Computational Physics, Vol. 83, No. 1, 1989,
pp. 32{78.
[20] Inoue, O. and Hattori, Y., \Sound Generation by Shock-Vortex Interactions," Journal
of Fluid Mechanics, Vol. 380, No. 3, 1999, pp. 81{116.
85