. ”v
1 if. . .
.t e6 4. \ .
t Li. , Hafiz. $13.5.
fir ratio; 3.23“

 

 

 

 

\ i211 \- )ll.v :1.

33. sgéfiwgfi . . gar 33%

.1 vol. .,

 

I
.u
13;.

N ’t‘ MICHIGAIJ‘sBRAmES
9 TATE UNIVERSITY
Com” 9 3 Q5 <9 v EAST LANSING, MICH 48824-1048

This is to certify that the
dissertation entitled

HIGH ORDER NUMERICAL METHODS FOR INVISCID
AND VISCOUS FLOWS ON UNSTRUCTURED GRIDS

presented by

YUZHI SUN

has been accepted towards fulfillment
of the requirements for the

PhD degree in Mechanical Enmeering

 

 

£04m pmficfi

 

Major Professor’s Signature

r2.[i7IZoo1-b.

 

Date

MSU is an Affinnative Action/Equal Opportunity Institution

 

PLACE IN RETURN BOX to remove this checkout from your record.
TO AVOID FINES return on or before date due.
MAY BE RECALLED with earlier due date if requested.

 

DATE DUE DATE DUE DATE DUE

 

 

 

 

 

 

 

 

 

 

 

 

 

 

6/01 c2/CIRC/DateDuep65‘p15

HIGH ORDER NUMERICAL METHODS FOR INVISCID
AND VISCOUS FLOWS ON UNSTRUCTURED GRIDS

By

Yuzhi Sun

A DISSERTATION

Submitted to
Michigan State University
in partial fulfillment of the requirements
for the degree of

DOCTOR OF PHILOSOPHY
Department of Mechanical Engineering

2004

ABSTRACT

HIGH ORDER NUMERICAL METHODS FOR INVISCID AND VISCOUS
FLOWS ON UNSTRUCTURED GRIDS

By

Yuzhi Sun

The spectral volume (SV) method is a newly developed high-order, conservative, and
efficient finite volume method for hyperbolic conservation laws on unstructured grids. It
has been successfully demonstrated for scalar conservation laws and multi-dimensional
Euler equations. In this study, the SV method is compared with another high-order
method for hyperbolic conservation laws capable of handling unstructured grids named
the discontinuous Galerkin (DG) method. Their overall performance in terms of the
efficiency, accuracy and memory requirement is evaluated using the scalar conservation
laws and the two-dimensional Euler equations. To measure their accuracy, problems with
analytical solutions are used. Both methods are also used to solve problems with strong
discontinuities to test their ability in discontinuity capturing. Both the DO and SV
methods are capable of achieving the formal order of accuracy while the DG method has
a lower error magnitude and takes more memory. They are also similar in efficiency. The
SV method appears to have a higher resolution for discontinuities because the data

limiting can be done at the sub-element level.

The SV method is also successfully extend to the Navier-Stokes equations. First, the SV

method is extended to and tested for the diffusion equation. In this study, three

different formulations named Naive SV, Local SV and Penalty SV for the diffusion
equation are presented. The Naive SV formulation yields an inconsistent and unstable
scheme, while the other two formulations are consistent, convergent and stable. A Fourier
type analysis is performed for all the formulations, and the analysis agrees well with the
numerical results. Second, the Local SV method is chosen to be extended to solve the
Navier-Stokes equations since it gives the optimum accuracy in solving the diffusion
equation. The formulation of the Local SV method for the two-dimensional compressible
Navier-Stokes equations is described. Accuracy studies are performed on the scalar
convection-diffusion and the Navier-Stokes equations using problems with analytical
solutions. It is shown that the designed order of accuracy is achieved for 15‘, 22nd and 3rd
order reconstructions. The solver is then used to solve other viscous laminar flow

problems to demonstrate its capability.

Copyright by
YUZHI SUN

2004

ACKNOWLEDGEMENTS

I wish to express my sincere gratitude to my advisor, Professor 2.]. Wang, for his
invaluable guidance throughout this study. He has taught me many things and this work

would not have been possible without his patience and encouragement.

I would also like to thank Professors Chichia Chiu, Farhad Jaben', Tom Shih, and Mei

Zhuan g for serving on my committee and for reviewing the work and the dissertation.

I would like to give special thanks to my husband Dongsheng Wu, and my daughter

Emily Wu, who have constantly supported and encouraged me throughout the study.

TABLE OF CONTENTS

LIST OF TABLES .................................................................................. viii

LIST OF FIGURES ................................................................................ xi

NOMENCLATURE .............................................................................. xiii
CHAPTER 1

INTRODUCTION ................................................................................. 1

1.1 Background of Computational Fluid Dynamics (CFD) .......................... l

1.2 Numerical Methods in CFD ......................................................... 2

1.3 Motivation and Objectives of This Study .......................................... 5

1.4 Outline of the Dissertation ............................................................ 6
CHAPTER 2

FRANEWORK OF DG AND SV METHODS ............................................... 8

2.1 DC Method ............................................................................ 8

2.1.1 Space Discretization ......................................................... 9

2.1.2 Grids and Data Reconstructions ............................................ 11

2.1.3 Time Integration .............................................................. 13

2.1.4 Monotonicity Limiter ........................................................ 13

2.2 SV Method .............................................................................. 15
CHAPTER 3

EVALUATION OF DG AND SV METHODS ................................................ 19

3.1 Number of Operations and Memory Requirement ............................... 19

3.1.1 DG Method .................................................................... 20

3.1.1.1 Numbers of operations .............................................. 20

3.1.1.2 Memory requirement ............................................... 22

3.1.2 SV Method ..................................................................... 23

3.1.2.1 Numbers of operations ................................................ 23

3.1.2.2 Memory requirement ............................................... 25

3.2 Accuracy and CPU Times ............................................................ 25

3.2.1 Scalar Conservation Laws ................................................... 25

3.2.1.1 2D linear wave equation ............................................ 25

3.2.1.2 2D Burger’s equation ............................................... 30

3.2.2 Euler Equations ............................................................... 33

3.2.2.1 Vortex propagation problem ....................................... 33

3.2.2.2 Double Mach reflection problem .................................. 38

3.3 Conclusions ........................................................................... 42

vi

CHAPTER 4

SV METHOD FOR THE DIFFUSION EQUATION ......................................... 43
4.1 Three Formulations ..................................................................... 43
4.1.1 Naive SV Formulation ......................................................... 43
4.1.2 Local SV Formulation ......................................................... 46
4.1.3 Penalty SV Formulation ...................................................... 48
4.2 Fourier Analysis .......................................................................... 49
4.2.1 Naive SV Formulation ........................................................ 53
4.2.2 Local SV Formulation ........................................................ 56
4.2.3 Penalty SV Formulation ...................................................... 59
4.3 Conclusions ............................................................................... 61
CHAPTER 5
EXTENSION OF SV METHOD TO CONVECTION DIFFUSION EQUATION S.... 63
5.1 Formulation for 1D Convection-Diffusion Equation ................................ 63
5.2 Formulation for 2D Linear Convection-Diffusion Equation. . . . . 66
5.3 Accuracy Study ........................................................................... 68
5.3.1 1D linear Convection-Diffusion Equation ................................... 68
5.3.2 1D Viscous Burger’s Equation ................................................. 69
5.3.3 1D Fully Nonlinear Case ...................................................... 69
5.3.4 2D Linear Convection- Diffusion Equation ................................. 71
5.4 Conclusions ............................................................................... 71
CHAPTER 6
EXTENSION OF SV METHOD TO THE NAVIER-STOKES EQUATIONS. . . . . 72
6.1 Formulation for Navier-Stokes Equations ............................................. 72
6.2 Numerical Experiments .................................................................. 76
6.2.1 The Couette Flow ............................................................... 76
6.2.2 Laminar Flow along a Flat Plate ............................................. 81
6.2.3 Subsonic Flow over 3 Circular Cylinder ..................................... 91
6.2.4 Laminar Subsonic Flow around the NACA0012 Airfoil .................. 98
6.3 Conclusions .............................................................................. 106
CHAPTER 7
SUMMARY AND FUTURE WORK .......................................................... 108
APPENDIX A
GAUSSIAN QUADRATURE FORMULAS .................................................. 112
APPENDIX B
PUBLICATIONS ASSOCIATED WITH THIS WORK .................................... 116
BIBLIOGRPHY ................................................................................... 117
VITA ................................................................................................ 123

vii

LIST OF TABLES

CHAPTER 3
Table 3.1 Number of operations for the DG method ........................................ 22
Table 3.2 Number of operations for the SV method ......................................... 24
Table 3.3 Errors and CPU times at t = 1 for a 2D linear wave equation using the DG
method on the regular mesh ......................................................... 28
Table 3.4 Errors and CPU times at t = I for a 2D linear wave equation using the SV
method on the regular mesh ......................................................... 28
Table 3.5 Errors and CPU times at t = I for a 2D linear equation using the DG method
on the irregular mesh ................................................................. 29
Table 3.6 Errors and CPU times at t = I for a 2D linear wave equation using the SV
method on the irregular mesh ....................................................... 29
Table 3.7 Errors and CPU time on the 2D Burger’s equation at t = 0.1 using DG on the
irregular mesh ......................................................................... 31
Table 3.8 Errors and CPU time on the 2D Burger’s equation at t = 0.1 using SV on the
irregular mesh ......................................................................... 31
Table 3.9 Errors and CPU time on the 2D Burger’s equation at t = 0.45 using DG on the

irregular mesh ......................................................................... 32

Table 3.10 Errors and CPU time on the 2D Burger’s equation at t = 0.45 using SV on

Table 3.11

Table 3.12

Table 3.13

Table 3.14

the irregular mesh ................................................................... 32

Errors and CPU times for the propagating vortex case at t = 10 using the DG
method on the regular mesh ....................................................... 36

Errors and CPU time for the propagating vortex case at t = 10 using the SV
method on the regular mesh ....................................................... 36

Errors and CPU times for propagating vortex case at t = 10 using the DG
method on the irregular mesh ..................................................... 37

Errors and CPU times for the propagating vortex case at t = 10 using the SV
method on the irregular mesh ..................................................... 37

viii

CHAPTER 4

Table 4.1 L1 and L.>0 errors and orders of accuracy based on the local SV formulation
for the diffusion equation ............................................................ 47

Table 4.2 L1 and L0,, errors and orders of accuracy based on the penalty SV formulation
for the diffusion equation ............................................................ 49

CHAPTER 5
Table 5.1 L1 and L0° errors and orders of accuracy on 1D linear convection-diffusion
equation (at t = 1.0) .................................................................. 68

Table 5.2 L1 and L,o errors and orders of accuracy on 1D viscous Burger’s equation
(at t = 1.0) .............................................................................. 70

Table 5.3 L1 and L0,, errors and orders of accuracy on the 1D fully nonlinear
convection-diffusion equation (I = 1.0) ........................................... 70

Table 5.4 L1 and Lo0 errors and orders of accuracy on the 2D linear convection-
diffusion equation (I = 1.0) ......................................................... 71

CHAPTER 6

Table 6.1 L1 and L0o errors and orders of accuracy on Couette flow
(Density in Estimation 1) ............................................................. 79

Table 6.2 L1 and Loo errors and orders of accuracy on Couette flow
(Temperature in Estimation l) ...................................................... 79

Table 6.3 L1 and LC,o errors and orders of accuracy on Couette flow
(Density in Estimation 2) ............................................................ 80

Table 6.4 L1 and L.>0 errors and orders of accuracy on Couette flow
(Temperature in Estimation 2) ....................................................... 80

Table 6.5 Pressure part (C d, p ) and viscous Part (Cdy ) of the drag coefficient for the
NACA0012 Airfoil ................................................................. 105

Table 6.6 Comparison of pressure pa1t (Cd,p ) and viscous part (C d", ) of the drag

coefficient for the NACA0012 airfoil computed by the present method with
various structured and unstructured solvers ...................................... 106

ix

APENDIX

Table A.l Christoffel weights W, and roots 5, ,i = l, 2,..., n , for Legendre polynomials of
degree 1 to 6 ......................................................................... 113

Table A2 Weights and evaluation points for integration on triangles ................... 115

LIST OF FIGURES

CHAPTER 2
Figure 2.1 Degrees of freedom for DO ........................................................ 12
Figure 2.2 Spectral volumes of various degrees .............................................. 15
Figure 2.3 SV partitions used in numerical simulations ..................................... 18
CHAPTER 3
Figure 3.1 Regular and irregular grids ......................................................... 27

Figure 3.2 Computational domain and boundaries for double Mach reflection
problem ................................................................................ 38

Figure 3.3 Density contours computed with second order DG scheme using a TVD
limiter ................................................................................... 40

Figure 3.4 Density contours computed with second order SV scheme using a TVD

limiter .................................................................................. 41
CHAPTER 4
Figure 4.1 1D linear spectral volumes ......................................................... 43

Figure 4.2 The numerical solutions versus the exact solution using the linear
reconstruction based on the Naive SV formulation for the diffusion
equation ............................................................................... 45

Figure 4.3 The numerical solutions versus the exact solution using the quadratic
reconstruction based on the Naive SV formulation for the diffusion

equation ............................................................................... 45
CHAPTER 5
Figure 5.1 One-dimensional spectral volumes ................................................ 64
CHAPTER 6
Figure 6.1 Computational grids for Couette flow case ...................................... 77

xi

Figure 6.2 Convergence trend of numerical solution to the steady exact solution in

Couette flow test ..................................................................... 78
Figure 6.3 The computational domain for the flat plate case ............................... 81
Figure 6.4 Meshes for the flat plate case ...................................................... 83
Figure 6.5 u-velocity profiles for the flat plate case ......................................... 84
Figure 6.6 Skin friction coefficients along a flat plate ...................................... 87
Figure 6.7 Residual history of continuity equation with cubic reconstruction ........... 90
Figure 6.8 Grid for the case of subsonic flow over a circular cylinder ................... 92

Figure 6.9 Instantaneous Mach contours for Ma 2 0.2 flow
over a circular cylinder at Re = 75 ................................................ 93

Figure 6.10 Instantaneous entropy contours for Ma = 0.2 flow
over a circular cylinder at Re = 75 ............................................... 94

Figure 6.11 Instantaneous vorticity contours for Ma = 0.2 flow
over a circular cylinder at Re = 75 ............................................... 94

Figure 6.12 Pressure histories for Ma = 0.2 flow
over a circular cylinder at Re = 75 .............................................. 95

Figure 6.13 Computational grid for the NACA0012 test case .............................. 98

Figure 6.14 Mach isolines around the NACA0012 airfoil
(Re = 5000, Ma = 0.5, a = 0°) .................................................... 100

Figure 6.15 Pressure coefficient distribution along the N ACA0012 airfoil
(Re = 5000, Ma = 0.5, a = 0°) ................................................... 103

Figure 6.16 Skin friction coefficient distribution along the NACA0012 airfoil
(Re = 5000, Ma = 0.5, a = 0°) .................................................... 104

Figure 6.17 Residual history of continuity equation with cubic reconstruction for
NACA0012 test case .............................................................. 105

xii

”U

‘I

{own mshfln<=booN-nm©

kfiQJ
{O

NOMENCLATURE

flow variable, or conservative flow vector
total energy

flux vector

flux in x direction

flux in y direction

density

velocity component in x direction
velocity component in y direction
pressure

absolute temperature

dynamic viscosity coefficient
ratio of specific heat

specific heat at constant pressure,

Prandle number

a space domain
boundary of domain £2
a test function, or a scalar limiter

order of interpolation polynomial

dimension of the interpolation polynomial space in 2D
number of control volumes in a spectral volume

upper bound of time

time

time step

position vector

gradient operator

the j-th control volume in the i-th spectral volume

volume of the control volume C ,3 j
number of faces in C ,3 j

the r-th face in C i, J-

shape function in a spectral volume
triangular coordinate

Lebesgue constant
the outward unit normal vector

xiii

CHAPTER 1

INTRODUCTION

1.1 Background of Computational Fluid Dynamics (CFD)

The physical aspects of transport phenomena in the macro-scale are governed by the
Newton’s laws of motion and the fundamental principles of mass, energy and species
conservation. The final objective of most engineering investigations is to obtain a
quantitative description of the physical problem by analytical, experimental or numerical

methods.

By the turn of the twentieth century, the development of closed form analytical solutions
for flow field problems had reached a highly mature stage and it was being realized that a
large class of problems still remained which were not amenable to exact analytical
solution methods. Experimental fluid dynamics has played an important role in validating
and delineating the limits of the various approximations to the governing equations. The
wind tunnel, as a piece of experiment equipment, provides an effective means of
simulating real flows. Traditionally this has provided a cost-effective alternative to full-
scale measurement. In‘the design of equipment that depends critically on the flow
behavior, e.g. aircraft design, full-scale measurement, as part of the design process is
economically unavailable. The steady improvement in the speed of computer and
memory size since 1950s has led to the emergence of computational fluid dynamics
(CFD) to study the characteristics of fluid dynamics using digital computers. CFD is

significantly cheaper than wind tunnel testing and will become even more so in the future.

1.2 Numerical Methods in CFD

The success of CFD is really dependent on two factors, i.e. improvement on computer
hardware and highly efficient computational algorithms. Therefore, there have been
intensive efforts to develop highly efficient and accurate numerical algorithms to seek

higher quality numerical solutions with less CPU time.

The basic issue of quality of numerical solutions in CFD simulation is fundamentally
important: how accurate are the numerical simulations and how does one obtain the most
accurate results given a fixed computational resource? These questions lie at the core of
modern numerical methods that aim to control the error in the computed solution and to

optimize the computational process.

Many methodologies have been developed to address this issue for the hyperbolic
systems in the last three decades. One of the most successful algorithms is the Godunov
method [24], which laid a solid foundation for the development of modern upwind
methods [20,25-26,46,60-61]. For example, van Leer [60-61] extended the first-order
Godunov method to second-order by using a piece-wise linear data reconstruction and a
limiter to remove spurious numerical oscillations near steep gradients. In addition, for
better efficiency the exact Riemann solver used in the Godunov method was replaced by
approximate Riemann solvers or flux-splitting procedures, such as the flux-vector
splitting [59] by Steger and Warming, the flux-difference splitting [52] by Roe, the
smoother flux vector splitting [62] by Van Leer, the differentiable approximate Riemann

solver [45] by Osher, and AUSM [40] by Liou, FUSS [63] by Wang, among many others.

One of the most popular schemes for obtaining solutions on unstructured meshes is the
discontinuous Galerkin finite element (DG) method, which was introduced in the early
1970’s for the numerical solution of first-order hyperbolic problems (see [6,13,15-18,23,
35,36,39,48-50]). Simultaneously, but independently, it was proposed as non-standard
schemes for the numerical approximation of second-order elliptic equations [2, 68]. In
recent years there has been renewed interest in the discontinuous Galerkin method due to
its favorable properties, such as a high degree of locality, stability in the absence of
streamline-diffusion stabilization for convection-dominated diffusion problem [29], and
the flexibility of locally varying the polynomial degree in hp-version approximations,
since no point wise continuity requirements are imposed at the element interfaces. Much
attention has been paid to the analysis of the DG method applied to non-linear hyperbolic
equations and hyperbolic systems [11,12,27], several other types of non-linear equations
(including the Hamilton-Jacobi equation [30], and non-linear Schrodinger equation [37],
and other non-linear problem [14]). Also, it was extended to the compressible Navier-

Stokes equations [7].

An alternative to the finite element method is the finite-volume method, in which the
governing equations are solved in integral form over the discrete volumes formed by the
cells of a mesh. Description of various finite—volume methods on unstructured meshes are
given by Barth and Jesperson [5], Whitaker, et a1. [69], Jameson, et a]. [32-34], and
Mariplis and Jameson [42]. Barth [3] presents a detailed account of the implementation of

finite volume schemes for the Euler and Navier-Stokes equations using efficient edge-

based data structures. Finite volume schemes generally solve for quantities averaged over
cells of the actual mesh in the case of cell-centered schemes or over cells of a dual mesh
in the case of vertex schemes. In any event, in order to evaluate the residual, a
polynomial data distribution must be reconstructed from these averaged quantities. To
achieve higher than second order accuracy, a higher order distribution must be
constructed in each cell, requiring information from more distant neighbors. This was
done by Barth and Frederickson [4] for quadratic reconstruction (and hence third order
accuracy). Hu and Shu [31] further devised a fourth order scheme without expanding the

third order stencil.

More recently, a high-order, conservative, yet efficient method named the spectral
volume (SV) method was presented by Wang [64] for hyperbolic conservation law. The
SV method is a finite volume method, in which the concept of a “spectral volume” is
' introduced to achieve high-order accuracy in an efficient manner similar to spectral
element and multidomain spectral methods. Each spectral volume is further subdivided
into control volumes, and cell-averageddata from these control volumes are used to
reconstruct a high-order approximation in the spectral volume. Then Riemann solvers are
used to compute the fluxes at spectral volume boundaries. Cell-average state variables in
the control volumes are updated independently. Furthermore, total variation diminishing
and total variation bounded limiters are introduced in the SV method to remove/reduce
spurious oscillations near discontinuities. Unlike spectral element and multidomain
spectral methods, the SV method can be applied to fully unstructured grids. A very

desirable feature of the SV method is that the reconstruction is carried out analytically,

and the reconstruction stencil is always nonsingular, in contrast to the memory and CPU-

intensive reconstruction in a hi gh-order k-exact finite volume method.

1.3 Motivation and Obiectives of This Study

The newly developed spectral (finite) volume method has been successfully
demonstrated for hyperbolic conservation laws including non-linear systems on
unstructured grids in a series of papers [65-67]. A framework has been established to
easily solve non-linear time-dependent hyperbolic systems of conservation laws using
explicit, non-linear Runge-Kutta time discretization [58] with approximate Riemann
solvers and TVB (total variation bounded) non-linear limiters [54]. One objective of this
study is to give a further numerical demonstration that the SV method is comparable to
other high order methods and also possesses some unique properties. To do so, we
evaluate the DG and SV methods on hyperbolic conservation laws, since the DG and SV
methods seem to be the most efficient among the high-order methods on unstructured

grids.

Ultimately, we wish to extend the SV method to the Navier-Stokes equations to perform
large eddy simulation and direct numerical simulation of turbulence flow for problems
with complex geometries. So, another objective of this study is to extend the SV method

[64-67] to the Navier-Stokes equations.

A key in the extension is to properly discretize the second order viscous terms. In a

second-order finite volume method, the solution gradients at an interface are computed

by averaging the gradients of the neighboring cells sharing the face, and were shown to
be adequate. For higher-order elements, special care has to be taken in computing the
solution gradients. For example, Cockbum and Shu developed the so-called local
discontinuous Galerkin method to treat the second order viscous terms and proved
stability and convergence with error estimates [l9] motivated by the successful numerical
experiments of Bassi and Rebay [7]. Baumann and Oden [8], Oden, Babuska and
Baumann [44] introduced a different discontinuous Galerkin method for the discretization
of the second order viscous terms. Riviere, Wheeler and Girault [51] analyzed three
discontinuous Galerkin approximations for solving elliptic problems in two or three
dimensions. More recently, Shu [57] summarized three different formulations of the
discontinuous Galerkin method for the diffusion equation, and Zhang and Shu [72]

performed a Fourier type analysis for these three formulations.

Motivated by the DG approach in handling the viscous term, three SV formulations for
pure diffusion equations will be presented in this research, and one of them will be
successfully applied to 1D and 2D scalar convection-diffusion equations, eventually to
viscous flows governed by the Navier-Stokes equations. The spatial convergence rate of
the SV method will be established for some scalar cases and Couette flow, and the

designed order of accuracy will be studies.

1.4 Outline of the Dissertation
The dissertation is arranged as follows. We first review the framework of the

discontinuous Galerkin (DG) spectral volume (SV) methods in Chapter 2. Then in

Chapter 3, we evaluate the DO and SV methods in terms of the number of operations,
memory requirement, accuracy and CPU times for inviscid flows. In Chapter 4, we
present three SV formulations for the diffusion equation. The extension of the SV method
to the viscous flow is presented in Chapter 5 and Chapter 6. Finally, a summary of the

present study and recommendations for further investigations are given in Chapter 7.

CHAPTER 2

FRAMEWORK OF DG AND SV METHODS

The DG method is a finite element method using discontinuous solution and test spaces
(usually piecewise polynomials of suitable degree), which means that the state variables
are not continuous across element boundaries. The fluxes through the element boundaries
are then computed using an approximate Riemann solver, mimicking the successful
Godunov finite volume method [24]. Due to the use of Riemann fluxes across element
boundaries, the DG method is fully conservative at the element level. The SV method
[64-67] is a finite volume method. For a given unstructured grid, each element (called a
spectral volume) is further partitioned into structured sub elements named control
volumes (CVs). Mean state-variables at the CVs inside a SV are employed to construct a
high-order polynomial within the element or SV, which is then utilized to update the
means at the CVs. The reconstruction problem can be solved analytically, and is identical
for all simplexes. Therefore a high-order SV method is much more efficient than a high-
order k-exact FV method, in which a reconstruction problem must be solved for each
control volume. The SV method is fully conservative at the sub-cell control volume level.

Both methods are reviewed next.

2.1 DG Method

Consider the following two-dimensional conservation laws

Q,+VoF=0, 52x(0,r) (2.1)

equipped with proper initial and boundary conditions. In Eq. (2.1), F = (f, g) is the flux
vector. Multiplying Eq. (2.1) by a test functiongo, integrating over the computational

domain £2, and performing integration by parts, we obtain the following weak statement

of the problem
j¢Q,dV+ §¢F(Q)-ndS—[V¢-F(Q)dv =0,\7’(o (2.2)
$2 89. .0

Note that the integral in Eq. (2.2) is understood to be performed in a component-wise

manner if Q is a column vector.

2.1.1 Space Discretization

Assume that the computational domain [2 is subdivided into N non-overlapping triangular
elements {Ti}. By applying Eq. (2.2) to each element T,, we can obtain the discrete
analogue of Eq. (2.2) on the computational grid. Let the solution and test function be
piece-wise polynomials in each element. Denote the polynomial basis

as §(r)={g‘l(r),~--,§n(r)}T. If the polynomial is of order k, the dimension of the

polynomial space in 2D is n = (k+1)(k+2)/2. The solution and the test function on

element T,- can be expressed as

0,-0.0 = ZQ/(IK‘J-(r), (0;. = 20%;“)- (2.3)
j=l

i=1
The expansion coefficients Qij denote the degrees of freedom (DOFs) of the numerical

solution on element Ti. Note that there is no global continuity requirement for Q,- , which

is generally discontinuous across the element boundaries. Using the solution and test

function, Eq. (2.2) on element T,- becomes

i IindV + §¢hF .nds — [wk . m = 0. (2.4)
dt
Ti 8T,- Ti

Equation (2.4) must be satisfied for any test function 01h. Since 5]- is the basis function

for (0,, , Eq. (2.4) is equivalent to the following system of n equations

d .
EjijidV+ §§jFondS— jvgj-dew, 131571. (2.5)
Ti 8T,- Ti
Because the approximate solution is discontinuous at the element boundaries, the

interface flux is not uniquely defined. It is at this stage the Riemann flux used in the

Godunov finite volume method [24] is borrowed. The interface flux function Fon is

replaced by a Riemann flux E(QL,QR,n), where QLand QR are the state variables at

the left and right side of the interface. In order to guarantee consistency and conservation,

the Riemann flux must satisfy

fi<Q,Q,n> =F<Q>-n, fi<QL.QR,n> =—fi<QR,QL.—n>. (2.6)
The surface and volume integrals in Eq. (2.5) can be computed with Gauss quadrature
formulas of suitable orders of accuracy, which are given in Appendix A. Following the
arguments given in [13], the surface integral must be exact for polynomials of degree 2k,
while the volume integral must be exact for polynomials of degree 2k-I , i.e.,

K
§tijondS = Z jij-ndS.

aTi r=1Ar

Ing.ndS z: Zwrsgjfivzs)fi(QL(rrs)aQR(rrs)rnr)A,-a
A, 5:1

10

[V5]. .de z ZwsV§j(rs)-F(Q,- (rs))V,-. (2.7)
T.

1 5:1

where K is the number of planar faces of Ti, ns is the number of quadrature points on a
planar face for the surface integral, nv is the number of quadrature points in the element

for the volume integral, Wrs and ws are the Gauss quadrature weights, rm and rs are the

Gauss quadrature points. Let Ui = {Q,~l,---,Q,-" }T be the DOFs for element Ti, and Wi

denote the mass matrix 61- fildV . Equation (2.5) can be further written as
Ti

i . _1

i(Iii-4W) [grands-jvgordv =0. (2.8)
t 8T,- Ti

By assembling together all the elemental contributions, a system of ordinary differential

equations that govern the evolution of the discrete solution can be written as

it: = R(U) , (2.9)
(it

where U is the global vector of DOFs, and R( U) is the global residual vector with the

element vector being

Ri(U)=—(Wi)_l [gr-ndS— jvg-FdV . (2.10)
8T,- Ti

2.1.2 Grids and Data Reconstructions
The degrees of freedom are chosen as the values at certain points in each element, which

are shown in Figure 2.1.

11

f
,r’ \ '
o- » é c \ o

 

.3

 

 

(a) (b) (C)

Figure 2.1 Degrees of freedom for DG,
(a) linear element, (b) quadratic element, (0) cubic element

The basis functions in terms of the triangular coordinates are given next.

Linear element:

{1:41,
52:42,
53:1—41—42-

Quadratic element:

951:41‘041—1),

§2=flq-(24/2-1),
53 ”1342434),
54:441'42,
5531/12/13,
56:443'41-

Cubic element.-
51:41'(341-1)'(341—2)/2.
52 =42 I342 -1)°(3/iq -2)/2,
53 ”1313/13 -l)'(3/13 -2)/2.
54:41'42'(341-1)'9/2,
§s=xizvi3-(3AQ—1)-9/2,
56 ”1341 '(343 —1)-9/2,

12

:7 wry/12 (3,12 —1)-9/2,
«fa =22 4313/13 —1)-9/2,
59:23.21-(3/11—1)-9/2,
510:2741'42'43,

where xij,j=l,2,3 are the triangular coordinates described in Appendix A,

11+12+A3 =1,/1j20.

2.1.3 Time Integration
An explicit multi-stage third-order TVD (total variation diminishing) Runge-Kutta
scheme is employed for time integration [55]. The Runge-Kutta scheme can be expressed
in the following form:

U“) = U" +AtR(U");

(1(2) 2%U" +iw“) +AtR(U(l))]; (2.11)

U"+1 = gr!" +§[U(2) +AtR(U(2))].

Where, U n is the global vector of DOFs at timet = tn , and U "+1 is the global vector of

DOFs at time! = t,,+1 = In + At.

2.1.4 Monotonicity Limiter
For the non-linear Euler equations, it is necessary to perform data limiting to maintain
stability if the solution contains discontinuities. There are two possible ways of applying

limiters in the system setting. One way is to apply a limiter to each characteristic variable.

13

The other is to apply a limiter to each of the conservative variables. In one dimension, the
former has the nice property of naturally degenerating to the scalar case if the hyperbolic
system is linear. In multiple dimensions, characteristic variables are defined in a
particular direction, e. g. in the face normal direction. In a fully unstructured grid, there is
no coordinate direction to define a characteristic variable. Therefore it is difficult to
design characteristics-based limiters in multiple dimensions. In this research, we choose
the component-wise approach for the limiter, which should also be much more efficient
than the characteristic approach. To this end, we first establish the following numerical

monotonicity criterion for each element

min —max
Qi —<- Q; (rs) S Q 9 (2.12)
where 2mm and aim“ are the minimum and maximum cell-averaged solutions among

all its neighboring elements sharing a face with T), and Qi(rs) is the solution at any of

the quadrature points. If the inequality (2.12) is violated for any quadrature point, then it
is assumed that the element is close to a discontinuity, and the solution in the element is

forced locally linear, i.e.,
Q,(r)=§,+VQ,-(r-r,-), Vrer}, (2.13)
where r,- is the position vector of the centroid of T,-. The magnitude of the solution

gradient is maximized subject to the monotonicity condition given in the inequality
(2.12). The original polynomial is used to compute an initial guess of the gradient at the

element centroid, i.e.,

922a]

VQi=[dx , 8y

rt

14

This gradient may not satisfy the inequality (2.12). Therefore it is limited by multiplying

a scalar limiter are [0, 1] so that the following solution satisfies the inequality (2.12)

0.10:5.- +¢VQ.--(r-r. ). (2.14)
The scalar limiter can be obtained by examining the numerical solutions at all the

quadrature points [65].

2.2 SV Method
In the SV method, the element T) is named a spectral volume, which is further partitioned

into subcells named control volumes (CVs), indicated by C as shown in Figure 2.2.

1.)»

Linear reconstruction

 

 

I) “A Quadratic reconstruction
(I

 

 

2) {1“ Cubic reconstruction

(l

 

Figure 2.2 Spectral volumes of various degrees

To approximate the solution as a polynomial of degree k in two dimensions (2D), we

need to partition the SV into n = (k+1)(k+2)/2 sub-cells. The degrees of freedom (DOFs)

in a SV are the volume-averaged mean variables -Q—,j at the n CVs. There are numerous

ways of partitioning a SV, and not every partition is admissible in the sense that the
partition may not be capable of producing a degree k polynomial. Once n mean solutions
in the CVs of an admissible SV are given, a unique polynomial reconstruction can be
obtained from
n
mm = ZLj-(r)§,-,j , (2.15)
j=l

where L j (r) are also degree k polynomials satisfying

j L", (r)dV = VLJ-ij, (2.16)
Ci,j

and VLJ-is the volume of C13!” This high-order polynomial reconstruction facilitates a

high-order update for the mean solution of each CV. Integrating Eq. (2.1) in each CV, we

obtain

rig—2,, j

K
TV,”- + Z jar-mars = 0, (2.17)

r=lAr
where K is the total number of faces in C13!” The flux integral in Eq. (2.17) is then

replaced by a Gauss-quadrature formula (see Appendix A) that is exact for polynomials

of degree k
ne
[ (F on)dS = Z w,,F(Q(r,, )) . n,A,, (2.18)
Ar 3:1

16

where ne is the number of quadrature points on the r-th face, wrs are the Gauss
quadrature weights, rrs are the Gauss quadrature points. Since the reconstructed

polynomials are piece-wise continuous, the solution is discontinuous across the
boundaries of a SV, although it is continuous across interior CV faces. The fluxes at the
interior faces can be computed directly based on the reconstructed solutions at the
quadrature points. The fluxes at the boundary faces of a SV are again computed using
approximate Riemann solvers given the left and right reconstructed solutions. The

Runge-Kutta scheme is again used for time integration.

The TVD limiter in the SV method [65] is very similar to the one described in the last
section. The main difference is that the limiter is applied for the sub-cell averaged state
variables, rather than for the averaged state variables of macro element, i.e., the SV. This
is possible because of the inherent local resolution in the SV method. In order to make an
objective comparison with the DG method, the limiters are implemented in a similar

fashion.

Remark: In Wang’s paper [65], the Lebesgue constant is employed to quantify the

quality of the reconstructions. Following the paper [65], the Lebesgue constant "I‘m" is

expressed as

n
F = L-
" nll 122131.24 10>

for a given partition of a triangle E. It was shown in paper [65] that the smaller the
Lebesgue constant, the better the interpolation polynomial. In this work, we use good

enough SV partitions so that the interpolation polynomial is convergent when the

17

computational grid is refined. Figure 2.3 shows the SV partitions used in our numerical

simulations for linear, quadratic, and cubic reconstructions. The order of grid nodes, faces,

control volumes and Gauss quadrature points (GQPs) are also shown in one spectral
volume.

 

 

 

 

  
  
  
    

 

No. of nodes
No. of face:
No- of CV:

No. of GOP:

 

 

 

 

 

 

(C)
Figure 2.3 SV partitions used in numerical simulations
(a) linear SV, (b)quadratic SV,(c) cubic SV
(this figure presented by color)

18

CHAPTER 3

EVALUATION OF DG AND SV METHODS

3.1 Number of Operations and Memory Requirement

In order to provide a reasonable estimate of the number of operations for both methods,
we need to specify the governing equation and the Riemann solver. Two equations are
considered in this research. One is the 2D scalar linear conservation law in which Q is a
scalar, and f = aQ and g = bQ with a and b being constants. The other equation is the
2D Euler equations. In both cases, the Rusanov flux [53] (also called local Lax-
Friedrich’s flux) is selected. The Rusanov flux for the scalar conservation laws takes the
following form

QL(anx+bny) if(anx+bny)>0

R (3.1)
Q (an x + bn y) otherwise

fi(QL.QR.n)={

Since modern computers can execute multiplications as fast as additions, 1 Operation is
defined to be one multiplication or one addition. Internal functions such as sqrt is
assumed to cost 10 operations. In addition, each if statement is also counted as 1

operation. In this case, this scalar Riemann solver costs M R = 5 operations (3 operations

to compute an x + buy, one if statement, and another multiplication). The analytical flux

takes M a operations.

For the Euler equations, the Rusanov Riemann flux [53] takes the following form

fi<QL.QR,n) = %{[F(QL)+ F(QR)]-n —a(QR -QL)} (3.2)

19

where a=|r7n|+E, 17,, is the average face normal velocity, and E the average speed of

sound at the interface. Given the vector of conservative variables, it is estimated that an

analytical flux evaluation costs M a = 24 operations, and a Riemann flux takes M R = 85

operations.

For simplicity we have not considered the cost of limiters in the number of operations.
We do believe that the limiters in the SV method are more expensive to be implemented
than those in the DG method because data limiting is carried out for each element in the

DG method, but for each sub-cell (CV) in the SV method.

3.1.1 DG Method
We consider linear, quadratic and cubic elements, which are expected to yield second,

third and fourth-order spatial accuracy respectively. The DOFs for these elements are

shown in Figure 2.1. Over each element Ti, the residual vector can be written as

[(F.v§1)dV— §E§1dS
. .-1 Ti 3T1"
R' (U) = -—(W') 5 . (3.3)
j(F-V.§,,)dV— [fiends
_T,' 8T;

 

 

3.1.1.1 Numbers of operations

 

The total number of operations can be roughly divided into three main parts,

corresponding to the cost for computing the state variables at all the Gauss quadrature

20

points (N 1), the number of operations to compute the fluxes (N2), and the cost to multiply

the mass matrix (N3).

There are a total of (nv+3 *ns) quadrature points that are used for surface and volume
integrals. We need It multiplications and n-1 additions to compute one state variable
given the DOFs. Assume Q has nc component. Then the total number of operations to
compute the solutions at all the quadrature points is then
N1=nc*(2*n-l)*(nv+3*ns). (3.4)
Note that we have ignored the number of operations to compute the limiter for simplicity.
To evaluate the volume integral, we need to compute the (analytical) fluxes at nv
quadrature points relating to n shape functions, while 3*ns Riemann fluxes are necessary
to evaluate the surface integral. However Riemann fluxes are shared between two
neighboring elements. Therefore we need to halve the number of operations for the
Riemann fluxes when evaluating the number of operations per element. We also need to
include the number of operations to carry out the Gauss quadrature formula. Thus we
obtain
N2 =n*nv*Ma +3*ns*MRl2+nc*n*(2*nv—l)+nc*n*(3*ns—l) (3.5)
N3 is simply the cost of a square matrix multiplying a vector, which is 2*n*n-n for one
component. For no components, we therefore have
N3 = nc*(2 *n *n-n). (3.6)

Note that N 3 = 0 if an orthogonal basis is used.

The total cost to compute the residual vector for a single element is then

21

NT=N1+N2+N3 (3.7)
The numbers of operations for the DG schemes of second to fourth orders are listed in

Table 3.1.

Table 3.1 Number of operations for the DG method

 

 

 

 

 

 

 

Equation k N nv Ns NT
1 3 3 2 141
2 6 6 3 512
Scalar Conservation Law 3 10 12 4 1496
1 3 2 831
2 6 6 3 2627
Euler Equations 3 10 12 4 7334

 

 

 

 

 

 

 

 

3.1.1.2 Memory requirement

The memory requirement for the DG method is estimated as the following:

0 Two solutions; one at the current time step, and the other at the last time step.

Residual.

0 Volume, centroid coordinates, face area, and face unit normal.
0 Coordinates of quadrature points (face, cell)
0 Gradient of shape function on quadrature points (face, cell).

0 Shape functions, and their gradients at the centroid of the elements.

22

The storage requirement is roughly 90 words per element for a second-order DG scheme,
221 words per element for a third-order DG scheme, and 512 words per element for a

fourth-order DG scheme for the 2D Euler equations.

3.1.2 SV Method

The degrees of freedom in the SV method are the mean state variables at the sub-cells,

i.e. control volumes. Over each spectral volume 7} , the residual can be expressed as

- ffids
. 3C”
R' (U) = s (3.8)
— §13dS
L 8C”,

 

 

3.1.2.1 Numbers of operations

There are two kinds of faces in a spectral volume. The faces that lie on the SV boundaries
are called Riemann faces, because the state variables are discontinuous across these faces.
The other faces that lie inside a SV are named continuous faces because the state
variables are continuous across these faces. Denote the total number of faces in a SV with
nf, and the number of Riemann faces nr. Then the number of continuous faces is then (nf
- nr). Let the number of quadrature points on each face (edge) be ne. Then the number of
operations to compute the state variables at all the quadrature points is nc*nf*ne*(2n-I).
In addition, a total of (nf - nr)*ne analytical fluxes need to be computed while nr*ne
Riemann fluxes must be computed. Since the Riemann faces are shared between two

neighboring SVs, the number of operations is again halved. We also include the number

23

of operations to carry out the Gauss quadrature formula (2 *ne-1)*nf*nc. Since the mass
matrix in the SV method is always the identity matrix, number of operations in the SV
method can be written as

NT =N1+N2+N3

where
N1=nc*nf*ne*(2n—l) (3.9)
N2 =(nf—nr)*ne*Ma +nr*ne*MR/2+(2*ne-1)*nf*nc (3.10)
N3=O (3.11)

The numbers of operations for the SV schemes of second to fourth orders are listed in

Table 3.2 for both the scalar and system conservation laws.

Table 3.2 Number of operations for the SV method

 

 

 

 

 

 

 

Equation k N ne nf nr NT
1 3 1 9 6 81
Scalar Conservation Law 2 6 2 15 9 468
3 10 2 27 12 1287
l 3 1 9 6 543
Euler Equations 2 6 2 15 9 2553
3 10 2 27 12 6168

 

 

 

 

 

 

 

 

 

24

3.1.2.2 Memory requirement

The permanent memory requirement is estimated as follows:

0 Two solutions, one at the current time step, and the other at the last time step;

0 The residual, volumes and centroid coordinates for the CVs, the face unit normal
vectors and areas for the sub-cell grid;

0 Face to cell and face to node connectivity for the sub-cell grid;

0 Coordinates of the sub-cell grid;

0 A connectivity linking each quadrature point on a face to a point of the local standard

SV to reconstruct the solution at the quadrature point.

The storage requirement for the 2D Euler equations is roughly 99 words per element for a
second-order SV scheme, 194 words per element for a third-order SV scheme, and 361

words per element for a fourth—order SV scheme.

3.2 Accuracy and CPU Tim
The following, we will do the comparison of DG and SV methods in terms of accuracy
and CPU time on both 2D scalar conservation laws and 2D inviscid flows governed by

Euler equations.

3.2.] Scalar Conservation Laws
3.2.1.1 2D linear wave equation
We first test the performance of both methods for the following linear scalar conservation

law

25

ut+ux+u =0; 0<x<2, 0<y<2

.Y
with u(x, y,0) = sin(7z(x+ )0) , and periodic boundary condition. The numerical
simulation was carried out until t = I on two different grids, one regular and one irregular
as shown in Figure 3.1. The finer meshes are produced recursively from the coarser
meshes by dividing each triangle into 4 smaller triangles. The third-order TVD Runge-
Kutta time integration scheme was used with a sufficiently small time step that the errors
are independent of the time step. The same time step was used in both the DG and SV
methods although larger time steps are permitted in the SV method for stability. The
errors are computed based on the cell-averaged state variable on the element or the SV.
No limiters were employed in the simulations since the problem is smooth. Tables 3.3
and 3.4 present the errors and CPU times using both methods on the regular mesh, while
Tables 3.5 and 3.6 display the errors and CPU times using both methods on the irregular
mesh. Note that on the regular mesh, both methods achieved the expected numerical

order of accuracy in both L1 and L0,, norms. However, the DG method consistently
produced L1 and Loo errors of smaller magnitude than the SV method. The SV method,

on the other hand, is about 15% - 140% faster than the DG method depending on the
order of accuracy of the scheme. On the irregular grid, the DG method is capable of

achieving the expected order of accuracy in both L1 and L0,, norms. Although the SV
method achieved the expected order of accuracy in L1 norm, the third-order SV scheme
showed a reduction of half an order in Loo norm. This may indicate the quality of the

quadratic SV partition can be further improved. Again, the SV method is consistently

faster than the DG method on the irregular grid.

26

 

(a) Regular (10x10x2)

 

 

 

 

 

 

(b) Irregular (10x10x2)

Figure 3.1 Regular and irregular grids

27

Table 3.3 Errors and CPU times at t = I for a 2D linear wave equation using the DG
method on the regular mesh

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Order of
accuracy Grid L1 error L1 order Loo error Loo order CPU (S)
10x 10x2 1. 146-02 --- 2.436-02 --- 3 .336-01
20x20x2 2.316-03 2.30 5.836-03 2.06 2.766+00
2 40x40x2 5.09e-04 2.19 1.42e-03 2.04 223er
80x80x2 l . 186-04 2.1 1 3 .496-04 2.02 1.816-l-02
160x 1 60x2 2.846-05 2.06 8.656-05 2.01 1.4564-03
10X 10x2 3.456-04 --- 7.656-04 --- 7.5906-01
20x20x2 4.276-05 3 .02 9.666-05 2.99 6.376+00
3 40x40x2 5.326-06 3 .00 1.216-05 3.00 5.096+01
80x80x2 6.65e-07 3 .00 1.51e-06 3.00 4.27e+02
160x l60x2 8.316-08 3 .00 1.896-07 3.00 3.376+03
10x 10x2 1.396-05 --- 2.436-05 --- 1.666+00
20x20x2 8.596-07 4.02 1.526-06 4.00 l.336+01
4 40x40x2 5.346-08 4.01 9.546-08 4.00 l.086+02
80x 80x2 3.336-09 4.00 5 .976-09 4.00 8.476+02
l60xl60x2 2.086-10 4.00 3.736-10 4.00 7.286+03

 

Table 3.4 Errors and CPU times at t = I for a 2D linear wave equation using the SV
method on the regular mesh

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Order of
accuracy Grid L1 error L1 order Loo error L0° order CPU (S)
10x 10x2 4.026-02 --- 5 .866-02 --- 1.216-01
20x20x2 1.066-02 1.92 1596-02 1.88 9.476-01
2 40x40x2 2.716-03 1.97 4.096-03 1.96 8.816+00
80X 80x2 6836-04 1.99 1.036-03 1.99 8.396+01
160x160x2 1.716-04 2.00 2.596-04 1.99 6.056+02
10X 10x2 3.736-03 -—- 5.216-03 --- 4.686-01
20x20x2 4.776-04 2.97 7.126-04 2.87 4.196+00
3 40x40x2 6.046-05 2.98 9.056-05 2.98 3.676+01
80x80x2 7.596-06 2.99 1.146-05 2.98 2.916+02
160x 160x2 9.516-07 3.00 1.436-06 2.99 2.216+03
10X 10x2 5.906-05 --- 8.406-05 --- 1.326-1-00
20X20X2 3.736-06 3.98 5.376-06 3.97 1356-1-01
4 40X40X2 2.356-07 3.99 3.346-07 4.01 8616-1-01
80x 80x2 1.486-08 3.99 2.096-08 4.00 6.906+02
160x 160X2 9.246-10 4.00 1.316-09 4.00 5.726+03

 

28

 

 

Table 3.5 Errors and CPU times at t = I for a 2D linear equation using the DG method

on the irregular mesh

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Order of
accuracy Grid L1 error L1 order Loo error L0,, order CPU (8)
10x10 2.17e-02 --— 6.05e-02 --- 3 .77e-01
20x20 4.676-03 2.22 1.646-02 1.88 3.266+00
2 40x40 1.076-03 2.12 4.176-03 1.97 2.656+01
80x 80 2.566-04 2.06 1.056-03 1.99 2.136+02
160x160 6.266-05 2.03 2.626-04 2.00 1.756+03
10x10 7.346-04 --- 2.566-03 --- 9.126-01
20x20 8.726-05 3.07 4.426-04 2.53 7.506+00
3 40x40 1.076-05 3.03 6.346-05 2.80 6.006+01
80x80 1.336-06 3.01 8.196-06 2.95 4.86-l-02
160x160 1.666-07 3.00 1.036-06 2.99 4.296-1-03
10x10 4.106-05 --- 1.806-04 --- 1.936+00
20x20 2.416-06 4.09 1.306-05 3.79 1.576+01
4 40x40 1 .476-07 4.04 8.546-07 3.92 1 .276+02
80x 80 9.086-09 4.01 5.726-08 3.90 1.026+03
160x160 5.656-10 4.01 3.726-09 3.94 8.876+03

 

Table 3.6 Errors and CPU times at t = 1 for a 2D linear wave equation using the SV

method on the irregular mesh

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Order of
accuracy Grid L1 error L1 order L0,, error L0,, order CPU (S)
10X10 6.716-02 --- 1.186-01 --- 1.386-01
20x20 1.836-02 1.87 3.406-02 1.80 1306-1-00
2 40x40 4.716-03 1.96 9.256-03 1.88 1.176+01
80x80 1 . 196-03 1.98 2.426-03 1.94 8.616+01
160x160 3.006-04 1.99 6.206-04 1.96 8.386+02
10x10 8.366-03 --- 1.686-02 --- 5.596-01
20x20 1.156-03 2.86 2.956-03 2.51 4.796+00
3 40x40 1.526-04 2.92 5.286-04 2.48 3.886+01
80X 80 2.016-05 2.91 1.316-04 2.01 3.276+02
160x160 2.646-06 2.93 2856-05 2.20 2716-1-03
10X 10 2.286-04 --- 7.396-04 --- 1.536+00
20x20 1.376-05 4.06 5.456-05 3.76 l.356+01
4 40x40 8.506-07 4.01 3.546-06 3.94 1.036+02
80x80 5.336-08 4.00 2.216-07 4.00 7.986+02
160x160 3.356-09 3.99 1.346-08 4.04 7.316+03

 

29

 

 

3.2.1.2 2D Burger’s equation
Consider the two-dimensional nonlinear wave equation:

ut+uux+uu =0, -l<x<l, —l<y<1

y
1 1 . . . . . .
u(x, y,0) = 4 + Esrn(7r(x + y)), wrth perrodrc boundary condition

We perform the simulation until t = 0.1 when the solution is still smooth. The errors and
CPU times are documented in Tables 3.7 and 3.8. We also test the performance of TVB
limiters [65], in which the simulation is performed until I = 0.45 when a shock wave
appeared in the solution. The solution errors in the smooth region {-02, 0.4] x [-0.2, 0.4]

are computed and presented in Tables 3.9 and 3.10.

30

Table 3.7 Errors and CPU time on the 2D Burger’s equation at t = 0.1 using DC on the

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

irregular mesh
Order of
accuracy Grid L1 error L1 order L0,, error L0,, order CPU (8)
10x10 1 .06e-02 --- 3.48e-02 --- 3.74e-02
20x20 2.75e-03 1.95 1.14e-02 1.61 3.22e~01
2 40x40 6.82e-04 2.01 3.21e-03 1.83 2.65e+00
80x80 1.70e-04 2.00 8.24e-O4 1.96 2.18e+01
160x160 4.24e-05 2.00 2.08e-04 1.98 l.74e+02
10x10 6.80e-04 --- 3. 17e-03 --- 1.80e-Ol
20x20 1.14e-04 2.57 8.32e-04 1.93 1.54e+00
3 40x40 1.79e-05 2.68 1.62e-04 2.36 l.29e+01
80x80 2.73e-06 2.71 3.45e-05 2.23 9.76e+01
160x160 4.08e-07 2.74 5.89e-06 2.55 7 .76e+02
10x10 6.01e-05 --- 4.58e-O4 --- 2.9le-01
20x20 3.68e-06 4.03 3.76e-05 3.61 2.34e+00
4 40x40 2.34e-07 3 .98 2.47e-06 3.93 l .93e+01
80x80 1.61e-08 3.86 2.07e-07 3.58 1.53e+02
160x160 1.20e-O9 3.75 1.95e-08 3.41 1.21e+03

 

Table 3.8 Errors and CPU time on the 2D Burger’s equation at t = 0.1 using SV on the

irregular mesh

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Order of

accuracy Grid L1 error L1 order L0,, error L0,, order CPU (S)
10X 10 5.796-03 --- 1.766-02 --- 1.346-02

20x20 1.466-03 1.99 4.916-03 1.84 1.296-01
2 40x40 3.676-04 1.99 1.416-03 1.80 1.176+00
80x80 9416-05 1.96 5.606-04 1.34 8.626+00

160x 160 2.396-05 1.97 2.686-04 1.06 6.926+01

10X 10 6.286-04 --- 2.156-03 --- 1.086-01

20x20 1.186-04 2.42 6.206-04 1.80 9.806-01
3 40x40 1916-05 2.62 1.386-04 2.17 7.996+00
80x80 3.026-06 2.66 3.286-05 2.07 6.836+01
160x 160 4.646-07 2.70 6.096-06 2.43 5.096+02

10x10 6.736-05 --- 4.436-04 --- 2.286-01
20x20 5.196-06 3.70 5.616-05 2.98 1.896+00

4 40x40 3 .936-07 3.72 4.466-06 3.65 1.506+01
80x80 2.956-08 3.74 3.106-07 3.85 1.216+02

160x160 2.386-09 3.63 3.136-08 3.31 1.136+03

 

31

 

 

Table 3.9 Errors and CPU time on the 2D Burger’s equation at t = 0.45 using DG on the

irregular mesh

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Order of

accuracy Grid L1 error L1 order Loo error Loo order CPU (S)
10x10 1.786-03 --- 6.1 16-03 --- 4.226-01
2 20x20 1.756-04 3.35 5.906-04 3.37 3.376+00
40x40 3.296-05 2.41 1.976-04 1.58 2.736+01
80x 80 6. 846-06 2.27 5.066-05 l .96 2.266+02

160x160 1.536-06 2.16 9.946-06 2.35 1.826+03

10x10 9.496-04 --- 4.856-03 --- 9.206-01
3 20x20 7.686-05 3.63 4.626-04 3.39 7.416+00
40x40 7.736-06 3.31 5.116-05 3.18 6.046+01
80x80 8.266-07 3.23 9.776-06 2.39 5.066+02

160x160 9.586-08 3.11 2.126-06 2.20 4.116+03
10x10 7. 186-04 --- 4.876-03 --- 1.916+00

4 20x20 8.806-06 6.35 1.146-04 5.42 1.556+01
40x40 2.636-07 5.06 2.226-05 2.36 1.266+02

80x 80 1.506-08 4.13 2.606-08 9.73 1.2064-03

160x 160 3.636-09 2.05 5.936-09 2.14 9.886-1-03

 

 

 

 

 

 

 

Table 3.10 Errors and CPU time on the 2D Burger’s equation at t = 0.45 using SV on the

irregular mesh

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Order of
accuracy Grid L1 error L1 order Loo error Loo order CPU (8)
10x10 1.996-03 --- 3 . 806-03 --- 2.056-01
2 20x20 4.846-04 2.04 1.536-03 1.31 1.676+00
40x40 1.146-04 2.09 3.926-04 1.97 1.816+01
80x 80 2.876-05 1.99 1.296-04 1.60 1576-1-02
160x160 7.116-06 2.01 3.316-05 1.97 1.216+03
10x10 1.056-03 --- 6.006-03 --- 6.736-01
3 20x20 8.036-05 3.71 6.206-04 3.27 5.856+00
40x40 7.306-06 3.46 3.896-05 3.99 5.826+01
80x80 8.686-07 3.07 7.306-06 2.41 4.536+02
160x160 1.116-07 2.97 1.466-06 2.32 3.616+03
10X 10 4.396-04 --- 3.256-03 --- 1.6 16+00
4 20x20 4.636-06 6.57 7.146-05 5.51 1.436+01
40x40 9.876-08 5.55 1.396-06 5.69 1.236+02
80x 80 1.556-08 2.67 3.566-08 5.28 1.026-1-03
160x160 3 .656-09 2.09 7.146-09 2.32 8.706+03

 

 

 

 

 

 

 

32

 

 

3.2.2 Euler Equations

We consider two—dimensional inviscid flow, and the governing equations are described as

follows,

0_Q +V o F = 0

a:
"p" ' pu ' ' pv '
pu puu + p puv

Q= .f= ,g= ,

pv puv pvv+p
L E _ Lu(E+ P)- _v(E+ p)]

 

 

 

 

 

 

I
Where, ,0 is the density, u and v are the velocity component in xand y directions, pis

the pressure, and E is the total energy. The pressure is related to the total energy by

1
E=—” +—(.0u2+,0v2)
y—l 2

with the ratio of specific heat y being a constant. F = [f , g] is the flux vector.

3.2.2.1 Vortex propagation problem
To compare the numerical accuracy by DG and SV methods, we test the vortex
propagation problem since it has analytical solution. This is an idealized problem for the

Euler equations in 2D, which was used by Shu [56].

The mean flow is {,0, u, v, p} = {1, l, 1, 1}. An isotropic vortex is then added to the mean

flow, i.e., with perturbations in u, v, and temperature T = p/p, and no perturbation in

entropy S = p/p’:

33

where (3r',52')=(x—5,y—5),r2 =i2+y2 , and the vortex strength 8 = 5. In the

numerical simulation, the computational domain is taken to be [0, 10] x [0, 10], with

characteristic inflow and outflow boundary conditions imposed on the boundaries.

It can be readily verified that the Euler equations with the above initial conditions admit
an exact solution that moves with the speed (1, 1) in the diagonal direction. Both the DG
and SV methods were employed to simulate this problem. The numerical simulation was

carried out until I = 10 on the two different grids shown in Figure 3.1.

The errors are computed based on the volume-averaged density on the element or the SV.
Table 3.11 and Table 3.12 present the errors and recorded CPU times of both methods on
the regular mesh, while Table 3.13 and Table 3.14 display the errors and CPU times on
the irregular mesh. Note that the DG method is more consistent in achieving the expected
order of accuracy than the SV method on both the regular and irregular grids. It is
interesting to see that the second-order SV method produced smaller error than the
second-order DG method for the Euler equations. The third and fourth-order DG schemes

appear to have smaller error magnitude than the corresponding SV schemes.

34

Based on the CPU times for the regular mesh, we note for per time step per element that

the DG method takes 43.47, 104.14, and 212.35 #5 for linear, quadratic, and cubic
elements respectively, while SV method spends 27.42, 99.88, and 237.69 [13 for the 2nd,

3rd and 4‘h order schemes. The SV method is faster than the DG method at 2'”d and 3rd

order, but is slightly slower at 4‘h order.

35

Table 3.11 Errors and CPU times for the propagating vortex case at t = 10 using the DG

method on the regular mesh

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Order of
accuracy Grid L1 error L1 order Loo error L0,, order CPU (8)
10x 10x2 7.74e-04 -- 4.02e-03 -- 8.9le+00
20x20x2 1.056-04 2.88 1.126-03 1.84 7.406+01
2 40x40x2 1.526-05 2.79 2736-04 2.04 6126402
80x80x2 2.396-06 2.67 1.216-04 1.17 4.706+03
10x10x2 2.86e-04 -- 2.22e-03 -- 2.11e-r-01
20x20x2 7.546-05 1.92 9.986-04 1.15 l.726+02
3 40x40x2 1.266-05 2.58 1.68604 2.57 l.366+03
80x 80x2 1.146-06 3.47 2.606-05 2.69 l.106+04
lele2 1 .206-04 -- 8. 126-04 -- 4.456+01
20x20x2 6.606-06 4.18 7.146-05 3.51 3.656+02
4 40x40x2 1 .476-07 5.49 3 .426-06 4.38 2.89e+03
80x80x2 3906—09 5.24 2.106-07 4.03 2.316+04

 

 

 

 

 

 

 

Table 3.12 Errors and CPU time for the propagating vortex case at t = 10 using the SV
method on the regular mesh

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Order of
accuracy Grid L1 error L1 order Loo error Loo order CPU (S)
lele2 9.50e-04 -- 6.19e-03 -- 3.77e+00
20x20x2 1.92e-04 2.31 1.61e-03 1.94 2.96e+01
2 40x40x2 4.14e—05 2.21 8.65e-04 0.90 2.38e+02
80x80x2 9.92e—06 2.06 2.96e-O4 1.55 1.91e+03
10x10x2 9.42e-04 -- 5.49e-03 -- 2.19e+01
20x20x2 9.20e-05 3.36 9.91e-04 2.47 1.68e+02
3 40x40x2 9.84e-06 3.22 2.45e-04 2.02 1.30e+03
80x80x2 1.1 1e-06 3.15 3.56e-05 2.78 1.02e+04
10x10x2 1.82e-04 -- 1.20e-03 —- 5.05e+01
20x20x2 1.01e-05 4.17 9036-05 3.73 3.91e+02
4 40x40x2 4.90e-07 4.37 1.01e-05 3.16 3.35e+03
80x 80x2 3.16e-08 3.95 5.8 le-07 4.12 2.46e+04

 

 

 

 

 

 

 

36

 

 

Table 3.13 Errors and CPU times for propagating vortex case at t = 10 using the DG
method on the irregular mesh

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Order of
accuracy Grid L1 error L1 order Loo error Loo order CPU (8)
10x10 6.296-04 -- 3 .326-03 -- 1.026+01
20x20 1.156-04 2.45 9.976-04 1.74 8806-1-01
2 40x40 2.786-05 2.05 3.056-04 1.71 7.076+02
80x80 4.206-06 2.73 8 .936-05 1.77 5 .626+03
10x10 1.276-04 -- 8.886-04 -- 2.416+01
20x20 1.836—05 2.79 3 .266-04 1.45 2.056+02
3 40x40 1.936-06 3.25 9.866-05 1.73 1.6064-03
80x80 1.526-07 3.67 1.616-05 2.61 1 .296-1-04
10x10 4.236-05 -- 3 .696-04 -- 5.126+01
20x20 2.566-06 4.05 6.176-05 2.58 4.946+02
4 40x40 9.806-08 4.71 3.086-06 4.32 3.386+03
80x80 2.676-09 5.20 2.126-07 3.86 2,696+04

 

Table 3.14 Errors and CPU times for the propagating vortex case at t = 10 using the SV
method on the irregular mesh

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Order of
accuracy Grid L1 error L1 order L0,, error Loo order CPU (S)
10x10 1.04e-03 -- 5.456-03 -- 4.38e+00
20x20 2.54e-04 2.03 1.926-03 1.51 3.426-1-01
2 40x40 8.98e-05 1.50 8.726-04 1.14 3.16e+02
80x80 2.346-05 1.94 2.906-04 1.59 2.22e+03
10x10 5.06e-04 -- 3.47e-03 -- 253er
20x20 7.44e-05 2.77 7.25e-04 2.26 l.93e+02
3 40x40 9.75e-06 2.93 1.636-04 2.15 1.51e+03
80x80 1.516-06 2.69 3.37e-05 2.27 1.19e+04
10x10 1.096-04 -- 5.576-04 -— 571er
20x20 7.36e-06 3.89 8766-05 2.67 4.45e+02
4 40x40 3.846-07 4.26 5.756-06 3.93 3.516+03
80x80 2.226-08 4.11 4.00e-07 3.85 2.81e+04

 

 

 

 

 

 

 

37

 

 

3.2.2.2 Double Mach reflection problem

This problem is also a standard test case [71] for high-resolution schemes, and has been
studied extensively by many researchers. The computational domain for this problem is
chosen to be [0, 4] x [0, 1], which is displayed in Figure 3.2. The reflecting wall lies at
the bottom of the computational domain starting from x=1/6. Initially a right-moving
Mach 10 shock is positioned at x=1/6, y=0 and makes a 60° angle with the x-axis. For the
bottom boundary, the exact post-shock condition is imposed for the region from x=0 to
x=l/6 and a solid wall boundary condition is used for the rest. For the top boundary of the
computational domain, the solution is set to describe the exact motion of the Mach 10
shock. The left boundary is set at the exact post-shock condition, while the right

boundary is set as an outflow boundary.

 

computational domain
[0.4114091]

 

 

Exact post-shock
Solid wall boundary

Exact motion of Mach 10 shock
Out now boundary

Figure 3.2 Computational domain and boundaries for double Mach reflection problem
(this figure presented by color)

38

The numerical simulation was carried out until I = 0.2. A mesh refinement study was
carried out on three different grids. The grids are generated from regular Cartesian
meshes by subdividing each Cartesian cell into two triangles. The coarse grid has
25*100*2 triangles, the medium grid 50*188*2 triangles, and the fine grid consists of
120*480*2 triangles. The density contours with 30 equally spaced contour lines from

p = 1.528 to ,0 = 20.863 are shown in Figure 3.3 and Figure 3.4 for the second order DG

scheme and SV scheme. Note that the “blown—up” region was also shown in those

figures.

Note that the SV method has a higher resolution than the DG method for the shock, slip
line and the other finer features near the triple point. The main reason is that the TVD
limiter in the SV method is applied for the sub-cells, but the limiter in the DG method is

applied for the elements (macro SVs).

39

 

 

 

 

°'5 " x-“T’I—‘i --_:—.....-\§: X .4
n I,":" " c: ‘ <
' " (<0
>' 0 ’ r r 1 1 1
x
(a) Coarse grid

    

 

 

.. ,0“ Ed 3

(b) Medium grid

 

 

 

(c) Fine grid

Figure 3.3 Density contours computed with second order DG scheme using a TVD
limiter (30 equally spaced contour lines from ,0 =1.528 to p = 20.863 ).

40

 

 

 

 

 

(c) Fine grid

Figure 3.4 Density contours computed with second order SV scheme using a TVD
lirrriter (30 equally spaced contour lines from p =1.528 to p = 20.863)

41

3.3 Conclusions

We have presented a comparison of the DG and SV methods for the 2D scalar
conservation laws and Euler equations. Generally speaking, both DG and SV method
have achieved the desired order of accuracy. The DG method has a lower error
magnitude than the SV method while SV is faster than DG. In the scalar case, the SV
schemes are consistently faster than the DG schemes of the same order of accuracy for
each residual evaluation. For the Euler equations, the 2“d-order SV scheme is faster than
the 2nd-order DG scheme. However, 3rd and 4th order SV schemes are quite similar to the
corresponding DG schemes in terms of efficiency (<12 % in difference). It is also clear
that the SV method has a higher resolution for discontinuities than the DG method
because of the sub-cell average based data limiting. We also confirm that the SV method
takes less memory and allows larger time steps than the DG method for both the 2D

scalar conservation laws and Euler equations.

42

CHAPTER 4

SV METHOD FOR THE DIFFUSION EQUATION

As a first-step towards extending the SV method to the Navier-Stokes e4quations, the SV
method is extended to and tested for the diffusion equation. In this chapter, we consider
the following 1D diffusion equation

u, =uxx, x6[0,27[] (4.1)

with periodic boundary conditions and initial condition u(x,0) = sin(x).

4.1 Three Formulations
4.1.1 Naive SV Formulation

Directly following the basic formulation described in [64] for the one-dimensional

hyperbolic conservation law, we integrate Eq. (4.1) in control volume C,- which is a

.j r
sub-cell of a spectral volume S,- =[x,-_1/2,x,-+1,2] depicted in Figure 4.1, replace the

flux by a

1____. 1 , 1

33,02 x5312 X1512

 

Figure 4.1 Linear spectral volume

numerical flux and obtain

dfii,j(t)_ 1
dt 11,31.

A

“x

A

i,j+1/2 ‘uxlr,j—1/2)=O'

 

(4.2)

 

43

Since there is no convection term in the diffusion equation, the first derivative is
“naturally” computed by taking a simple average of the derivatives from the two

neighboring CVs, i.e.,

A

1 _
"1 i.j+l/2 :§((ux):j+l/2 +(“INJH/Zl- (4-3)

For time integration, we employ the third order TVD Runge-Kutta method [19]. This
formulation was used to compute a numerical solution for Eq. (4.1) at t =0.7. Two
different grids were used in the simulation. In Figures 4.2 and 4.3, the numerical
solutions with 40 and 320 SVs are compared with the exact solution using linear and
quadratic reconstructions. It seems this formulation leads to a seemingly converged, but

wrong solution. Note that the numerical solutions have an 0(1) error, which does not

decrease with grid refinement. The same phenomenon was reported by Zhang and Shu

[72] for DG.

44

 

1.0 V I V V I Y Y Y Y T Y Vf fir V 0 Y 1’ ‘ 771' fi I I T I 1 r Y I’ I

    

 

 

 

”'3 3 —— Exact Solution 3
3_ 0 Linear SV (40 Cells) 3
0" : 09009:; —--- Linear sv (320 Cells) :
0.4 I -:
02 i -j
=1 0.0c €
.02 E a
-0.4 i
2 9 6,6 i
-0.6 E 9090 i
-0.8 E i
'11,. 14111L111111111111114111111L11111‘
0.0 1.0 2.0 3.0 4.0 5.0 6.0 7.0

X

Figure 4.2 The numerical solutions versus the exact solution using the linear
reconstruction based on the Naive SV formulation for the diffusion equation

 

1.0 h I V V V I V V T Y I 7 r T r r V T T Y7 I r I T I l I V V V I V V V V
p

"-3 — Exact Solution
t o Quadrtic sv (40 Cells)
- - - - Quadratic SV (320 Cells)

LlLLLLAiLLL‘

0.6 ,-
0.4 :
0.2
=1 0.0 c

-0.2 i

-0.4 ’

 

-0.6 :— i
E i
-0.8 f “j
I 43 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 i
0.0 1.0 2.0 3.0 4.0 5.0 6.0 7.0
X

 

 

-1.0

 

Figure 4.3 The numerical solutions versus the exact solution using the quadratic
reconstruction based on the Naive SV formulation for the diffusion equation

45

4.1.2 Local SV Formulation

The second formulation is obtained by mimicking the local discontinuous Galerkin

method [19]. A new variable q is introduced, which is equal to u x, the gradient of u.

Then the diffusion equation becomes the following system
_ = 0
{u’ q" . (4.4)

The spectral volume method is then applied to this system directly. Integrating Eq. (4.4)

over each control volume, we obtain

lama)” 1 . 1
d, "E:(qii,j+1/2’qli,j—1/2):O

< (4.5)

 

A

u

 

_ 1 1
Lqi’j-h- .(ui.j+l/2— li.j—1/2)"O'
1,]

 

The numerical fluxes are chosen as following [57]

+

i.j+I/2 = ”li.j+1/2 (4‘6)

Lil

qii,j+l/2 : qli,j+l/2’ (47)

i.e. we alternatively take the downwind value for u and upwind value for q (we could of
course also take the opposite pattern). Let m be the degree of the reconstruction
polynomial. Numerical solutions are computed at t = 1.0 for the three cases m = 1 (liner
reconstruction), m = 2 (quadratic reconstruction) and m = 3 (cubic reconstruction). The
L1 and Loo errors and numerically observed orders of accuracy are presented in Table

tlr

4.1, from which we note that a (m+1) order of accuracy is achieved for a degree m

polynomial reconstruction.

46

Table 4.1 L1 and Loo errors and orders of accuracy based on the local SV formulation
for the diffusion equation

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Order of
accuracy h L1 error L1 order L,<> error Loo order

272/10 2.356-02 3.616-02
272/20 6.006-03 1.97 9.446-03 1.94
2 272/40 1.516-03 1.99 2.376-03 1.99
(Linear SV) 2 72/80 3.786-04 2.00 5.946-04 2.00
272/160 9.456-05 2.00 1.496-04 2.00

272/320 2.366-05 2.00 3712-05 2.01

272/10 1.156-03 1.78e-03
272/20 1.426-04 3.02 2.226-04 3.00
272/40 1.76e-05 3.01 2.77e-05 3.00
3 272/80 2.206-06 3.00 3.46e-O6 3.00
(Quadratic SV) 272/160 2.756-07 3.00 4.326-07 3.00
272/320 3.446-08 3.00 5.406-08 3.00

272/10 6.996-05 1.076-04
272/20 4.36e-06 4.00 6.82e-O6 3.97
272/40 2.72e-07 4.00 4.27e-07 4.00
4 272/80 1.706-08 4.00 2.66e-08 4.00
(Cubic SV) 272/160 1.05e-09 4.02 1.656-09 4.01
272 /320 5.38e-ll 4.29 8.456-11 4.29

 

47

 

4.1.3 Penalty SV Formulation

In order to remedy the first formulation, Baumann and Oden [8], also Oden, Babuska,
and Baumann [44], and Riviere, Wheeler and Girault [51] introduced a penalty term to
the numerical flux in the DG implementation. However if the formulation of Baumann
and Oden [8] is used directly in the SV method, the penalty term vanishes because the
weighting function is piece-wise constant in the SV method. Therefore the Baumann and
Oden formulation for the SV method is identical to the first formulation. Instead, a
penalty-like term in the following form is added to the numerical flux for the SV method

at the interface. The SV scheme then becomes

 

 

(1‘72“ ‘(0 1

.1 1 1 _

dt — hid- 0‘" i.j+l/2 _l‘xit'.)'-1/2)’0 (4.8)
A _1 + _ 8 + — 49
“4 i,j+1/2 “5((“x)t,j+1/2+(“x)i.j+l/2)+h__f(“lr,j+1/2‘ulr,j+1/2)’ ( - )

1,}
where a is a constant. A Fourier analysis is performed for this formulation in the case of

m = l, and it is found that 6‘ must be one to preserve 2nd order accuracy. Furthermore,
numerical simulations have showed that this formulation can achieve 2nd order accuracy

for linear and quadratic reconstructions, and 4th order accuracy for cubic reconstructions.
Table 4.2 shows the L1 and Leo errors and the numerically observed orders of accuracy

att= 1.0.

48

Table 4.2 L1 and Leo errors and orders of accuracy based on the penalty SV formulation

for the diffusion equation

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Order of
accuracy h L1 error L1 order L0,, error Loo order

272 /10 6.056-03 9.356-03
272 /20 1.51e-03 2.00 2.34e-03 2.00
2 272 /40 3.786-04 2.00 5.92e-04 1.98
(Linear SV) 272 /80 9.46e-05 2.00 1.48e-04 2.00
272 /160 2.366-05 2.00 3.71e-05 2.00
272/320 5.916-06 2.00 9.286-06 2.00

272 /10 2776-03 4.286-03
272 /20 6.77e-04 2.03 1.056-03 2.03
3 272 /40 1.68e—04 2.01 2.636-04 2.00
(Quadratic SV) 272 /80 4.206-05 2.00 6.606-05 1.99
272/160 1.05e-05 2.00 1.65e-05 2.00
2 72/320 2.636-06 2.00 4.13e-06 2.00

2 72/10 6.47e-05 1.00e-04
272 /20 3.996-06 4.02 6.16e-06 4.02
4 272 /40 2.486-07 4.01 3.886-07 3.99
(Cubic SV) 2 72/80 1.55e-08 4.00 2.436-08 4.00
272/160 9.70e-10 4.00 1.52e-09 4.00
272/320 6.40e-11 3.92 1.01e-10 3.91

 

 

4.2 Fourier Analysis
In this analysis, we follow a technique described by Zhang and Shu [72], and focus on the

linear reconstruction only. In this case, a SV is partitioned into two equal CVs, as shown

in Figure 4.1, and CV-averaged mean solutions are 17 17,1 and it" 132' Assuming that the

. . 2n h . .
mesh rs uniform, we have h = h j = N’h 171 = h 132 = 5. The linear reconstructron can be

expressed as
P} (x) = L1 (1)5},1 + 100017132 . (4.10)

with

49

(4.11)

2 l
L1(x)=-Z(x—xj,1/2 +x—xj,5/2)+Z(x—xj,1,2 +x—xj,3/2),

‘4 j-1,2

1
lq(x)=;(x-xj,1/2 +x—xj’3/2). (4.12)
The derivative of pj (x) is constant in Sj , i e ,
, 2 __ 2 _
p,- (x)=-;u,-,1 +7713. (4.13)
All the three SV formulations can be cast in the following form
17- 17-_ 17- 17-
.” =A- __1 1’1 +B- _J’l +C- _j+l’l , (4714)
'4 j,2 “141,2

1
dt u 132
where A, B and C are constant matrices. We seek general solutions of the following form

u(x,t) = tik (t)e'.k“‘r ,
where k is the index of modes (k = 1, 2, ...) representing the wave number, and I

represents the units of the imaginary number. Obviously, the analytical solution for (4.1)

.' — 2 . .
"U k t . In addrtron, we have

 

is u(x,!)ze
P xj,3/2 i
-:-iik (t) J‘e'kgdfi

[171310) = Xj,1/2

ii720) 2 1175/2

_~ .22
hit/20) je d5
4933/2 4

 

Assuming 77 = g — x133”, we obtain

50

 

 

 

 

 

 

_xj,3/2 q — 0 — — 0 -
je df I810? 413/2)” jet/277d”
fill/2 —h/2 -h/2 .
_ _. lkits/2
_ -— e .
xj,5/2 2 ik( . h/2.
- 77+x ,3/2) k
Ie'kgdfi I e 1 d7] Ie' ”d7;
[1933/2 _ r 0 ~ - 0 -
Therefore, the solutions we are looking for can be expressed as
17- (t) 12 (I) 762-
_1’1 = f" 2' 13/2, (4.15)
141-’20) uk,2(t)

where

 

 

 

[2 , 0 7. '
. —-uk(t) [2' ”217;
“12.10) h —/7/2
4220) m.

377,,(7) je'Mdn
lh
0 a
The initial condition can be computed from
, “3,2. _
;' jelxdx
L7 131(0) Xj,1/2
= , (4.16)
17,;2 (0) H.512
71" Ie'xdx
L xj,3/2 j

 

by taking the imaginary part. Note that this analysis depends on the assumption of

uniform mesh size and periodic boundary conditions. Substituting Eq. (4.15) into Eq.

(4.14), we obtain the following evolution equation

51

117.1(1) : G(k, h).[‘:k,1(t)], (4.17)
uk’2(t “1.20)

where the amplification matrix is given by

G(k,h) = (”"1 -A+ B +e"‘" C. (4.18)

The following equation can be used to find the initial condition for (4.17),

u—j

5423/2
Ieixdx
101(0) , 1931/2

=2..e-'kxj.3/2 . . (4.19)
1712,2(0) h xj,5/2
J'eixdx
_4 133/2

 

 

In particular for the low frequency mode k = 1, we have

P —

0 n
1 Ie'xdx
“1.1“” 2 -/7/2
=—- =—- . 4.20)
* h h/2 ih 77/2 (
2712(0) [“2172 e. _1
e

lo -

Let 711 and 7172 be the two eigenvalues of the amplification matrix G(k,h) , V1 and V2 be

l_e—”1/2

 

 

the corresponding eigenvectors of G(k,h) , the general solution of Eq. (4.17) can be

expressed as

[Li/2,10)

_ 7111 7121
1 —a'-e V + -e V. (4.21)
141130)] 1 73 2

By studying the properties of this general solution at the lowest mode (k =1), we can

obtain consistency and convergence results; by investigating the boundedness of the

52

general solution at the high modes (large k), we can establish the stability of the

formulations.

4.2.1

Naive SV Formulation

The naive SV formulation can be expressed as

 

 

 

{417,1 2_
=—“'—11"
dt [,2 J ’
4
(15,32 2_
. dt

ii—
112

2 _ _ _
=—u- ——u- ——u- +—-u-
1,1 1,2 j+1,1 1+l,2
h2 h2 h2 h2

uj--1,2 ’

2

The coefficient matrices A, B, and C are

43. __E_
h2 h2
A:
0 0
. 1

 

 

-11

h2

2

 

1'23

2
h2

2

h2
2
hzy

 

and the amplification matrix G(k,h) is given by

p—

112
G(k,h) =
2

l h2

 

The eigenvalues and eigenvectors of the amplification matrix G(k,h) are

2 _.
e rkh_

ikh +

£1

h2
2

112

711 = ——4— (l — cos(kh)) ,
h2

-—iklr
V]: e ,
l

2

h2
2

112

~ikh + 22
h

ikh _ _2_

h2

22:0,

53

17 +-—2 if
131 1.2
hz

—

 

 

‘

 

(4.22)

(4.23)

(424)

(4.25)

(4.26)

Note that G(k,h) has a zero eigenvalue, which may cause a weak instability for this

semi-discrete system.

We first study the lowest mode, i.e., k =1. Applying the initial condition Eq. (4.20), the

coefficients a and ,8 can be computed as

a = 4(1— cos(h / 2))

ih(e—ih -1) '

’ fl : ih(e"" — 1)

(4.27)

 

 

We therefore have the explicit solution of the SV scheme (4.22) with the initial condition

(4.18). For example

7714(7) =(a-e41’ .24” + 62/12’)-2""‘J'v3”2 (4.28)
Applying a Taylor expansion assuming small h, we obtain the imaginary part of if j,1 (t)
to be

1+ e_2’

 

1m{77j,1(7)} =( )-sin(xj,1)+0(h),

where xj,1=%(xj,1,2+xj’3,2). The solution is about 0.6233sin(xj,1) at t=0.7 ,

which agrees very well with the numerical solution shown in Figure 4.2. We also clearly

see that the scheme is not consistent, i.e. the numerical solution does not converge to the

solution of the PDE (which equals to sin(x)e-t ).

Next we study the stability of the Naive SV formulation by considering the high modes

(large k). When cos(kh) = l , the amplification matrix G(k,h) = 0. Therefore the solution

to Eq. (4.17) remains to be the initial solution given by Eq. (4.21). When cos(kh) ¢ 1, the

54

amplification matrix G(k,h) is diagonalizable with the matrix composed of the

eigenvectors

-khi
R: e 1 , (4.29)
1 1

which has an inverse when cos(kh) at l

-1 1 1 -1

We can obtain explicitly the solution as

l2U“) G(k h): ‘7ka)
,1 = ’ 1. , 4.31
[111(20):] e uk,2 (0) ( )

with

7111’
ecu/2,17): =R- 6 0 -R'1.
0 1

G(k,h) t

The L2 norm of the matrix e can be computed explicitly, which reveals the

stability property. The two eigenvalues of the symmetric matrix (em/“’0 ’ )H (eG(k"’) ’)

(where AH 2 (KY , AT denotes the transpose of A, and A is the conjugate of A) are

found to be

 

A1 ={(1—27)2 +2222—(1—22)\[(1—22)2 +2/10’}/0',

 

A2 =1<1—/7>2 +ua+(1—u>J(1-u>2 +2/zarla,

55

with ,u=e’11‘ , a=l—cos(kh). Then the [72 norm of emk’h)’ is given by

 

eG(k,h) t

H 6602,11)!" = Jmax(IA1|,|A2|) . If is uniformly bounded with respect. to k and

 

 

 

 

h, (4.32) is said to be stable. However, if we take a = h2 / t , we can easily get

 

/\1={1—2-e‘4 +e‘8 -(l—e_4)\/l+e-8 -2-e‘4 +2-e‘4h2/7 +e‘4h2/7}-2/h2,

 

A2 = {1—2-e‘4 +278 +(1—.«.»‘4)\/1+e'8 -2-e-4 +2.e‘4h2 /7 +2'4h2 /t}-t/h2,
and obviously

"2ka"I = 0(1/17),

 

 

which is unbounded when h ——> 0. Therefore system (4.17) for the naive SV scheme is

unstable.

4.2.2 Local SV Formulation

The local SV formulation can be written as

f

 

 

 

duj,1_1_ 5_ 12- 4- 3- 1—
217 —;,—u,-_1,1+h—2uj_1,2 ”7.7"“ + hz “/2 + hz “1+IJ‘ hz “1+”
] . (4.32)
dl7' 2 6 6 2
1.2 — - - —
1 7. = 7.2“41’73923—29m"Writ

The corresponding coefficient matrices A, B, and C are

_L 3_ -12 LL 3___L
h2 h2 h2 h2 h2 h2
A: , B: , C: . m3»
0 () LL -1i __._ji
L _ _ h2 h2_ _h2 h2_

 

 

 

 

 

 

The amplification matrix, its eigenvalues and eigenvectors are

56

”Le—ikhWLieikIuE _5_e—ikh__l_eikh 4 -

 

 

G(k,h) = (4.34)
£621.12 +_2_ __2_eikh_£
L h2 172 172 172 -
711 = —%, 7172 = ——23(1-cos(kh)) (4.35)
h h
_ ikh .
V - ——5+e- V — 1937”“) 436
1

Clearly both eigenvalues are real and negative. To study the accuracy and consistency,
we again examine the lowest mode k =1. By applying the initial condition Eq. (4.20),

the coefficients a: and ,6 in Eq. (4.21) are found to be

a: 7. 27. (4.37)
(1+14el +e ' )ih

 

 

,8 = . . (4.38)
(1+ 142'” + eZ'h )ih
We thus have the following explicit solution for scheme (4.32)
_ ih . - .
271(2) =(a-241’ —5—+—e.— + 6842910?“ + 1)) «“123”- . (4.39)
J, 1+ 3e"' 2

Applying a Taylor expansion assuming small h, the imaginary part of 17 131 (t) is found to
be
Im{rTJ-,l (2)} = sin(xj,1)e_t + 0072).

Clearly, the numerical solution converges to the exact solution with second order

accuracy.

57

The matrix composed of the eigenvectors of G(k,h) is

—5+eikh 1

—'kh
_ ——,— —(e ' +1)
R - 1+ Belkh 2 9 (4'40)
1 1
with its inverse given by
R_1= 1 -1—3-e""" (1+e—ikh)~(l+3-eikh)/2 (441)
7 +COS(kh) 1+3-e‘k" 5—eikh ' '

In order to prove the stability of scheme (4.32), it is sufficient to show that the Z72 norms

of R and RT1 are uniformly bounded with respect to k and h since both eigenvalues of

G(k,h) are negative. In fact, the eigenvalues of RH R are

A 1 = (51 + 10 . cos(kh) + 3 - (cos(kh))2 -

 

J641— 716 - cos(kh) + 30 - (cos(kh))2 + 36 - (cos(kh))3 + 9 - (cos(kh))4 )/(2(10 + 6cos(kh)))

A2 = (51+10-cos(kh) +3-(cos(kh))2 +

 

J641— 716 - cos(kh) + 30 . (cos(kh))2 + 36 - (cos(kh))3 + 9 . (cos(kh))4 )/(2(10 + 6cos(kh))),

and the eigenvalues of (RTl)H (R—l) are

2771 = (1225 + 1330 - cos(kh) + 452 - (cos(kh))2 + 62- (cos(kh))3 + 3-(cos(12h))4

 

- (7 + cos(kh))2 - J 745 — 288 - cos(kh) - 306 - (cos(kh))2 + 96- (cos(kh))3 + 9 - (cos(kh))4 )/ D

(02 = (1225 +1330-cos(kh) + 452 - (cos(kh))2 + 62 - (cos(kh))3 + 3 - (cos(kh))4

 

+ (7 + cos(kh))z . J 745 - 288 - cos(kh) — 306 - (cos(kh))2 + 96- (cos(kh))3 + 9 - (cos(kh))4 )/ o

58

where D = 2 - (2401 + 1372 - cos(kh) + 294 - (cos(kh))2 + 28 2 (cos(kh))3 + (cos(kh))4)

 

 

Hence "R" = JmadilHAZD and "R71" = Jmaxdwlllalz I) . It is easy to see that both “R"

and [IR-l" are uniformly bounded with respect to kh. Thus the stability of scheme (4.32)

is established.

4.2.3 Penalty SV Formulation

Finally we turn to the analysis of the third formulation, for which we obtain the following

scheme based on the linear reconstruction

 

d17-
},l 2 _ 2 _ 2 _ 2 _
=—l-£u-_ ——1—3£u-_ -— 1+3Eu- +—1+£u~
dt h2( ) J 1,1 h2( ) 1 1,2 h2( ) 1,1 h2( ) 1,2
4 (4.42)
417m 2

 

 

2 2 2
=—1+817- ——1+3£u- ——l-3cu‘- +—l—£u’-

k

Choosing 8 =1, the scheme (4.42) reduces to

 

 

 

rdfijJ 4 _ 8 _ 4 _
dt = pal-4,2 —h—2-uj,1+;2—uj,2
4 . (4.43)
21171-2 4 _ 8 _ 4 _
dt =h—zuj,l ‘22—“);2 +Fuj+hl

The coefficient matrices are

— I1 I— —l _ .-

 

 

 

 

4 8 4
° 7.7 72 772' 0
A: B: (4.44)
0 0 i -1 i.
_ L h?” hzd Lhz

The amplification matrix G(k,h) is

59

 

 

r _i i—ikh+4
h2 h2 h2
(sanh)= . (445)
_1_2u. :1. __§_
_h2 h2 h2 _

 

 

The two eigenvalues of G(k,h) are
8 8
21 = —-—2 (1 + cos(kh / 2)) , 712 = ——7 (l — cos(kh/ 2)) . (4.46)
h h

Clearly both eigenvalues are real and negative. The corresponding eigenvectors are

_ e-ikh/Z e—ikh/Z
m: , n: . MM)
1 1

The coefficients a and ,6 in Eq. (4.21) are

277/2 _
aza 6:39—7—11 MA&
1h
We thus have the explicit solutions of scheme (4.42). For example
l—l-j,1(t) : (a.e/111 . (_e—ih/2) + fl ezizt oe~lh/2) _ eixj,3/2 . (4.49)

Using a Taylor expansion, we obtain the imaginary part of 171310) to be

Im{17j,1(t)} = smut/(1)2" + 0(h2).
Clearly, the scheme is consistent and second order accurate. Note that we may take
8 >0(h) for consistency, but we only have Im{u‘j,1(t)} =sin(xj,1)e"’ +0(h) ifs #1,

which is why we suggest 6‘ = 1.

The matrix composed of the eigenvectors (4.47) of G(k,h) is

60

_ —rk/7/2 -ikh/2
R=[ e 6’ ] (4.50)
l 1

with its inverse given by

'kh/2
_1 l-e‘ l
R =— . . 4.51
2!:etkh/2 1] ( )

Finally the IQ norms of R and R”1 can be computed. They take the following form

"R"=\/2- and "R-1||=\/2/2.

It is clear that both "R“ and "RT1 II are uniformly bounded with respect to kh. Thus the

stability of scheme (4.42) is established.

4.3 Conclusions

Three different formulations of the spectral volume method are presented for the
diffusion equation. Numerical tests and analysis are performed for these formulations.
We have found that both the local SV and penalty-like SV formulations are stable and
consistent while the naive SV formulation is neither consistent nor stable. Numerical

results agree well with the analysis. It appears that the local SV formulation achieves

(m +1)’h order accuracy with a degree m polynomial reconstruction, while the penalty

th It

SV formulation achieves (m +1) order accuracy if m is odd, but 271' order accuracy if

m is even.

61

Since the local SV formulation achieves the highest accuracy for a given polynomial
reconstruction, we will extended it to the 1D/2D convection-diffusion equations in

Chapter 5, to the Navier-Stokes equations in Chapter 6.

62

CHAPTER 5

EXTENSION OF SV METHOD TO CONVECTION-DIFFUSION
EQUATIONS

5.1 Formulation for 1D Convection-Diffusion Equation

Consider the following 1D convection-diffusion equation with fixed boundary conditions

$14-53; f(u)—i g[u,%%)=0 in (22,17)ch (5.1a)
u(a,t) = ua , u(b,t) = ”b (5.1b)

Bu . . . . .
where f (u) and g(u,-5—) rnrght be lrner or nonlinear contrnuous functrons, ua and ab are
x

constant.

Introducing an auxiliary variable v = g_u , equation (5.1) becomes
x

v = a_u in (a,b) (5.2a)
ax

Bu 8 a .

E+$f(u)——a—;g(u,v)—O 1n (a,b) (5.2b)

u(a,t) = ua , u(b,t) = “b- (5.2c)

Integrating equations (5.2a-5.2b) on each control volume C,- which is a sub-cell of a

.j’
spectral volume S,- = [xi_1/2,xi+1/2] depicted in Figure 5.1, and replacing the flux by

the numerical flux, we obtain the following equations for til-d- , 172', j ,

63

 

 

 

 

 

' _ 1 . .
V13} _E(uii,j+l/2 _u|i,j-1/2)_ 0 (5-3-1)
i
dil- - 1 1
1.} “Fl—(fl _fl )T—L(§i' ° 1/2_§i' ' 1/2):O (53-2)
(1! hm- i.j+l/2 i,j—1/2 hi.j 171+ 1.1-
' l ' Linear sv
79,172 x1372 11,572
I I I ' Quadratic SV
79,172 32,372 x1312 15,772
' ' l l I Cubic SV
X117 2 323/2 X1572 XI)?” 32972

Figure 5.1 One-dimensional spectral volumes

The numerical fluxes are chosen as follows
At internal interfaces:

7. +

“i.j+1/2 =“2.j+1/2 (in (5'3“)

64

. _ f(u|;j+l/2) ifaf/au>0

. _. 3.. =g(27.. ,v|.‘. )(in(5.3.2))
lj+l/2 f(ulsz/Z) ifdf/du<0 rJ+1/2 ll,]+1/2 l,j+1/2

OT

i,j+1/2

A

u (in (53.1))

ii,j+1/2 =“1

l _ f(u|;j+l,2) ifaf/au>0

. — + (in (53.2))
l7‘+1/2 f(ulinIZ) ifdf/au<0

1 — +
gliJ-H/z = g(uii,j+1/2’vli,j+1/2)

— +
Where’ 17ii,j+1/2 ‘ (“12,141/2 +u|i,j+1/2)/2

At x = a:

2|a =ua (in (53.1))

fl“ = f(ua) 2|, = 2(u.,v|;) (in (53.2))
At x = b:

2|b=ub (in the(5.3.1))
flb =f<ub> 2|, =2<ub.v|;) (in (5.3.2»

If the boundary is not fixed, for example, 31 = va is given instead of u x: = ua ,
x=a

then the calculation of numerical flux at x = a becomes

= 24; (in (53.1))

ila = 7:04;) 2|, = 2(u|;‘.v,) (in (5.3.2»

The flux calculation at internal interfaces is the same as before.

65

5.2 Formulation for 2D Linear Convection-Diffusion Equation

Consider a model second-order convection-diffusion equation with proper boundary

conditions
u,+V-(6-u)—Vo(AVu)=0 11122ch (5.4a)
u = f on To (5.4b)
(AVu) 0 n = g on I‘N (5.40)

411 412

] is a matrix, and n = (nx,ny)is the unit
A21 422

where ,6 = (7317/32) is a vector, A =[

outward normal of the boundary of the computational domain 9..

By introducing auxiliary unknowns v = Bu / ax , and w = Bu / 0y , Eq. (5.4) becomes

vzau/ax in o: R2 (5.5a)
w=au/ay in 12: R2 (5.5b)
u,+Vo(,6-u)-V-(A(v,w)T)=0 in 52:122 (5.5c)
u = f on PD (5.521)
(A(v,w)T)-n = g on 1“,, (5.5e)

Integrating equations (5.5a-5.50) on each control volume C i, J- , we obtain the following

integral form for the CV-averaged mean for 27,3]- , Vi, j 7872‘, j ,

K
_ 1
7,3]. _V—(Z Iu-ndi) =0 (5.6a)
i,j r=1Ar

66

K
W . L(Zju-nyd71)=0 (5.6b)

 

1,] Vial r=lAr
daij 1 K K T
+——{Z I[(fl1u)x-nx+(fl2u)y~ny]dA-Z [(716,247) )ondA}=0 (5.6c)
‘1‘ V2.7 7:11, 7:11,

Note that we use the same reconstruction for the auxiliary unknowns v,w as that for the

original unknown u . We use the following upwind numerical flux

u-nx z (u ~nx)R (5.7)

u-ny z (u-ny)R (5.8)
((flruix-n.+(flzu> -n )L.iffl°n>0

(flu>.-n,+(flu> '27 = ,. ,, (5.9)
l 2 , , {((fllu)x'"x+(flzu)y'"y)R7if/3'n<0

(41v. w)T ) on - ((A(v.w>T)- 22)" (5.10)

and for boundary faces

. , 1‘
u-nxz u nxLon D (5.11)
(u-nx) ,onFN

~ u-ny,on FD
u ny ~{(u.ny)l‘,on TN (512)
«Aux-Murry) -n )LJffl-wo
(flu)x'nx+(fl Ll) 'n z y y (513)
l 2 y y 1((fllu)x .nx +(fl2u)y 'ny)Rll:ffl.n<O
T L
(A(V,w)T)onz{((A(v’w; ).n) ,OnrD (5.14)
(A(v,w) )Onzg,onFN

67

5.3 Accuracy study

Extensive accuracy studies were carried out for both 1D and 2D convection and diffusion
equations. In 1D, both linear and non-linear equations are employed, and problems with
exact solutions are designed to test the local spectral volume (LSV) approach. These

accuracy studies are presented next.

5.3.1 1D linear Convection-Diffusion Equation

The following linear equation is solved with SV schemes of various orders

it, +c-ux —,u-uxx =0
subject to the initial condition of u(x,0) = sin(x) and periodic boundary condition. The

2
computational domain is [-7r,7r]. The exact solution is u(x,t) = e": I“ sin(rr(x—ct)).
The numerical simulation was carried out until I =1 . The L1 and Loo errors are presented
in Table 5.1 for c = 1.0, ,u = 1.0. Note that the LSV approach is capable of achieving the

optimum orders of accuracy in all cases.

Table 5.1 L1 and Leo errors and orders of accuracy on 1D linear convection-diffusion
equation (at t = 1.0)

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Order of
accuracy nTCell L1 error L1 order L0,, error L0,, order
10 7.606-03 -- 9.976-03 --
20 2.076-03 1.88 2.896-03 1.79
2 40 5.466-04 1.92 7.866-04 1.88
80 1.406-04 1.96 2.046-04 1.95
10 3.836-04 -- 5 .376-04 --
20 4.556-05 3.07 6.646-05 3.02
3 40 5.576-06 3.03 8.166-06 3.02
80 6.896-07 3.02 1.016-06 3.01
10 1.246-05 -- 1.856-05 --
4 20 7.926-07 3.97 1.176-06 3.98
40 5.016-08 3.98 7.436-08 3.98
80 3.156-09 3.99 4.696-09 3.99

 

 

68

5.3.2 1D Viscous Burger’s Equation

Consider
u, +u-ux —,qux =0, [1:01, x6 (0,1)

with the following initial condition u(x,0) = —tanh(2i) and boundary condition
,u

u(0,t) = O, u(l,t) = —tanh[—l—].
2!!

The problem has the following exact solution
u(x,t) = —tanh[i].
2!!

The simulation is conducted until t = 1.0 with various SV schemes. The L1 and Loo

errors are presented in Table 5.2. Note that the LSV approach is again capable of

achieving the optimum orders of accuracy in all cases.

5.3.3 1D Fully Nonlinear Case

Further, we consider the following fully nonlinear equation

u,+u-ux—%(u-ux)x=0 xE(0,1)

with initial and boundary conditions as u(x,0) = e"; u(0,t) =1, u(l,t) = e. The problem

has the following exact solution u(x,t) = ex .
The simulation was conducted until t = 1.0 with various SV schemes. The L1 and L0,,

errors are presented in Table 5.3. Note that the LSV approach is again capable of

achieving the optimum order of accuracy in all cases.

69

Table 5.2 L1 and Loo errors and orders of accuracy on 1D viscous Burger’s equation

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

(at t = 1.0)
Order of
accuracy nTCell L1 error L1 order L0,, error Loo order
10 2.126-03 -- 1.226-02 --
20 6.166-04 1.78 4.306-03 1.50
2 40 1.736-04 1.83 1.336-03 1.69
80 4.766-05 1.86 3 .666-04 1.86
10 4.846-04 -- 3.436-03 --
20 7.626-05 2.67 5.256-04 2.71
3 4O 1.046-05 2.87 7066-05 2.89
80 1.366-06 2.93 9.126-06 2.95
10 3.176-05 -- 2.166-04 --
4 20 1.496-06 4.41 9.536-06 4.50
40 7.316-08 4.35 8.626-07 3.47
80 4.226-09 4.1 1 6.206-08 3.80

 

Table 5.3 L1 and L0,, errors and orders of accuracy on the 1D fully nonlinear
convection-diffusion equation (I = 1.0)

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Order of
accuracy nTCell L1 error Ll order Loo error Loo order
10 9.526-04 -- 2.16-03 --
20 2.536-04 1.91 6.126-04 1.78
2 40 6.526-05 1.96 1.656—04 1.89
80 1.656-05 1.98 4.276-05 1.95
10 9.326-06 -- 2.926-05 --
20 1.226-06 2.93 3.996-06 2.87
3 40 1.566-07 2.97 5.216-07 2.94
80 1.976-08 2.99 6.656-08 2.97
10 6.306-08 -- 2.376-07 --
4 20 4.076-09 3.95 1.646-08 3.85
40 2.596-10 3.97 1.076-09 3.94
80 1.636-11 3.99 6.816-11 3.97

 

70

 

 

5.3.4 2D Linear Convection-Diffusion Equation
We also tested the LSV method on a 2D linear equation written as
u, +c(ux +uy)—,u(uxx +uyy) = 0, (x,y)E(—1,1)X(—1,l);c =1, ,u = 0.01
with the initial condition u(x, y,0) = sin(7r(x+ y)) and periodic boundary condition. The

2
‘2” 4“ sin(rz(x+ y - 2ct)) . The recorded L1 and L0,, errors

exact solution is u(x, y,t) = e
in Table 5.4 again show that the LSV approach is capable of achieving the optimum

orders of accuracy in all cases.

Table 5.4 L1 and Loo errors and orders of accuracy on the 2D linear convection-
diffusion equation (I = 1.0)

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Order of
accuracy Grid L1 error L1 order Loo error L0,, order
10x10x2 3.736-02 -- 5.356-02 --
20x20x2 8606—03 2.12 1.286-02 2.06
2 40x40x2 2.166-03 1.99 3.256-03 1.98
80x80x2 5.366-04 2.01 8.116-04 2.00
10x10x2 2.796-03 -- 3.916-03 --
20x20x2 3406—04 3.04 5.106-04 2.94
3 40x40x2 4.166-05 3.03 6.316-05 3.01
80x80x2 5.216-06 3.00 7.976-06 2.98
lele2 4.14e-05 -- 5.72e-05 --
20x20x2 2.476-06 4.07 3 .376-06 4.09
4 40x40x2 1476-07 4.07 1.956-07 4.1 1
80x80x2 8.456-09 4.12 1.146-08 4.10

 

 

5.4 Conclusions

We’ve extended the local spectral volume method to the 1D and 2D scalar convection-
diffusion equations. Accuracy studies with 1D linear, viscous Burger’s, fully nonlinear
and 2D linear cases have been carried out, and the order of accuracy claim has been
numerically verified, i.e. the local spectral volume method achieve 2“, 3rd and 4th order

of accuracy for linear, quadratic and cubic reconstructions.

71

CHAPTER 6
EXTENSION OF SV METHOD TO THE NAVIER-STOKES
EQUATIONS

6.1 Formulation for Navier-Stokes Equations

We consider the two-dimensional Navier-Stokes equations written in conservation form

%?—+VoF8(Q)—Vof,(o,vo)=0 (6.1a)

in the computational domain 52 subject to suitable initial and boundary conditions. The

conservative variables Q and the Cartesian components fe (Q) and ge(Q) of the

inviscid flux vector Fe (Q) are given by

 

 

 

 

p l P“ i . W l
2
Q= p“. f.<Q)=<"“ fl"). 2.<Q)=< ‘2‘” l. (64b)
10V puv 10" +17
E ,ME+pi |uE+py

Here ,0 is the density, u and v are the velocity components in x and y directions, p is the

pressure, and E is the total energy. The pressure is related to the total energy by

E=—”—I+%p(u2+v2) (6.1c)

with ratio of specific heats y, which is taken to be 1.4 in all the simulations in this

research. The Cartesian components fv(Q,VQ) and gv(Q,VQ) of the viscous flux

vector Ev (Q,VQ) are given by

72

r O 1

2ux +71(ux +vy)

fv(Q7VQ)=.u'i v, +77, i (6.1d)
Cp

)+-———Tx
P

r

u[2ux +49% +vy)]+v(vx +uy

 

 

0
vx +uy
g,(Q,VQ)=,7-J 2vy+71(ux+vy) (6.1e)

 

 

C
P
ku(vx +uy)+v[2vy +71.(uJr +vy )]+-B,—Ty

where ,u is the dynamic viscosity, C p is the specific heat at constant pressure, P, is the

Prandle number, T is the temperature and using the Stokes hypothesis, 71 = —2/ 3.

Integrating Eq. (6.1a) in C,- we obtain the following integral equation for the CV-

 

.j ’
averaged mean
dQ- . K _.
"J + 1 Z [F(Q,VQ)-ndA=0 (6.2)
(It Vi’j r=1Ar

where Q,j is the vector of the CV-averaged conservative variables in C ,- K is the

.j’

number of faces in C and A, represents the r-th face of C,-

l,j ’ ,j 9
f(QNQ): fie(Q)_13v(Q,VQ) _ We treat the gradient of the conservative variables
VQ = G(Q) as auxiliary unknowns of the Navier-Stokes equations following [7,10,19],

which are therefore reformulated as the following coupled system for the unknowns

Gand Q,

73

G—VQ=O mso

%%+Vofie(Q)-V-F,(Q,G)=o (63b)

Integrating system (6.3) over each control volume C we obtain the following integral

2,),

form for the CV-averaged mean G,- j and Q,- j ,

 

K
.,.-27’—.<219-nd4>=0 (6'4”
19.] f=1Ar
d5! , 1 K _. K _.
7.1 + (Z IFe(Q).ndA_ Z IFV(Q,G)ondA)=0. (64b)
dt V571 r=1Ar r=lAr

Note that at the SV boundaries, reconstructed solutions Q and G are not continuous. It is
necessary to substitute those fluxes with interface numerical fluxes. In this research, we

use the following numerical fluxes for internal interfaces:

QoanLon (6.53)
17;, (Q) . n z Roe flux splitting (6.5b)
fi,(Q,o)-nzi,('Q‘,GR)-n (6.56)

where Q =(QL +QR)/ 2. Combining the Gauss quadrature formulae and numerical

fluxes, all the integrals appearing in equations (6.4a-6.4b) can be evaluated. An important
issue is that the reconstructions for the auxiliary variable G have the same structure as the
one employed for the original conservative variable Q. For boundary faces, we borrow

the idea of the reference [7], but the form is different.

74

For inviscid flux:

— —

pVon

Film-n: ”“V’M’m‘ (6.6)

vaon+pny
_.o(E+p)l7-n,

 

 

On both inviscid and no-slip surfaces, the normal velocity V 0n vanishes. The flux is
equal to the pressure contribution of the inviscid flux function in the normal direction to
the surface, with the pressure p being taken form the internal boundary state.

At inflow (outflow), E}, (Q) on = Ee(Qb) on , Qb is computed by imposing the available
data and the Riemann invariant associated to outgoing characteristics.

For the flux in the auxiliary equation:

Q o n ==: Qb o n
where Qb has the same value as that for the inviscid flux.
For the viscous boundary flux:

0n the no-slip surface, Ev (Q,G)0n = Fv(Qb,Gb)on, where Qbis that for the inviscid

flux; Gb = GL if there are no boundary conditions on VQ on. When such conditions are

instead prescribed, the value of Gb is modified accordingly.

At inflow (outflow), Ev (Q,G)0n = Ev (QL,GL)On

75

6.2 Numerical Experiments
6.2.1 The Couette Flow

To verify the formal accuracy of the local SV method on Navier-Stokes equations, we
consider the Couette flow between two parallel walls. The lower is stationary with

temperature To and the upper is moving at speed of U with temperature T1. The distance

between the two walls is H.

The steady analytic solution is

U
u=— , v=0
Hy

 

2
y #-U y y
T=T +—- T—T + -——- 1——
0 H (1 0) 2k H( H)

P
= constant, = —
P 70 R . T

The parameters were chosen as U =1.0, H = 2.0, To = 0.8, T1 = 0.85, ,1! = 0.01 . The

formulation described above was tested with this problem on the two dimensional domain

[0,2]x[0,2] . The computational grid is displayed as Figure 6.1.

Convergence test: When the simulation was started with the following initial condition:
it =0,v=0,p=l,T=l.

Figure 6.2 shows that the numerical solution converges to the steady analytical solution.

Error estimation 1: The simulation was started with the steady solution. We record the

L1 and L0,, errors of the numerical solution away from the steady solution until the

76

residuals approach to the machine zero. Tables 6.1 and 6.2 demonstrate the L1 and L0,,

errors for various reconstructions with mesh refinements in density and temperature in

this sense.

Error estimation 2: The simulation was also stared with the steady solution. We record
the L1 and L0,, errors of the numerical solution away from the steady solution until the

physical time t=1.0 for all reconstructions and grids. The numerical results are

displayed in Tables 6.3 and 6.4.

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Figure 6.1 Computational grid for Couette flow case

77

u—velocity Profile
( order = 3, slmpleBxBDTF )

 

      

 

 

2.0

1.5 -
—-- initial
O—Ostep=100

in 1.0 - H step=500 -

A—Astep=1000
Hstep=2000
D—D step=5000
k-A step=8000

05 - Hatep=10000 r
‘i——l" 8kp=m
-— steadysohltion 1

00 4 n l l L

 

 

' —0.2 0.0 0.2 0.4 0.6 0.8 1.0 1.2

Temperature Profile

( order = 3, simple8x8.DTF)
2.0 . 1 . , . .

 

- - - initial

1.5 F 0—0 step = 100
H step = 500

, A—-A step = 1000
H step = 2000
D—Cl step = 5000

>’ 1'0 _ k-A step=8000
H step = 10000

> -I— -+- step = 20000
— steady sohltion

I
r
l
I
I
t
l
O

 

 

0.0 * 1
0.6 0.65 0.7

 

 

Figure 6.2 Convergence trend of numerical solution to the steady exact solution in
Couette flow test

78

Table 6.1 L1 and L0,, errors and orders of accuracy on Couette flow

(Density in Estimation 1)

 

Accuracy of

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

order Grid L1 error L1 order L0,, error Loo order
2x2x2 1.52796-02 --- 2.42636-02 ---
2 4x4x2 4.00456-03 1.93 8.68406-03 1.48
8x8x2 1.07676-03 1.90 3.25616-03 1.42
16x 16x2 3.00546-04 1.84 1.32376-03 1.30
2x2x2 1.54936-04 --- 6.32436-04 ---
3 4x4x2 3.83106-05 2.02 1.27016-04 2.32
8x8x2 7.59356-06 2.33 2.18126-05 2.54
16x 16x2 1.30646-06 2.54 3.46976-06 2.65
2x2x2 4.9086e-05 -- 9.4826e-05 --—
4 4x4x2 4.01506-06 3.61 8.12476-06 3.54
8x8x2 3.25356-07 3.63 7.67646-07 3.40
16x16x2 2.64446-08 3.62 1.07716-07 2.83

 

Table 6.2 L1 and L0,, errors and orders of accuracy on Couette flow

(Temperature in Estimation 1)

 

Accuracy of

 

 

 

 

 

 

 

 

 

 

 

 

 

order Grid L1 error L1 order Loo error Loo order
2x2x2 3.21546-03 --- 9.1 1846-03 ---
2 4x4x2 8.61936-04 1.90 2.28786-03 1.99
8x8x2 2.13496-04 2.01 8.95396-04 1.35
16x 16x2 5.15776-05 2.05 4.61196-04 0.96
2x2x2 4.06876-05 --- 1 .77 876-04 ---
3 4x4x2 5.47956-06 2.89 1.92636-05 3.21
8x8x2 8.22046-07 2.74 2.52006-06 2.93
16x16x2 1.13576-07 2.86 4.20906—07 2.58
2x2x2 5.51086-06 --- 1.96156-05 ---
4 4x4x2 3.95746-07 3.80 1.49416-06 3.71
8x8x2 2.69446-08 3.88 1.32526-07 3.49
16x16x2 2.10466-09 3.68 1.15506-07 0.2

 

 

 

 

 

 

79

 

 

Table 6.3 L1 and Loo errors and orders of accuracy on Couette flow

(Density in Estimation 2)

 

Accuracy of

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

order Grid L1 error L1 order L0,, error L0,, order
2x2x2 5 .65066-03 --- l .22766-02 ---
2 4x4x2 1.79486-03 1.65 4.35196-03 1.50
8x8x2 4.71206-04 1.93 1.41616-03 1.62
16x16x2 l. 15856-04 2.02 5.46566-04 1.37
2x2x2 9.18196-05 --- 4.36846-04 ---
3 4x4x2 1.73556-05 2.40 9.27206-05 2.24
8x8x2 2.82736-06 2.62 1.74546-05 2.41
16x16x2 4.08876-07 2.79 2.81766-06 2.63
2x2x2 2.09696-05 --- 5.46576-05 ---
4 4x4x2 1.55426-06 3.75 5.16306-06 3.40
8x8x2 1.04836-07 3.89 4.52256-07 3.51
16x 16x2 6.78746-09 3.95 4.47486-08 3.34

 

Table 6.4 L1 and Leo errors and orders of accuracy on Couette flow

(Temperature in Estimation 2)

 

Accuracy of

 

 

 

 

 

 

 

 

 

 

 

 

 

order Grid L1 error L1 order L0,, error Loo order
2x2x2 1.73036-03 --- 5.09136-03 ---
2 4x4x2 5.26096-04 1.72 1.69076-03 1.59
8x8x2 1.45766-04 1.85 4.28636-04 1.98
16x 16x2 3.74756-05 1.96 1.55686-04 1.46
2x2x2 2.63046-05 --- 7.99466-05 ---
3 4x4x2 3.56496-06 2.88 1.36436-05 2.55
8x8x2 5.24626-07 2.76 1.84146-06 2.89
16x 16x2 7.56206-08 2.79 2.64956-07 2.80
2x2x2 4.78286-06 --- 1.59846-05 ---
4 4x4x2 3 .37936-07 3.82 1.27036-06 3.65
8x8x2 2.39366-08 3.82 9.67066-08 3.72
16x 16x2 1.59906-09 3.90 6.69686-09 3.85

 

 

 

 

 

 

80

 

 

6.2.2 Laminar Flow along a Flat Plate

We consider the laminar flow on an adiabatic flat pate characterized by a free stream
Mach number Ma = 0.3 and by a Reynolds number based on the free stream condition
and on the plate length Re = 10000. The length of the plate is set to L :10 as shown in

Figure 6.3. At x = 1.0, based on the Blasius solution, the thickness of the boundary layer

 

 

5 is
6| =10=5 ”' =5-—5———=0.05
x . poo 'uoo x:1.0 1|R6|x=130

Therefore, the size of computational domain in y-direction is chosen to be 20 times of the
boundary layer thickness at x = 1.0 such that the flow at the top boundary is nearly
inviscid. The range of computational domain in x-direction is [-l, 1]. Figure 6.3 shows
the computational domain. The red point (0, 0) indicates the leading edge of the flat plate.
Adiabatic wall boundary is used on the plate surface. From (-1, 0) to (0, 0), a symmetry
boundary condition is applied. At the inlet, the free stream condition is imposed. At the

top and exit boundaries, the static pressure is fixed.

 

Computational Domain
I-l.l]X[0.ll

 

 

 

 

-_ .‘ll'h'u. -.-. -«v:c ‘-.1-._L_‘ "LIA—HF

a flat plate

Figure 6.3 Computational domains for the flat plate case
(this figure presented by color)

81

Since the singularity around the leading edge of the flat plate, we employ cluster meshes
both in x-direction and y-direction. The computations have been performed on three
triangular meshes, coarse mesh (208 cells) with 8 cells along the plate, medium mesh
(832 cells) with 16 cells along the plate, and fine mesh (3328 cells) with 32 cells along

the plate.

In y-direction of the coarse mesh, we set two cells within the boundary layer, and do a
proper transmission to the rest. In x-direction we set the proper clustering factor such that
the grids almost have the same length in x and y direction around the leading edge of the

plate. Figure 6.4 shows three different meshes used in our numerical simulations.

82

 

(a) Coarse mesh

 

(b) Medium mesh

 

(c) Fine mesh

Figure 6.4 Meshes for the flat plate case

83

Figure 6.5 shows the profiles based on the numerical solution and Blasius solution.

Where u is the velocity in x-direction, and y"= = y- ’M .
,u-x

 

 
  
  
   

 

 

 

 

 

a: 5 — Blasius solutlon
>0 —-B—-Coarse mesh
4 - __e_- Medium mesh
3 —-q—-Flne mesh
2
1
O l
0 0 2 0.4 0 6 0 8 1

Figure 6.5a u—velocity profiles for the flat plate case with linear SV

84

 

Blasius solution
— -E|— — Coarse mesh
_ .0... — Medium mesh
-— +— - Fine mesh

 

 

 

 

 

u/U

Figure 6.5b u-velocity profiles for the flat plate case with quadratic SV

 

Blasius solution
_ _a_ _ Coarse mesh
_ Medium mesh

 

 

 

 

 

 

Figure 6.5c u-velocity profiles for the flat plate case with cubic SV

85

Furthermore, we compute skin friction coefficient C f by both numerical solution and the

well-known Blasius formula for the C f distribution along a flat plate in the case of

incompressible flow in Figure 6.6. The computed results show a good agreement with the
Blasius solution. Also, the residual history plot in Figure 6.7 shows that we achieve the
steady state solution numerically.

Numerical solution:

 

2'
C f — w
L . 2
2 p... U...
Blasius solution:
Cf : 0.664
Re x
u poo -U°o -x .
where, shear stress 1w = ,u-— , Reynolds number Re I =————, U 00 IS the
FD ,u

speed of the free stream, and x is the distance away from the leading edge of the plate.

86

0.2 -

 

 

0.15 -
Blasiussoiution
—o-E|—-— Coersemesh
-..0_.- Mediummesh
230.1 ‘ -..p..- Finemesh

 

 

 

 

 

Figure 6.63 Skin friction coefficient along a flat plate with linear SV

87

0.2

0.05

Figure 6.6b Skin friction coefficient along a flat plate with quadratic SV

 

Blasius solution
_._a_ . .. Coarse mesh

-..Q_... Medium mesh
-..p..- Finemesh

 

 

 

 

 

 

88

0.2

 

0.15 "
Biasiussolutlon
-..g..-Coerse mesh
.._ . -..0..-Medium mesh
0 0.1 . -.+.- Flnemesh

 

 

 

0.05 '

 

 

 

Figure 6.6c Skin friction coefficient along a flat plate with cubic SV

89

    

 

 

:10'25

.9

*5 -

3 10'35-

O' E

0 :
£1045-

: E

.E Z
2105?

3 E
"510'“;

To E

3 c
2104:-

m E

d) E
I"I110”:- . 1

0 1 2 3

Iteration step/100000

Figure 6.7 Residual history of continuity equation with quadratic SV

90

6.2.3 Subsonic Flow over a Circular Cylinder

When a fluid flows over an isolated cylindrical solid barrier and the Reynolds number is
great than about 50, vortices are shed on the down streamside. The vortices trail behind
the cylinder in two rolls, alternatively from the top or the bottom of the cylinder. This
vortex trail is called the Von Karman vortex street or Karmn street after Von Karman’s
1912 mathematical description of the phenomenon. Since then, many numerical and
experimental studies have focused on the dynamics of vortex street formation in the near
wake. Measured value of the Strouhal number of the vortex shedding frequency can be
found in Reference [70]. Spectral solutions of cylinder wake flows include those of

[21,28,38]. A recent detailed review of the problem can be found in Reference [9].

fv'Dc

 

Where, the Strouhal number S = , Dcis the diameter of the cylinder, fv is the

frequency of vortex shedding, and V is the flow velocity. To demonstrate the capability
of LSV method in dealing with complex geometry (with curve wall effect), and complex

flow properties, we choose this as one of our test cases.

The computations have been performed on the following unstructured grid shown in

Figure 6.8. The center of cylinder is located at the origin and its diameter is equal to 1.

The computational domain are (-10, 16) x (-10, 10).

91

 

 

 

 

 

 

 

 

 

 

(a)

\

 
 
  
  
    
  
  
  
 

4A 4 <5 A
4) «ENE ‘Hfi‘w a:
Faaeawt‘ewaéiéram
4) ‘9‘ "‘V 4) vuv
,Qzexa'eEgAsxattets'a‘exaau
422219;?“ ‘ ttggga uuuua
N‘ 5’55?
‘Kxfi‘x‘t‘x .rasséss‘flk‘fifi
‘9‘ews- mtwsxv55‘V‘V“
\»§J5aa16v¢w> «my; Q'AVAV
$fi'EEaWtéifi‘yvyfifi9s
‘Vfi'fizwléflhyéé'é'fln

(b)

Figure 6.8 Grid for the case of subsonic flow over a circular cylinder
(a) global grid, (b) grid near the cylinder

92

We consider a subsonic flow over an adiabatic circular cylinder at a small angle, free
stream Mach numberMa = 0.2 , Reynolds number based on the free stream condition and
on the cylinder diameter Re =75. We fix pressure P on the right boundary of the

rectangular and fix every thing on the rest boundary of the rectangular.

After a sufficiently long time, the effects of the initial condition propagate out of the
computational domain, and the periodic shedding of vortices is observed. Instantaneous
contours of the Mach number, entropy, and vorticity showing the Von Karman vortex
street generated by the cylinder are presented in Figures 6.9, 6.10 and 6.11. The periodic
nature of the flow is show in Figure 6.12, which plot as a function of time the pressure
calculated at three different locations, i.e. at points (1,1), (5,1), and (10,1). The period of
the oscillations corresponds to a Strouhal number of 0.151, which agrees with Reference

[63, 65, 66] very well.

 

Figure 6.9 Instantaneous Mach contours for Ma = 0.2 flow
over a circular cylinder at Re = 75

93

 

Figure 6.10 Instantaneous entropy contours for Ma = 0.2 flow
over a circular cylinder at Re = 75

 

Figure 6.11 Instantaneous vorticity conto'urs for Ma = 0.2 flow
over a circular cylinder at Re = 75

94

Pressure

Pressure

0.99

0.988

0.986

0.984

0.982

 

0.98
o I
Ti me

0.99

0.988

0.986

0.984

 

0.982
130 140 150

Tlme

Figure 6.12a Pressure history at (x, y) = (1,1) for Ma = 0.2 flow
over a circular cylinder at Re = 75

95

Pressure

Pressure

‘LOO4

‘1002

0998

0996

0994

0992

099

1004

1002

0998

0996

0994

0992

099

130

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

: l T ' j j ' z 7
‘
z n s 1 s z s s
at“; 3 “L .§.i;up§ . i -
~- 1- 1 1 ~ ‘ 1' 111 1 1~
@- L--+-.. g. A ..J i 4.; -1
3 1
2 1- 1 1 -~
~- 1 «11- f
,.. .é. J ...é... 1 .é... 1 —.
3 2 f
i i '
’1 - i
50 100 150
Time
. . r T I .r . T .
—- 5 I -—
_. 3 3 .1
_.-.§,__ ; . , ;.- , , , , , =, -
1 ‘1 %

 

 

140
TI me

150

Figure 6.12b. Pressure history at (x, y) = (5,1) for Ma 2 0.2 flow

96

over a circular cylinder at Re = 75

1 .002

0.998

Pressure

0.996

0.994

 

0.992

Time

 

1.002

 

 

0.998 \ ‘ I

0.996 — ..... .. \ I

Pressure

 

 

0.994 \. \V/ V

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

0.992 =
130 140 150

Time

 

Figure 6.12c Pressure history at (x, y) = (10,1) for Ma = 0.2 flow
over a circular cylinder at Re = 75

97

6.2.4 Laminar Subsonic Flow around the NACA0012 Airfoil
This test has been considered as validation test cases for shock capturing Navier-Stokes

codes in a GAMM work shop, and is very well documented in the literature [7 , 41, 43, 47,

and 73].

The computations have been performed on the relative coarse grid shown in Figure 6.13,
which is an unstructured triangulation of a 64x16 O-grid. The larger number refers to the
number of cells distributed along the airfoil surface and smaller one to the number of
cells in the radial direction. The grid extents about 20 chords away form the airfoil. The

computations have been performed using linear, quadratic and cubic SVs

1

I;

AA”
Z‘Wmézliia

      
    
  

 

    

\\ Vb vmv
s igltiaatz§5§§$mt§z§ffihi ,
Ea‘fiesss'm AV (;
sssg ‘ V ' -1‘.1-$Y2?§§§§v:fz¢Afl
Ar 7" 5:" z‘iig$%1$~

            

1
‘\
‘A
11'

1,.
i1
11

A
AV
4 .
5’1.
1'11"
1
‘1 l
‘1‘
I“

a aaa%:;§‘€; ¢«> \V‘

mar; amuwexsfi'fifiae» «»

hvfifiiifimflg€ififi\$
auuuwnfltbҤ\

\
s.

s

Figure 6.13 Computational grid for the NACA0012 test case

98

We run a laminar subsonic flow at an angle of attack a=0°, free stream Mach number
Ma =05 , and Reynolds number Re=5000. In the test, the wall is adiabatic. The
- Reynolds number is near the upper limit for steady laminar flow. A distinguishing feature
of this test case is the separation region of the flow occurring near the trailing edge,
which causes the formation of a small recirculation bubble that extend in the near-wake

region of the airfoil.

Figure 6-14 shows the Mach isolines computed with linear, quadratic and cubic SVs. It is
obvious that the solution is getting smoother and smoother with the increasing of the

order of polynomial reconstruction.

Figures 6.15 and 6.16 show the pressure coefficient (C p) and skin friction coefficient
(C f) distributions along the airfoil computed with linear, quadratic, and cubic SVs,

where

7w

1 . 2
21)... U...

C =____p_p°° Cf:

P 1 2
Table 6.5 reports the pressure part of the drag coefficient (C d, p) and stress part of the

drag coefficient (C d,v ) computed With linear, quadratic and cubic SVs, respectively. The

comparison of the values of the drag coefficients computed with cubic SVs and those
obtained by other authors with both structured and unstructured solvers [7, 41, 43, 47] is

given in Table 6.6, where

99

it’d/i

§rwdA
‘1—.2._ -

Cd _——————
V 1pm-Ui-A

9

Cd,p

A is a characteristic area of the object. In addition the residual converge history is

displayed in Figure 6.17, which demonstrates that we have reached the steady solution.

 

Figure 6.14a Mach isolines around the NACA0012 airfoil with linear
sv (Re = 5000, Ma = 0.5, a = 0°)

 

Figure 6.14b Mach isolines around the NACA0012 airfoil with quadratic
SV(Re = 5000, Ma = 0.5, a = 0°)

101

 

Figure 6.14c Mach isolines around the NACA0012 airfoil with cubic
sv (Re = 5000, Ma = 0.5, a = 0° )

102

 

 

 

 

 

 

 

Figure 6.15 Pressure coefficient distribution along the N ACA0012
airfoil (Re = 5000, Ma = 0.5, a = 0° )

103

 

 

 

 

 

 

012 ‘0" o 2ndorderscheme
1%“ a 3rdorderscheme
d“ 4thorderscheme
0.06
0 '3 5:7 = ‘-
5 O . ; ‘z- 2 523‘?“ ' '
-0.06 l
V‘f
1'. '
-0.12 '
_ ' 1 r . l . r 1 r . l
0'180 0.2 0.4 0.6 0.8 1

Figure 6.16 Skin friction coefficient distribution along the NACA0012 airfoil
(Re = 5000, Ma = 0.5, a' = 0°)

104

IO“
.3
°.

to

ty equat

10"

10"

Residual of oontinui

_L

C
.1.
o

 

 

l I l I l l l I l I l
0 1 2 3 4 5

Iteration step/5000

Figure 6.17 Residual history of continuity equation with linear reconstruction
for NACA0012 test case

Table 6.5 Pressure Part (C d, p ) and Viscous Part (C d“, ) of the Drag Coefficient for the

 

 

 

 

 

NACA0012 Airfoil
SV—type Accuracy C d, p C d 3,
Linear 2“d 2.325e-02 3.311e-02
Quadratic 3rd 2.246e-02 3.275e-02
Cubic 4‘h 2.231e-02 3.302e-02

 

 

 

 

105

 

Table 6.6 Comparisons of Pressure Part (Cd,p ) and Viscous Part (CdJ, ) of the Drag
Coefficient for the NACA0012 Airfoil Computed by the Present Method with Various

Structured and Unstructured Solvers

 

 

 

 

 

 

 

 

 

 

 

 

 

Method Grid Degree C d, p C d",
of freedom

Cubic SV 64 x16 10240 0.02231 0.03302
Cubic DG [36] 64 x16 10240 0.02208 0.03303
Triangle scheme from Ref. [43] 320x64 20480 0.0229 0.0332
Cell-center scheme form Ref. [41] 320x64 20480 0.0219 0.0337
Cell-vertex scheme form Ref. [47] 256x64 16384 0.0227 0.0327
Cell-center scheme form Ref. [47] 256x64 16384 0.02256 0.03301
Cell-center scheme form Ref. [47] 512x128 65536 0.02235 0.03299

 

 

6.3 Conclusions

The SV method is successfully extended to the Navier-Stokes equations by following a
mixed formulation named the local discontinuous Galerkin approach originally
developed for the DG method. The approach, which is named the local SV (LSV)
approach, has been tested extensively for 1D and 2D convection-diffusion equations
using a serious of accuracy studies. All tests indicated that the formulation is capable of

achieving the formal optimum order of accuracy in both L1 and Loo norms. The LSV

approach is then implemented and tested for the Navier-Stokes equations, and was able to
achieve the formal order of accuracy for the compressible Couette flow problem. Also,
we have test on more complex problems such as laminar viscous flow along a flat plate,
subsonic flow over a circular cylinder, and laminar viscous flow around NACA0012
airfoil. The case of laminar flow over a flat plate was simulated successfully with good
agreement with the Blasius solution. The numerical results based on LSV approach, for

the case of subsonic flow over a circular cylinder, match with those in experiments and

106

some other method. In the test case of NACA0012 airfoil, the physical solutions have
been achieved, and the numerical solution is getting smoother and smoother with the

increasing of degree of reconstruction polynomials.

107

CHAPTER 7
SUMMARY AND FUTURE WORK

7.1 Summary

Two major research areas have been undertaken in this study: the evaluation of the DG
and SV methods, and the extension of the SV method to the Navier-Stokes equations. In
the first research area, the DG and SV methods have been evaluated for the scalar
conservation laws and the Euler equations in both 1D and 2D. The overall performance
of both methods in terms of the efficiency, accuracy and memory requirement has been
evaluated. In the second research area, several different algorithms have been suggested
and tested for second-order derivatives. Fourier analysis was employed to analyze the
accuracy, consistency and stability of the proposed algorithms for the 1D heat equation.
The analysis has been numerically verified. Then the best performing algorithm, the local
SV approach, has been extended to multi-dimensional scalar convection and diffusion
equation and to the Navier—Stokes equations. Extensive numerical tests have been carried

out to test the overall performance.

Generally speaking, both the DG and SV methods are capable of achieving the formal
order of accuracy, i.e. they both achieve 2'”, 3rd, and 4th order accuracy for the
corresponding linear, quadratic and cubic reconstructions respectively, while the DG
method usually has a lower error magnitude and takes more memory. In the scalar case,
the SV schemes are consistently faster than the DG schemes of the same order of
accuracy for each residual evaluation. For the Euler equations, the 2“d-order SV scheme

is faster than the 2"d-order DG scheme. However, 3rd and 4'h order SV schemes are quite

108

similar to the corresponding DG schemes in terms of efficiency (<12 % in difference). It
is also clear that the SV method has a higher resolution for discontinuities than the DG
method because of the sub-cell average based data limiting. We also confirm that the SV
method takes less memory and allows larger time steps than the DG method for both the

2D scalar conservation laws and the Euler equations.

In order to test the SV method for second-order derivative terms, the 1D heat equation
was first employed to evaluate the SV method in Chapter 4, in which three formulations
are presented, i.e. Naive SV formulation, Local SV formulation and Penalty SV
formulation. A Fourier type analysis has been performed and it has been proved that the
local SV method and penalty SV method can achieve the formal order of accuracy for
linear reconstructions on pure diffusion equation. Second, the local SV method has been
applied to scalar convection-diffusion equations in Chapter 5 based on the optimal
accuracy it has achieved for the 1D heat equation. Extensive accuracy studies were
carried out for both 1D and 2D convection-diffusion equations. The numerical results
have demonstrated that the local SV method can achieve the optimal order of accuracy
for convection-diffusion equation. Finally, the local SV method has been implemented
for Navier—Stokes equations in Chapter 6. We have successfully solved several viscous
flow problems such as the Couette flow, laminar flow along a flat plate, unsteady
subsonic flow over a circular cylinder and laminar subsonic flow around the NACA0012
airfoil, by using local SV method. For the Couette flow problem, the numerical solution
converged to the steady analytical solution and achieved the formal order of accuracy.
For the case of laminar flow along a flat plate, it was demonstrate that the numerical

solution agrees well with the Blasius solution. For the case of subsonic flow over a

109

circular cylinder, the numerical experiment has shown that an unsteady, stable vortex
street formed behind the cylinder at a Reynolds number of Re = 75 . The Strouhal
number of the vortex shedding frequency based on the local SV method matches those
based experimental measurement and a staggered-grid multidomain spectral method. For
the NACA0012 airfoil case, the numerical solution was getting smoother and smoother
with increasing order of reconstructions. In addition, the drag based on the local SV
method, due to pressure and shear stress, was quite similar to that computed with several
2"d order finite volume methods. However, the total number of degrees of freedom is

much less than those used in the 2ud order finite volume methods.

In summary, the SV method has some unique properties comparing with other high-order
methods such as the DG method, and is capable of achieving hi gh-order accuracy for the
Navier-Stokes equations. Accurate numerical results can be computed with much coarser
meshes than those used in second-order finite volume methods. For steady state
computations, however, the explicit high-order SV schemes are still not competitive to
implicit second-order finite volume method because of the time step limit and the slow
convergence of the method. Much more efficient time integration algorithms must be

developed for these hi gh-order methods to be used routinely in engineering design.

7.2 Future Work
Although the feasibility of the SV method for the Navier—Stokes equations has been

successfully demonstrated in this study, major obstacles still remain to apply the method

110

for “real world” engineering flow problems. In order to achieve that goal, the following
activities are planned:
0 Test the SV method for more complex viscous flow problems involving more
complex geometries, higher Reynolds numbers;
0 Develop an efficient time marching algorithm for the SV method. Two
approaches are possible: implicit solution approach and multi-grid approach;
0 Parallelize the SV method on distributed memory machines using domain
decomposition and message passing;
0 Finally the method must be extended to three-dimensions since all real flow

problems are 3D.

lll

APPENDIX A

GAUSS QUADRATURE FORMULAS

We consider integrals and quadrature rules of the form

I = [mad/1 -—= Zwifu’.) (11.1)
A

i=1
where, A is integral domain, W,- are the quadrature rule’s weights and P,- are the
evaluation points, i =1, 2, ..., n. The integration rule (A.l) is called exact to order q if it
is exact when f is any polynomial of degree q or less. In References [1,22], we can see

the detail.
A-l. 1D Quadrature Rule

Considering the one-dimensional quadrature rules on the canonical [-1,l] element
1 n
I = [f(§)d§ = Ewe.) (A2)
_1 i=1
Gaussian quadrature is preferred for numerical integration because they have fewer
evaluation points for a given order. Actually, if we choose n evaluation points, then we
can get an integration rule be exact to order 211-]. The flowing Table A.l gives the
evaluation points 5,, i =1, 2, ..., n, which are the roots of the Legendre polynomial of

degree n. The weights Wi , i =1, 2, ..., n , called Christoffel weights are also shown for n

ranging from 1 to 6.

112

Table A.l Christoffel weights W, and roots 6,- ,i = 1, 2,. ..,n , for Legendre polynomials

of degree 1 to 6

 

 

 

 

 

 

 

 

 

i 51' Wt
0.00000 00000 00000 2.00000 00000 00000
0.57735 02691 89626 1.00000 00000 00000
0.00000 00000 00000 0.88888 88888 88889
0.77459 66692 41483 0.55555 55555 55556
0.33998 10435 84856 0.65214 51548 62546
0.86113 63115 94053 0.34785 48451 37454
0.00000 00000 00000 0.56888 88888 88889
0.53846 93101 05683 0.47862 86704 99366
0.90617 98459 38664 0.23692 68850 56189
0.23861 91860 83197 0.46791 39345 72691
0.66120 93864 66265 0.36076 15730 48139
0.93246 95142 03152 0.17132 44923 79170

 

 

Example. Consider evaluating the integral

b
1 = [ g(x)dx

a

 

 

b—a
2

1:

 

b-a
2

[2:

 

-a

2

 

1
[f(odé ..
-l

by Gauss quadrature. Let us transform the integral to [-1,1] using the mapping

 

 

2W1 f(gi)
i=1

where , we could find Wi and 4‘,- in Table A.l for a given n ranging from 1 to 6, and

f (91') = g(x,~ ), i = 1, 2, ...n . 4f,- and x,- are related each other by relation (A.3). So we get

211/1800)
i=1

 

A-2. 2D Quadrature Rule
Two-dimensional integrals on triangles are conveniently expressed in terms of triangular

coordinates as

jjfrx. mm 2 11.2mm," .53. 5;)
£2 i=1

e
where ({1' , {£410 are the triangular coordinates of evaluation point i and Ae is the area

of triangle e. The relation between the coordinate of (x, y) and the triangular coordinate is
described as follows:

P;

A .,

P:
Let (x,y),(x,~,y,-) denote the triangular coordinate of points P,P,-,i =l,2,3 , then the

triangular coordinates ({1, {2 , {3 ) of point P can be described as

areao APP P l
{1: f 2 3 = .1 x2 y2, (A.4.l)
area 0f AP1P2P3 2A8

 

 

1‘3 Y3
APP P 1 x y
area 0
2 = f 3 1 = '1 X3 Y3 . (A.4.2)
area of AP1P2P3 2A8
x1 1’1
1 x y

areao APPP l
93 = f 1 2 = '1 x1 y1 . (A43)
area of AP1P2P3 2A8
9‘2 Y2

 

114

Note that there is an identity {I + {2 + {3 E 1. Where A, denotes the area of AP1P2P3,

| * | denotes a determinant. Therefore (A.4.l-A.4.3) give a relation between the

coordinate of (x, y) and the coordinate of (41,4343). We also can express (x, y) by

using ({1,§2,§3)as

x=x3+(x1—x3)'§1+(x2—x3)'§2 (Ail)
y=)’3+(}’1“Y3)'§1+(Y2—Y3)°{2 (A-5-2)
We list some quadrature rules in Table A.2. A multiplication factor M indicates the
number of permutations associated with an evaluation point having a weight W. The

factor p indicates the order of the quadrature rule, i.e. the quadrature formula is exact

integration for any polynomial of degree p or less. The error between the exact

integration and quadrature formula is 0(hp+l)where h is the maximum edge of the

triangle.

Table A.2 Weights and evaluation points for integration on triangles

 

W1

{1

{2

 

1 .000000000000000

0.333333333333333

0.333333333333333

 

0.333333333333333

0.666666666666667

0.166666666666667

 

Abbv— 3

-0.562500000000000
0.520833333333333

0.333333333333333
0.600000000000000

0.333333333333333
0.200000000000000

wig—“o

 

0.109951743655322
0.223381589678011

0.816847572980459
0.108103018168070

0.091576213509771
0.445948490915965

 

0.225000000000000
0.125939180544827
0.132394152788506

0.333333333333333
0.797426985353087
0.059715871789770

0.333333333333333
0.101286507323456
0.470142064105115

 

12

0.050844906370207
0.116786275726379
0.082851075618374

0.873821971016996
0.501426509658179
0.636502499121399

0.063089014491502
0.249286745170910
0.310352451033785

 

 

13

-0. 149570044467670
0.175615257433204
0.053347235608839
0.0771 13760890257

0.333333333333333
0.479308067841923
0.869739794195568
0.638444188569809

 

0.333333333333333
0.260345966079038
0.065130102902216
0.312865496004875

 

C‘beJr-‘GUJWWUJHUJUJWU—‘Wv— z

 

 

 

 

115

 

APPENDIX B
PUBLICATIONS ASSOCIATED WITH THIS WORK
Sun, Yuzhi and Wang, Z.J., “Evaluation of Discontinuous Galerkin and Spectral Volume
Methods for Scalar and System Conservation Laws on Unstructured Grids,” International
Journal for Numerical Methods in Fluids, Vol. 45, pp819-838, 2004.
Sun, Yuzhi and Wang, Z.J., “Formulations and Analysis of the Spectral Volume Method
for the Diffusion Equation,” Communications in Numerical Methods in Engineering, Vol.

20, Pp927-937, 2004.

Sun, Yuzhi and Wang, Z.J., “Extension of Spectral Volume Method to the Navier-Stokes
Equations on Unstructured Grids,” to be submitted.

Sun, Yuzhi and Wang, Z.J., “Evaluation of Discontinuous Galerkin and Spectral Volume
Methods for Conservation Laws on Unstructured Grids,” AIAA-2003-0253.

Sun, Yuzhi and Wang, Z.J., “Evaluation of Discontinuous Galerkin and Spectral Volume
Methods for 2D Euler Equations on Unstructured Grids,” AIAA-2003-3680.

Sun, Yuzhi and Wang, Z.J., “High-Order Spectral Volume Method for the Navier—Stokes
Equations on Unstructured Grids,” AIAA-2004-2133.

116

10.

11.

12.

BIBLIOGRAPHY

Abromowitz, M. and Stegun I. A., Handbook of Mathematical Functions, volume 55
of Applied Mathematics Series. National Bureau of Standards, Gathersburg, 1964.

Arnold, D.N., An interior penalty finite element method with discontinuous
elements, SIAM J. Numer. Anal. , 19, 742-760 (1982).

Barth, T.J., Numerical aspects of computing viscous high Reynolds number flows
on unstructured meshes, AIAA 29th Aerospace sciences meeting, AIAA 91-0721,
January 1991.

Barth, T.J. and Frederickson, P.O., High order solution of the Euler equations on
unstructured grids using quadratic reconstruction, AIAA 90-0013, January 1990.

Barth, T.J. and Jesperson, DC, The design and application of upwind schemes on
unstructured meshes, AIAA 27th Aerospace sciences meeting, AIAA 89-0366,
January 1989.

Bassi, F. and Rebay, 8., Hi gh-order accurate discontinuous finite element solution of
the 2D Euler equations, J. Comput. Phys, 138, 251(1997).

Bassi F. and Rebay, S., A High-order accurate discontinuous finite element method
for the numerical solution of the compressible Navier-Stokes equations, J. Comput.
Phys. 131, 267-279 (1997).

Baumann, GE. and Oden, IT, A discontinuous hp finite element method for
convection-diffusion problems. Comput. Methods Appl. Mech. Engrg. 175, 311-341
(1999).

Beaudan, P. and Moin, P., Numerical experiments on the flow past a circular
cylinder at sub-critical Reynolds Number, Technical Report, TF-62, Department of
Mech. E., Stanford University, 1994.

Castilo, Paul, An optimal estimate for the local discontinuous Galerkin method, in
Discontinuous Galerkin Methods, B. Cockbum, G.E. Kamiadakis, and C. -W. Shu,
Lecture Notes in Computational Science and Engineering, 11, Springer, 285-290
(2000).

Cockbum, B., Devising discontinuous Galerkin methods for non-linear hyperbolic
conservation laws, J. Comput. Appl. Math., 128,187-204 (2001). Numerical analysis
2000, Vol. VII, Partial differential equations.

Cockbum, B., Gremaud, P. —A., and Yang, J.X., A priori error estimates for
nonlinear scalar conservation laws, in hyperbolic problems: theory numerics,

117

13.

14.

15.

16.

17.

18.

19.

20.

21.

22.

23.

applications, Vol. I (Ziirich, 1998), vol. 129 of Inemat. Ser. Numer. Math.,
Birkh'auser, Basel, 1999, pp.167-l76.

Cockbum, B., Hon, 8. and Shu, C.-W., TVB Runge-Kutta local projection
discontinuous Galerkin finite element method for conservation laws IV: the

multidimensional case, Mathematics of Computation 54, 545-581(l990).

Cockbum, B., Kamiadakis, G. E., and C. -W. Shu, The development of
discontinuous Galerkin methods, in Discontinuous Galerkin methods (Newport, RI,
1999), vol. 11 of Lee. Notes Comput. Sci Eng., Springer, Berlin, 2000, pp.3-50.

Cockbum, B., Liu, S. -Y. and Shu, C. -W., TVB Runge-Kutta local projection
discontinuous Galerkin finite element method for conservation laws III: one-
dimensional systems, J. Comput. Phys. 84,90-113 (1989).

Cockbum, B. and Shu, C. -W., TVB Runge-Kutta local projection discontinuous
Galerkin finite element method for conservation laws II: general framework,
Mathematics of Computation 52, 411-435 (1989).

Cockbum, B. and Shu, C.-W., The Runge-Kutta local projection P1-discontinuous-
Galerkin finite element method for scalar conservation laws, RAIRO Modél. Math.
Anal. Numér. , 25, 337-361 (1991).

Cockbum, B. and Shu, C. -W., The Runge-Kutta discontinuous Galerkin method for
conservation laws. V. Multidimensional systems, J. Comput. Phys, 141, 199-224
(1998)

Cockbum, B. and Shu, C. -W., The local discontinuous Galerkin method for time-
dependent convection diffusion system. SIAM J. Numer. Anal. , 35, 2440-
2463(1998).

Colella, P. and Woodward, P., The piecewise parabolic method for gas-dynamical
simulations, J. Comput. Phys. 54, (1984).

Don, W. —S. and Gottlieb, D., Spectral simulation of unsteady compressible flow
past a circular cylinder, Comput. Methods Appl. Mech. Eng. 80, 39 (1990).

Dunavant, D.A., High degree efficient symmetrical Gaussian quadrature rules for
the triangle. International Journal of Numerical Methods in Engineering, 21, 1129-
1 148(1985).

Falkkk, R. S. and Richter, G.R., Local error estimates for a finite element method

for hyperbolic and convection-diffusion equations, SIAM J. Numer. Anal, 29,730-
754. (1992).

118

24.

25.

26.

27.

28.

29.

30.

31.

32.

33.

34.

35.

36.

Godunov, S.K., A finite-difference method for the numerical computation of
discontinuous solutions of the equations of fluid dynamics, Mat. Sb. 47,271(1959).

Harten, A., High resolution schemes for hyperbolic conservation laws, J. of Comput.
Phys. 49, 357-393(1983).

Haren, A., Engquist, B, Osher, S., and Chakravarthy, S., Uniformly high order
essentially non-oscillatory schemes III, J. Comput. Phys. 71, 231 (1987).

Hartmann, R. and Houston, P., Adaptive discontinuous Galerkin finite element
methods for nonlinear hyperbolic conservation laws, SIAM J. Sci. Comput., 24, 979-
1004 (2002).

HesthavenJ. S., A stable penalty method for the compressible Navier-Stokes
equations. ii. One dimensional domain decomposition schemes, SIAM J. Sci.
Comput. 18,658 (1997)

Houston, P., Schwab, C., and Sfili, E., Discontinuous hp-finite element methods for
convection-diffusion-reaction problems, SIAM J. Numer. Anal., 39, 2133-2163.
(2002)

Hu, C. and Shu, C.-W., A discontinuous Galerkin finite element method for
Hamilton-Jacobi equations, SIAM J. Sci. Comput., 21, 666-690 (1999).

Hu, C. and Shu, C.-W., Weighted essentially non-oscillatory schemes on triangular
meshes, J. Comput. Phys, 150,97 (1999).

Jameson, A. and Baker, T., Improvements to the aircraft Euler method, AIAA 25‘h
Aerospace sciences meeting, AIAA 87-0452, January 1987.

Jameson, A. and Mavriplis, D., Finite volume solution in the two-dimensional Euler
equations on a regular triangular mesh, AIAA 23rd Aerospace sciences meeting,
AIAA 85-0435, January 1985.

Jameson, A., Baker, T., and Weatherill, N., Calculation of inviscid transonic flow
over a complete aircraft, AIAA 24‘h Aerospace sciences meeting, AIAA 86-0103,
January 1986.

Johnson, C., N'avert, U. and PITKARANTA, J ., Finite element methods for linear
hyperbolic problems, Comput. Methods Appl. Mech. Engrg. , 45, 285-312 (1984).

Johnson, C., Navert, U. and PlTKARANT A, J., An analysis of the discontinuous
Galerkin method for scalar hyperbolic equation, Math. Comp. , 46,1-26 (1986).

119

37.

38.

39.

40.

41.

42.

43.

44.

45.

46.

47.

48.

49.

50

Karakashian, O. and Makridakis, C., A space-time finite element method for the
nonlinear Schrdinger equation: the discontinuous Galerkin method, Math. Comp.,
67, 479-499 (1998).

Kopriva, D. A,. A staggered-grid multidomain spectral method for the compressible
Navier—Stokes equations. J. Comput. Phys, 143, 125-158 (1998).

Lesaint, P. and Raviart, P. A., On a finite element method to solve the neutron
transport equation, in Mathematical Aspects of Finite Elements in Partial
Differential Equations, edited by C. de. Boor (Academic Press, New York, 1974),
P.89.

Liu, M. —S., Mass flux schemes and connection to shock instability, J Comput. Phys.
160, 623-648 (2000).

Martinelli, L., Calculations of viscous flows with a multigrid method, Ph.D thesis,
Dept. of Mechanical and Aerospace Engineering, Princeton University, 1987.

Mavriplis,D. and Jameson, A., Multigrid solution of the Navier-Stokes equations on
triangular meshes, AIAA Journal, vol. 28, pp. 1415-1425, January 1986.

Mavriplis, D. J. and Jameson, A., Multigrid solution of the Navier-Stokes equations
on triangular meshes, AIAA J. 28, No. 8,1415(1990).

Oden, J.T., Babuska, Ivo and Baumann, C.E., A discontinuous hp finite element
method for diffusion problems, J. Comput. Phys, 146, 491-519(l998).

Osher, S., Riemann solvers, the entropy condition, and difference approximations,
SIAM J. on Numerical Analysis 21, 217-235 (1984).

Osher, S. and Chakravarthy, S., Upwind schemes and boundary conditions with
applications to Euler equation in general geometries, J. Comput. Phys. 50 447-
481(1983).

Radespiel, R. and Swanson, R. C., An investigation of cell centered and cell vertex
multi grid schemes for the Navier-Stokes equations, AIAA paper No. 89-0453, 1989.

Reed, W. H. and Hill, T. R., Triangular mesh methods for the neutron transport
equation, Los Alamos Scientific Laboratory Report LA-UR-73-479,
1973(unpublished).

Richter, G. R., An optimal-order error estimate for the discontinuous Galerkin
method, Math. Comp., 50,75-88 (1988).

Richter, G. R., The discontinuous Galerkin method with diffusion, Math. Comp., 58,
631-643 (1992).

120

51.

52.

53.

54.

55.

56.

57.

58.

59.

60

61.

62.

63.

Riviere, B., Wheeler, M. and Girault, V., Apriori error estimates for finite element
methods based on discontinuous approximation spaces for elliptic problems, SIAM
J. Numer. Anal., , 39, 902-931 (2001).

Roe, P.L., Approximate Riemann solvers, parameter vectors, and difference
schemes, J. Comput. Phys. 43 357-372 (1981).

Rusanov, VV, Calculation of interaction of non-steady shock waves with obstacles,
J. Comput. Math. Phys. USSR 1, 267-279 (1961).

Shu, C. -W., TVB uniformly high-order schemes for conservation laws.
Mathematics of computation, 49,105-121(1987).

Shu, C. -W., Total-variation-diminishing time discretization. SIAM Journal on
Scientific and Statistical Computing, 9, 1073-1084 (1988).

Shu, C. -W, Essentially non-oscillatory and weighted essentially non—oscillatory
schemes for hyperbolic conservation laws, in Advanced Numerical Approximation
of Nonlinear Hyperbolic Equations, B. Cockbum, C. Johnson, C.-W. Shu and
E.Tadmor, Lecture Notes in Mathematics, volume 1697, Springer, 325-432 (1998).

Shu, C. -W., Different formulations of the discontinuous Galerkin method for the
viscous terms. In Advances in Scientific Computing, Z. -C. Shi, M. Mu, W. Xue and
J, Zou, editors, Science Press, 2001; pp. 144-145.

Shu, C. —W. and Osher, 8., Efficient implementation of essentially non-oscillatory
shock capturing schemes. Journal of Computation Physics, 77, 439-471(1988).

Steger, J.L. and Warming, R.F., Flux vector splitting of the inviscid gasdynamics
equations with application to finite difference methods, J. Comput. Phys. 40
263(1981).

van Leer, B., Towards the ultimate conservative difference scheme 11. Monotonicity
and conservation combined in a second-order scheme, J. Comput. Phys. 14,
361(1974).

van Leer, B., Towards the ultimate conservative difference scheme V. a second
order sequel to Godunov’s method, J. Comput. Phys. 32, 101-136(l979).

van Leer, B., Flux-vector splitting for the Euler equations, Lecture Notes in Physics,
170, 507 (1982).

Wang, Z.J., A fast flux-splitting for all speed flow, in Proceedings of Fifteen

International Conference on Numerical Methods in Fluid Dynamics, pp141-146,
1997.

121

65.

66.

67.

68.

69.

70.

71.

72.

73.

Wang, Z.J., Spectral (finite) volume method for conservation laws on unstructured
grids: basic formulation. J. Comput. Phys. 178, 210 (2002).

Wang, 2.]. and Liu, Yen, Spectral (finite) volume method for conservation laws on
unstructured grids H. Extension to two-dimensional scalar equation. J. Comput.
Phys. 179, 665-697 (2002).

Wang, 2.]. and Liu, Yen, Spectral (finite) volume method for conservation laws on
unstructured grids HI: one-dimensional systems and partition optimization. J.
Scientific Computing, 20,137-157 (2004).

Wang, Z.J., Zhang, L. and Liu, Y., Spectral (Finite) Volume Method for
Conservation Laws on Unstructured Grids IV: Extension to Two-Dimensional Euler
Equations, J. of Comput. Physics, 194, 716-741 (2004).

Wheeler, M. E, An elliptic collocation-finite element method with interior penalties,
SIAM J. Numer. Anal., 15,152-161 (1978).

Whitaker, D.L., Slack, DC, and Walters, R.W., Solution algorithms for the two-
dimensional Euler Equations on Unstructured meshes, AIAA 90-0697, 1990.

Williamson, C. H. K., Oblique and parallel modes of vortex shedding in the wake of
a cylinder at low Reynolds number, J. Fluid Mech. 206, 579 (1989).

Woodward, P. and Colella, P., The numerical simulation of two-dimensional fluid
flow with strong shocks, J. Comput. Phys, 54, 115-173 (1984).

Zhang, M. and Shu, C. -W., An analysis of three different formulations of the
discontinuous Galerkin method for diffusion equations. Mathematical Models and
Methods in Applied Sciences, 13, 395-413 (2003).

“Numerical simulation of compressible Navier-Stokes equations-extemal 2D flows
around a NACA0012 airfoil,” in GAMM Workshop, December 4-6, 1985, Nice,
France (Edt. INRIA, Centre de Rocquefort, de Renes et de Sophia-Antipolis, 1986).

122

VITA

Yuzhi Sun, she received a M. Sc. in Mathematics in June of 1989, and worked as

an assistant professor from 1989 to 1997 in the Mathematics Department at Hebei
University in China. She began working with her major professor Z.J. Wang at CFD lab
as a research assistant in Department of Mechanical Engineering at Michigan State
University from the fall, 2001. She obtained a M.Sc. in Mechanical Engineering at
Michigan State University in December of 2002. The following is her publications in this
period. Yuzhi will be a post-doctoral research associate with Professor Z.J. Wang in the
Department of Aerospace Engineering at Iowa State University.

Journal Papers:

Wang, Z.J. and Sun, Yuzhi, “Curvature-Based Wall Boundary Condition for the
Euler Equations on Unstructured Grids,” AIAA Journal, Vol. 41 , pp27-33, 2003.
Sun, Yuzhi and Wang, Z.J., “Evaluation of Discontinuous Galerkin and Spectral
Volume Methods for Scalar and System Conservation Laws on Unstructured
Grids,” International Journal for Numerical Methods in Fluids, Vol. 45, pp819-
838, 2004.

Sun, Yuzhi and Wang, Z.J., “Formulations and Analysis of the Spectra] Volume
Method for the Diffusion Equation,” Communications in Numerical Methods in
Engineering, Vol. 20, pp927-93 7, 2004.

Sun, Yuzhi and Wichman, Indrek S. “On Transient Heat Conduction in One-
Dimensional Composite Slab,” International Journal of Heat and Mass Transfer,
Vol. 47, pp1555-1559, 2004.

Conference Papers:

Sun, Yuzhi and Wang, Z.J. “High-Order Spectral Volume Method for the Navier-
Stokes Equations on Unstructured Grids,” AIAA-2004-2133.

Sun, Yuzhi and Wang, Z.J. “Evaluation of Discontinuous Galerkin and Spectral
Volume Methods for 2D Euler Equations on Unstructured Grids,” AIAA-2003-
3680.

Sun, Yuzhi and Wang, Z.J., “Evaluation of Discontinuous Galerkin and Spectral
Volume Methods for Conservation Laws on Unstructured Grids,” AIAA-2003-
0253.

Wang, Z.J. and Sun, Yuzhi, “A Curvature-Based Wall Boundary Condition for
the Euler Equations on Unstructured Grids,” AIAA-2002-0966.

123

   
 
 

ITY BRAR

lllijlilljjjlijijjjjlllwill

043