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