 

 

 

 

 

 

1.3.
(«'35.- ‘z Pg. SUE-r»
4‘ ON t: v
. NM. . .
it.
a .5 t

.r . 53

v ¢ A '1‘.

,v . .
a: .Luamzwhmwahw
. .

II y
‘ a 1.5
c ... .

km a
I] L, 3
... 1%
a.

.
‘-§ A
2...». 3
i2}...
ﬂan z. ..

é?

AA:
.65..

H..-

7.

.939
hate: 1.! :1
5.51.31? 5......
x 9.9.. 2 nmmnHNHZF:

 

 

meals

( ﬂog)

IHIIHNIHHIIIHIllll‘lllll(I‘lll’lllllllllllllllllll

31293 01771 0850

 

LIBRARY
Michigan State
University

 

 

 

This is to certify that the

dissertation entitled

Study on Optimized Upwind Schemes for
Computational Areoacoustics and Extension
of a Fully Conservative Chimera to Finite

Difference Schemes

presented by

Rangfu Chen

has been accepted towards fulﬁllment
of the requirements for

Ph.D. degreein Mech. Engr.

 

/"4 4’7 / C/‘-/

M ajor professor

Date November 19, 1998

 

MSU Ls an Afﬁrmative Action /Equal Opportunity Institution 0 12771

 

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

 

DATE DUE DATE DUE DATE DUE

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

STUDY ON OPTIMIZED UPWIND SCHEMES FOR ‘
COMPUTATIONAL AREOACOUSTICS AND EXTENSION
OF A FULLY CONSERVATIVE CHIMERA TO FINITE
DIFFERENCE SCHEMES

By

Rangfu Chen

A DISSERTATION

Submitted to
Michigan State University
in partial fulﬁllment of the requirements
for the degree of

DOCTOR OF PHILOSOPY
Department of Mechanical Engineering

1998

ABSTRACT

STUDY ON OPTIMIZED UPWIND SCHEMES FOR COMPUTATIONAL
AREOACOUSTICS AND EXTENSION OF A FULLY CONSERVATIVE
CHIMERA TO FINITE DIFFERENCE SCHEMES

By

Rangfu Chen

A computational methodology has been developed in the ﬁrst part of the thesis for the
simulations of acoustic radiation, propagation and reﬂection. The developed methodology
is high order accurate, uses less grid points per wave length comparing to standard high
order accurate numerical methods, and automatically damps out spurious short waves.
Furthermore, the methodology can be applied to acoustic problems in the presence of
objects with curved geometries. To achieve these results, high order accurate optimized
upwind schemes, which are applied to discretize spatial derivatives on interior grid points,
have been developed. High order accurate optimized one-side biased schemes, which are
only applied to discretize the spatial derivatives on grid points near computational bound-
aries, have also been constructed. The developed schemes are combined with a time differ-
ence scheme to fully discretize acoustic ﬁeld equations in multi-dimension in arbitrary
curvilinear coordinates. Numerical boundary conditions are investigated and intuitively
illustrated. Applications of the developed methodology to a sequence of one-dimensional
and multi-dimensional acoustic problems are performed. The numerical results have vali-

dated the developed methodology and demonstrated advantages of the methodology over

central-difference Dispersion-Relation-Preserving method. Numerical results have also
shown that the optimized upwind schemes minimize not only the disserpation error but
also the dissipation error, while retaining the numerical stability.

The second part of the thesis deals with a fully conservative Chimera methodology.
The fully conservative Chimera was originally developed based on a ﬁnite volume
approach. A ﬁnite difference scheme is shown to be identical to a ﬁnite volume scheme
with proper deﬁnition of control volumes and metrics. The fully conservative Chimera has
been successfully extended to ﬁnite difference schemes for viscous ﬂows including turbu-
lence models and successfully implemented into NASA’s widely used code OVERFLOW.
In the implementation Roe’s numerical ﬂuxes are employed to compute the inviscid ﬂuxes
across the patch interfaces. The thin layer approximation is used in the calculating of vis-
cous ﬂuxes. The implementation also includes implicit treatment of the patch interfaces,
which is shown to be extremely important to in stabilizing the overall time-marching pro-
cedure. The implementation of the fully conservative Chimera into the OVERFLOW
allows the error in the traditional Chimera approach due to non-conservative data interpo-
lations to be evaluated through direct comparisons. Tests on several inviscid, viscous lam-
inar and turbulent, exterior ﬂows with complex geometries are performed to validate the

implementations and demonstrate superiority of the conservative Chimera methodology.

ACKNOLOGEMEN T

I would like to express my sincere appreciations to the guidance committee, especially
my advisor Dr. Mei Zhuang, for their valuable discussions and encouragement. I would
like to thank Dr. Z. J. Wang of CFD Research Corporation (CFDRC) for his advice, help
and encouragement on the second part of the thesis. With his recommendation and help, I
have the opportunity to work in CFDRC where I gained a great deal of knowledge and
experiences in Computational Fluid Dynamics (CFD), computer languages, graphics and
operating systems. The work in the second part of the thesis was funded by NASA Ames
Research Center under Contract NASZ- 14226, with Dr. lPieter Buning being the technical
monitor, Dr. Z. J. Wang being the principal investigator. Helps from many colleagues from
CFDRC are very appreciated. The use of CFD-FASTRAN, CFD-GEOM, CFD-VIEW
from CFDRC are gratefully acknowledged. Prof. Christopher K. W. Tam of Florida State
University provided a sample code of the Dispersion-Relation-Preserving (DRP) method.
His contribution is also gratefully acknowledged.

Special thanks go to my wife Ms. Fengxia Ji for taking good care of me and our ener-
getic daughter Stephanie Chen, all the time. Without her great support, I would not have

been able to work on and ﬁnish the thesis.

iv

TABLE OF CONTENTS

LIST OF TABLES ............................................................................................... viii
LIST OF FIGURES ............................................................................................... ix
LIST OF SYMBOLS ............................................................................................. xv
INTRODUCTION ................................................................................................... 1
1.1 Background of Computational Aeroacoustics ............................................. l
1.2 Background of Conservative Chimera ......................................................... 6
1.3 Outline of the Thesis .................................................................................... 8
ACOUSTIC FIELD EQUATIONS ........................................................................ 11
2.1 Equations in Vector Form .......................................................................... 1 l
2.2 Equations in Cartesian Coordinates ........................................................... 15

2.2.1 Three-Dimensional Equations ....................................................... 15

2.2.2 Two-Dimensional Equations .......................................................... 17
2.3 Equations in General ‘Curvilinear Coordinates .......................................... 17

2.3.1 Three-Dimensional Equations ....................................................... 17

2.3.2 Two-Dimensional Equations .......................................................... 22
2.4 Wave Decomposition ................................................................................. 24
OPTIMIZED UPWIND SCHEMES ...................................................................... 27
3.1 Model Wave Equation ................................................................................ 27
3.2 Optimized Upwind Spatial Discretizations ................................................ 31

3.3 Extension of Optimized Upwind Schemes to Acoustic Field Equations
in Multi-Dimension .................................................................................... 48

3.4 Time Discretizations .................................................................................. 52

3.5 Boundary Conditions ................................................................................. 55
3.5.1 Optimized One-Side Biased Spatial Discretizations Near

Boundaries .................................................................................... 57

3.5.2 Wall Boundary Condition .............................................................. 60

3.5.3 Radiation Boundary Condition ...................................................... 63

3.5.4 Outﬂow Boundary Condition ......................................................... 66

APPLICATION OF OPTIMIZED UPWIND SCHEMES ..................................... 68

4.1 One-Dimensional Model Wave Equation .................................................. 68

4.2 Acoustic Radiation from an Oscillating Piston .......................................... 72

4.3 Acoustic Propagation in a Two-Dimensional Uniform Flow .................... 79

4.4 Acoustic Reﬂections .................................................................................. 83

4.4.1 Reﬂection of an Acoustic Pulse off a Straight Wall in a Two-

Dimensional Uniform Flow ........................................................... 83

4.4.2 Scattering of an Acoustic Pulse off a Cylinder .............................. 85

4.4.3 Reﬂection of an Acoustic Pulse off a Sphere ................................. 90

FIN ITE DIFFERENCE VERSUS FINITE VOLUME DISCRETIZATIONS ...... 93
5.1 Navier—Stokes Equations ............................................................................ 93
5.2 Finite Difference Versus Finite Volume Discretizations ............................ 97
5.3 Finite Difference Grid and Finite Volume Grid ...................................... 100
IMPLEMENTATION OF CONSERVATIVE CHIMERA IN OVERFLOW ....... 102
6.1 Review of the Chimera Approach ............................................................ 102
6.2 Numerical Fluxes through the Patch Interface ......................................... 104
6.3 Implicit Treatment of the Patch Interface ................................................ 109

vi

VALIDATION AND DEMONSTRATION OF THE FULLUY

CONSERVATIVE CHIMERA IN OVERFLOW ................................................. 115
7.1 Transonic Flow Through a Channel with a 10% Bump ........................... 115
7.2 Laminar and Turbulent Boundary Layers over a Flat Plane .................... 120
7.3 Grid Reﬁnement Study with Turbulent Flow over a Two-Element

Airfoil ...................................................................................................... 123
7.4 Turbulent Flow over a Three-Element Airfoil ......................................... 130
7.5 Summary .................................................................................................. 136
DEMONSTRATION OF THE FULLY CONSERVATIVE CHIMERA IN
OVERFLOW ........................................................................................................ 137
8.1 Hypersonic Flow Around a Launch Rocket ............................................. 137
8.2 Demonstration of Local Reﬁnement with Hypersonic Flow Around a

Launch Rocket ......................................................................................... 139
8.3 Flow over a Wing and Store Combination ............................................... 143
8.4 Transonic Flow Around a Combination of Wing and Missiles ............... 147
CONCLUSIONS .................................................................................................. 152
9.1 Summary .................................................................................................. 152
9.2 Suggestions to Future Work ..................................................................... 156
REFERENCE ....................................................................................................... 159

vii

LIST OF TABLES

Table 1 Coefﬁcients for the 7-stencil optimized upwind schemes ........................ 45

Table 2 Coefﬁcients for one-side biased schemes used in boundaries .................. 58

viii

Figure 1

Figure 2

Figure 3

Figure 4

Figure 5

Figure 6

Figure 7

Figure 8

Figure 9

Figure 10

Figure 1 1

LIST OF FIGURES

Real part (a) and imaginary part (b) of effective numerical wavenumber
versus wavenumber for 7-stencil schemes: Optimized upwind scheme
(N=4, M=2, P=4); Tam &Webb’s central DRP (N=3, M=3, P=4); Tam
&Webb’s backward DRP used near boundary (N=4, M=2, P=4);

Standard six-order accurate upwind scheme (N=4, M=2, P=6). ............... 47

Schematically showing the mapping of a physical space to a
computational space. .................................................................................. 57

Effective numerical wavenumber versus wavenumber for one-side biased
7-stencil schemes: Present optimized scheme used near boundary (N=5,
M=1, P=4); Tam &Webb’s backward DRP used near boundary (N=5,
M=1, P=4); Standard six-order accurate scheme (N=5, M=1, P=6). ......... 59

Schematically showing acoustic, vorticity and entropy wave sources in a
uniform ﬂow. .............................................................................................. 64

Comparison of exact and computed solutions at t=400 with different
schemes for the one-dimensional modal wave equation with an initial
Gaussian pulse. .......................................................................................... 70

Comparison of computational errors by different schemes for the one—
dimensional modal wave equation with a initial Gaussian pulse ............... 70

Comparison of exact solution and computed results at t=200 with (a): the
optimized upwind schemes, (b): Tam and Webb’s central DRP without
damping, for the one-dimensional modal wave equation with an initial
discontinuity. ............................................................................................. 7 l

Schematically showing the computational domain for acoustic radiation
from an oscillating piston ........................................................................... 72

Pressure distribution computed with the optimized upwind scheme at the
beginning of a cycle. .................................................................................. 75

Pressure contours at the beginning of a cycle computed with the
optimized upwind scheme .......................................................................... 76

Pressure contours at the beginning of a cycle computed with Tam and
Webb’s central DRP method without damping. ......................................... 76

ix

Figure 12

Figure 13

Figure 14

Figure 15

Figure 16

Figure 17

Figure 18

Figure 19

Figure 20

Figure 21

Figure 22

Figure 23

Figure 24

Figure 25

Figure 26

Comparison of pressure distributions along the axis of symmetry at the
beginning of a cycle. .................................................................................. 77

Comparison of computational errors for pressure along the axis of
symmetry at the beginning of a circle. ....................................................... 78

Schematically showing the computational domain, boundary conditions
and initial disturbances for acoustic propagation in a two-dimensional

uniform ﬂow. .............................................................................................. 80

Density distribution computed with the optimized upwind scheme
at t=80. ...................................................................................................... 80

The comparison of density distributions along the line y=x at t = 80 ........ 82

The comparison of computational errors for density along the line
y=x at t = 80. .............................................................................................. 82

Pressure distribution computed with the optimized upwind scheme _
at t=80. ...................................................................................................... 83

The comparison of computational errors for pressure along the wall
at t = 80. ..................................................................................................... 84

The comparison of computational errors for pressure along the line x=0
at t=80. ....................................................................................................... 85

Schematic diagram showing the computational domain and boundaries of
acoustic scattering problem ........................................................................ 86

Pressure distribution computed with the optimized upwind scheme at t = 7
for the ﬁrst case in which an acoustic pulse is released initially ................ 88

Time history of pressure ﬂuctuation at points A, B and C as indicated in
Figure 21. ................................................................................................... 88

Pressure wave pattern at t = 40 simulated by the optimized upwind
scheme ........................................................................................................ 89

Directivity of radiated sound computed by the optimized upwind
scheme. .................................................................................................... 89

Pressure patterns on the surface of the sphere and the plane of z = 0
at t = 2.5. .................................................................................................... 92

Figure 27

Figure 28

Figure 29

Figure 30

Figure 31

Figure 32

Figure 33

Figure 34

Figure 35

Figure 36

Figure 37

Figure 38

Figure 39

Figure 40

Figure 41

Figure 42

Figure 43

Pressure patterns on the surface of the sphere and the planes of x=0, y=0
and 2:0 at t=2.5 ......................................................................................... 92

Coordinate transformation and the grid cell around grid point (i, j, k) ...... 95

Schematically showing ﬁnite volume grid (solid line) and mutual ﬁnite
difference grid (dash line). ....................................................................... 101

Schematically showing cut cells, smallest faces and patch interface. ....103

Two-zonal overlapped grid system for the transonic ﬂow over a 10%
bump. ...................................................................................................... l 16

Pressure contours computed using OVERFLOW with the conservative
Chimera. ................................................................................................... l 17

Pressure contours computed using OVERFLOW with the original
Chimera. ................................................................................................... l 17

Comparison of pressure coefﬁcients proﬁles on lower wall of
the channel. ............................................................................................. 118

Comparison of convergence histories in using OVERFLOW with the
conservative Chimera and the original Chimera. ..................................... 118

Grids after local reﬁnement for the simulation of ﬂow over a bump by

OVERFLOW with the conservative Chimera. ......................................... 119
Pressure contours computed by OVERFLOW with the conservative

Chimera and with local grid reﬁnement ................................................... 119
Grids for the study of boundary layer growth over a ﬂat plane ............... 120

Velocity distributions cross interfaces between major and minor grids
inside the turbulent boundary layer over a ﬂat plane. .............................. 121

Velocity contours cross interfaces between major and minor grids inside

the laminar boundary layer over a ﬂat plane. ........................................... 122
Velocity contours cross interfaces between major and minor grids inside
the turbulent boundary layer over a ﬂat plane .......................................... 122
Velocity contours cross interfaces between major and minor grids inside
the laminar boundary layer over a ﬂat plane. ........................................... 122
The coarse grid for the turbulent ﬂow over the two-element airfoil ........ 124

xi

Figure 44
Figure 45

Figure 46

Figure 47

Figure 48

Figure 49

Figure 50

Figure 51

Figure 52

Figure 53

Figure 54

Figure 55

Figure 56

Figure 57

Figure 58

Figure 59

Patch interface for the conservative Chimera in the coarse grid .............. 124
Holes cut for the original Chimera in the coarse grid .............................. 125

Mach contours around two-element airfoil computed by OVERFLOW
with the conservative Chimera. ................................................................ 127

Mach contours around two-element airfoil computed by OVERFLOW
with the original Chimera. ....................................................................... 127

Pressure contours around two-element airfoil computed by OVERFLOW
with the conservative Chimera. ................................................................ 128

Pressure contours around two-element airfoil computed by OVERFLOW
with the original Chimera. ....................................................................... 128

Convergence histories for OVERFLOW with conservative Chimera and
original Chimera in the ﬁne grids. ........................................................... 129

Pressure distributions along the airfoil computed by OVERFLOW with the
original and the conservative Chimera in different levels of grids. ......... 129

Pressure distributions along the rear airfoil computed by OVERFLOW
with the original and the conservative Chimera in different levels of

grids. ........................................................................................................ 130
Medium grids and patch interfaces for the thee-element airfoil .............. 131

Pressure contours around the three-element airfoil computed using
OVERFLOW with the conservative Chimera on the ﬁnest grid. ............. 132

Pressure contours around the three-element airfoil computed using
OVERFLOW with the original Chimera on the ﬁnest grid. .................... 132

Pressure distributions along airfoil surfaces computed using OVERFLOW
with the original Chimera. ....................................................................... 134

Pressure distributions along airfoil surfaces computed using OVERFLOW
with the conservative Chimera. ................................................................ 134

Comparison of pressure distributions along the airfoil surfaces computed
using OVERFLOW with the conservative and the original Chimera. ..... 135

Comparison of pressure distributions for the inviscid ﬂow along the airfoil
surfaces computed with OVERFLOW. .................................................... 135

xii

Figure 60

Figure 61

Figure 62
Figure 63

Figure 64

Figure 65

Figure 66

Figure 67

Figure 68
Figure 69

Figure 70

Figure 71

Figure 72

Figure 73
Figure 74
Figure 75

Figure 76

Four zone Chimera grid for hypersonic ﬂow around a launch rocket ...... 138

Mach distributions on the surfaces of the rocket and the boosters and on
the last cross-streamwise plane computed using OVERFLOW with the

conservative Chimera. .............................................................................. 138
Mach contours on the surfaces of the rocket and the boosters. ............... 139
Grids around the rocket and the boosters with local reﬁnement .............. 140

Mach contours on the body surfaces and on the last cross stream plane
simulated with OVERFLOW with the conservative Chimera. ................ 141

Density contours on the body surfaces computed with OVERFLOW
with the conservative Chimera. ................................................................ 142

Mach distributions on the body surfaces and some sections simulated

with CFD-FASTRAN with local grid reﬁnement. ................................... 142
Mach contours on the body surfaces simulated with CFD-FASTRAN with

local grid reﬁnement. ............................................................................... 143
Grids and the conﬁguration of the wing-store combination. ................... 144
Patch interfaces, smallest faces and holes cut by the ZIG. ...................... 144

Pressure contours on the surfaces of the wing and the store and on one
spanwise section computed using OVERFLOW with the conservative
Chimera. ................................................................................................... 145

Mach contours on the surfaces of the wing and the store and on a
spanwise section computed using OVERFLOW with the conservative

Chimera. ................................................................................................... 146
Residual history for the simulation of ﬂow over wing and store

combination by OVEFLOW with conservative chimera. ........................ 146
Conﬁguration of the combination of wing and two missiles. .................. 147
Grids for the combination of wing and two missiles. .............................. 148
Holes automatically cut in the ZIG. ......................................................... 148

Pressure contours on the surfaces and two streamwise sections
simulated using OVEFLOW with the conservative Chimera. ................. 150

xiii

Figure 77

Figure 78

Figure 79

Pressure contours on the surfaces and two streamwise sections
simulated using CFD-FASTRAN with the conservative Chimera ........... 150

Pressure contours on the surfaces and two spanwise sections simulated
using OVEFLOW with the conservative Chimera. .................................. 151

Pressure contours on the surfaces and two spanwise sections simulated
using CFD-FASTRAN with the conservative Chimera. .......................... 151

xiv

A, B,C

A0 30 C0

ha

~

,B,C

>1

A0, in, 50

E if, 3
31%., 3.,
r

i

K.

matter;

(I:

LIST OF SYMBOLS

Roman Symbols
Jacobian matrices with respect to the conservative variable vector in x, y, 2
directions respectively

J acobian matrices with respect to the mean ﬂow primitive variable vector in
x, y, 2 directions respectively

Jacobian matrix for the smallest face

J acobian matrices with respect to conservative variable vector in the gener-
alized curvilinear coordinate §, 1], C directions respectively

Jacobian matrices with respect to the mean ﬂow primitive variable vector in
the generalized curvilinear coordinate 7;, n, C directions respectively

Inviscid ﬂux vectors in x, y, 2 directions respectively

Viscous ﬂux vectors in x, y, 2 directions respectively

Flux through the smallest face

Numerical approximation to the ﬂux through the smallest face
Determinant of the Jacobian matrix of generalized curvilinear transform
Stencils from the right and left

Conservative ﬂow variable vector

Speciﬁc gas constant, right cell, curvature radius of the wall

Vector of right hand side of discretized equations, magnitude vector in fun-
damental solutions

Cell face, face area

Right hand side of acoustic ﬁeld equations, source term vector

XV

Tg, TTI’ Tg

<¢<

<1

aNM
j’ j

i.j.k,l

Q9 ”U :v

‘IV

Temperature, eigenvector matrix in one-dimension
Eigenvector matrices corresponding to A, B, C or A0, B0, C0

Magnitude of velocity, cell volume

Velocity vector

Relative speed of sound to the mean ﬂow

Coefﬁcients in ﬁnite different schemes to discretize spatial derivatives
Coefﬁcients in ﬁnite different schemes to discretize temporal derivatives
Speed of sound, propagating speed of scalar wave equation

Speciﬁc heat at constant volume

Internal energy per unit mass

Total energy per unit mass

Grid point indices

The pure imaginary number

Face normal

Pressure

Primitive variable vector

Position vector

Time

Velocity components in x, y, 2 directions respectively

Wave phase velocity

Wave group velocity

Cartesian coordinates

xvi

9|

Bo

Y

At
Ax, Ay, Az

Ali. An, AC

Ag, V,
6..

U

5:, 51v 5r.

82’ ’34

AE’ An’ AC

95’ Am 9C

”L' “T
émé

Greek Symbols

Wavenumber

Effective numerical wavenumber
Parameter in optimization procedure
Speciﬁc heat ratio

Time step

Space steps in x, y, 2 directions respectively

Space steps in the generalized curvilinear coordinate 1;, 11,2; directions
respectively

Forward and backward-difference operators in 1; direction
Kronecker delta function

Central difference operators in the generalized curvilinear coordinate
E, 11, C directions respectively

Second order and fourth order dissipation coefﬁcients respectively
Angle

Generic mean ﬂow velocity component in the generalized curvilinear coor-
dinates

Thermal conductivity

Diagonal eigenvalue matrices corresponding to matrices A0, B0, C0 in the
generalized curvilinear coordinate g, n, I; directions respectively

Diagonal eigenvalue matrices corresponding to matrices A, B, C in the
generalized curvilinear coordinate §, 1], C directions respectively

Eigenvalue, weighting coefﬁcient in optimization procedure
Laminar viscosity and turbulence eddy viscosity respectively

General curvilinear coordinates

xvii

‘Ps

8|

()

()
()0
()L
()R
(),~,,~,k

( )'
( )'
( )‘"’

Density

Parameter in optimization procedure

Fundamental solution vector
Phase in the fundamental solution for scalar model wave equation

Angles

Phase in the fundamental solution for semi-discrete scalar model wave
equation

Fundamental solution for scalar model wave equation
Fundamental solution for semi-discrete scalar model wave equation
Frequency

Effective numerical frequency
Ancillary Symbols
Fourier transform, numerical approximation, quantity in the generalized
curvilinear coordinate
Numerical approximations to the inviscid and viscous ﬂuxes
Mean ﬂow quantity
Quantity evaluated at left cell L
Quantity evaluated at right cell R
Quantity evaluated at cell (i, j, k)
Matrix with only non-negative eigenvalues in matrix splitting procedure
Matrix with only non-positive eigenvalues in matrix splitting procedure

Disturbance with respect to the mean ﬂow quantity

Quantity at time level n

xviii

1 INTRODUCTION

The thesis contains two relatively independent parts. The ﬁrst part is about the study
and applications of optimized upwind schemes for Computational Aeroacoustics (CAA).
The second part is on extending a fully conservative Chimera algorithm to ﬁnite difference

schemes and implementing it in NASA’s widely used OVERFLOW code.

1.1 Background of Computational Aeroacoustics

Aeroacoustic problems arise in many engineering applications such as jet noise, auto-
mobile noise, fan and turbomachinery noise, propeller and helicopter noise, duct acous-
tics, and especially aircraft noise (Eversman, 1991). As air trafﬁc continues to increase
quickly, the Federal Aviation Regulation (FAR) is expected to enforce stricter rules which
require 34 dB reduction in current aircraft noise levels (Mecharn et al., 1994). The ability
to predict and reduce noise levels efﬁciently has become increasingly important to both
aircraft industry and government agencies such as NASA. As a ﬁeld of study on the gener-
ation and propagation of sound in aerodynamic ﬂow ﬁelds, aeroacoustics is getting more
and more attentions. Current methods used in aeroacoustics, however, rely on simple theo-
retical models which are not able to well characterize acoustic waves, or experimental
tests which are usually expensive and time-consuming to carry out (Motsinger et al.,
1991). CAA emerged as a relatively new area of research which focuses on numerically
predicting the generation and propagation of sound in aerodynamic ﬂow ﬁelds (Hardin et

al., 1995). With continuous advancement in computer hardware and software, CAA will

gradually reach its goal to accurately compute sound generating mechanisms, directivity,
spectral features and the propagation of the sound to the far ﬁeld. CAA has the potential to
ﬁnally provide design engineers a useful tool in a design environment to fully and rapidly
asses the effect that changes in the product or product environment have on the sound
ﬁeld.

I Aeroacoustic problems are governed by the same equations as those in aerodynamics,
namely the N avier Stokes equations. Aeroacoustic problems, however, have their own
nature, characteristics and objectives, which are distinctly different from those commonly
encountered in aerodynamics (Tam, 1995). First of all, aeroacoustic noises are inherently
unsteady which implies that accurate time-dependent calculations should be performed.
The need to accurately resolve the temporal ﬂuctuations in the calculation requires that a
sufﬁcient number of time steps are included in each period resolved. In order to perform
spectral analysis of the computed solution, a long-time computation is necessary. Second
of all, acoustic waves can cover a broad frequency range. To cover the range of human
hearing, for example, frequencies from 50 Hz to 20 kHz wound need to be resolved. This
would require both a large domain to resolve at least one wavelength of the longest wave
and a ﬁne mesh to accurately resolve the shortest wave. Third of all, ﬂuctuations caused by
acoustic waves are typically very small. The ﬂuctuation velocities associated with the jet
noise, for example, are four orders of magnitude smaller than the mean ﬂow, which could
even be smaller than the numerical errors incurred in the computation of mean ﬂow.
Finally, boundary conditions should be placed in the far ﬁeld such that acoustic waves can
go out of the boundary without reﬂections. These characteristics for aeroacoustic prob-

lems demand numerical schemes be almost free of numerical dispersions and dissipations

to enforce that the obtained numerical waves travel over a long distance with accurate
speed and magnitude, and make computation of aeroacoustic waves challenging.

Current computer hardwares are not capable of predicting sound generation and prop-
agation by directly performing three-dimensional, time dependent computation of ﬂow
past or through arbitrary complex geometries. The memory and CPU requirements by the
direct numerical simulation of the full Navier-Stokes equations are prohibitive, although
very high order accurate computational methods may be used so that the numerical error is
order smaller than acoustic ﬂuctuations. Alternative less expensive computational
approaches must be used. Non-linearity, however, is crucial to the noise generation pro-
cess which is related to turbulence in the ﬂow. Currently used turbulence modeling is a
barrier to resolve the generation process. Turbulence models, which are computationally
cheap, but model the physics necessary to adequately resolve the mechanisms of sound
generation, are needed.

Nevertheless, different approaches with some simpliﬁcations had been developed and
are currently being researched. The approach which combines near- and far-ﬁeld analysis
numerically solves the Navier-Stokes with turbulence models or large eddy simulation
(LES) near source region (Duque et al., 1995, 1996, Lim et al., 1993). The acoustic far
ﬁeld is determined with the Acoustic Analogy or Kirchhoff method. The Acoustic Anal-
ogy originally introduced by Lighthill (1952) is most developed method and widely in use
in the aircraft industry. The Kirchhoff method originally suggested by Hawkings is cur-
rently under development and gaining popularity as the higher resolution aerodynamics
and more powerful computer are getting available (Farassat, 1996). Another approach is to

ﬁrst solve the Navier-Stokes to obtain mean ﬂow solution. The set of linear equations gov-

erning the acoustic ﬁeld obtained by perturbing the N avier-Stokes equations about the
mean ﬂow are then solved. The effect of viscosity on the propagation of sound is neglected
A set of numerical schemes with least numerical dispersion and dissipation must be devel-
oped to solve the linear equations for acoustic waves. This approach is current trend of
Computational Fluid Dynamics based (CFD-based) CAA, which is under development
and is the least mature compared to the Acoustic Analogy and Kirchhoff methods

CFD has made rapid progress over the past two decades in many aspects including
grid generation and algorithms. Grid generation techniques, for example, had advanced
from very beginning Cartesian uniform grid to body-ﬁtted arbitrary curvilinear grids, from
single zonal grid to multi-zonal Chimera grids, from structured to unstructured polyhedral
grids (Mavriplis, 1995). High resolution schemes such as Total Variation Diminishing
(TVD) and essentially non-oscillatory (EN 0) schemes were developed in dealing with the
discontinuities such as shock waves in the ﬂow ﬁeld (Chakravarthy er al., 1987, Harten et
al., 1987, Yee, 1989). Various implicit schemes have been proposed and applied to
unstructured grids to accelerate convergence to the steady state (Blanco et al., 1998, Ven-
katakrishnan, 1995). CFD is being increasingly employed as a valuable tool for design in
many industries, such as aircraft industry for aircraft load prediction, and automotive for
underhood convection heat transfer, and etc. Current CFD commercial codes are able to
provide an integrate tool on the geometry deﬁnition, grid generation, physical model deﬁ-
nition, ﬂow solver and ﬁeld visualization (CFD Research Corporation, 1993, 1998). How-
ever, they are using standard low-order methods (usually second order) and are designed
for the computation of mean ﬂows (Buning et al 1991). Their applications to CAA are not

truly feasible because of the enormous computer memory and time requirements.

Many CFD schemes such as the MacCormarck scheme, standard upwind schemes and
ENO etc. have been extended to higher order using more stencil points and applied to the
computation of acoustic problems (Hardin et, al., 1995, Hixon, 1997). Many compact and
non-compact optimized schemes including Tam and Webb’s Dispersion-Relation-Preserv-
ing (DRP) schemes were designed for the linear acoustic waves (Kim et al., 1997, Lele,
1992, Lin et al., 1995, 1997, Tam and Webb, 1993, Tam et al., 1993, Zhuang et al., 1997,
1998, Zigg et al., 1993). It has been shown that for waves with high wavenumbers (short
waves) optimized schemes require less grid points per wavelength (PPW) than traditional
high order CFD schemes. The property of requiring less PPW is essential in CAA, since a
large computational domain is usually required. Most optimized schemes, however, are
restricted to central difference algorithms (Lockard et al., 1995, Zigg, 1997). This restric-
tion inevitably leads to stability problems. Spurious oscillations in numerical solutions are
generated, when there are discontinuities in ﬂow ﬁeld which are often encountered near
objects. The elimination of the spurious oscillations must be dealt with through the use of
ﬁlters or explicit dissipation terms. Although ﬁlters or dissipation terms are proved quite
successful in many acoustic problems, they are problem-dependent and require the prior
knowledge of the problems.

Upwind schemes have been widely used in CFD and shown very efﬁcient and robust
(Anderson et al., 1984, Barth, 1990, Warming et al., 1975, Li, 1997). Upwind schemes
ensure that waves propagate in correct physical directions. With the built-in dissipation,
upwind schemes automatically damp out high wavenumbers component in the solution.
Since the dissipation does not distinguish the spurious waves from actual acoustic waves,

upwind schemes need to be optimized so that the effects of dissipation on a large range of

acoustic waves are minimized. It is our objective of the thesis to investigate, develop and

apply high order, optimized, upwind schemes for CAA.

1.2 Background of Conservative Chimera

The generation of a computational grid may require up to eighty percent of total work
for a CFD ﬂow simulation with complex geometries. There is a great demand of automatic
grid generation to reduce the extensive requirement of the man power in grid generation.
Unstructured grid approach seems to be the most promising for the automatic grid genera-
tion and is getting dominant in commercial CFD ﬂow solvers nowadays. Unstructured
grid, however, imposes difﬁculties on the control of qualities of grid cells. Stability and
accuracy associated to numerical schemes are usually sensitive to the qualities of grid
cells. This is particularly true when dealing with viscous laminar and turbulent ﬂows.
Severely screwed cells certainly reduce the order of accuracy of numerical schemes and
bring out the instability to the numerical results. Furthermore, unstructured grid requires
more memory and CPU overhead because the connectivity of grid cells must be explicitly
described. In a structured grid the connectivity of grid cells is simply implied in grid indi-
ces. The structured grid, however, is very difﬁcult to handle complex geometries. Multi-
zonal overlapped structured grid approach (Chimera) was ﬁrst introduced about a decade
ago by Steger et al and Benek et al (Benek et al., 1983, Steger et al., 1983, 1985). Chimera
approach divides the computational domain into a number of geometrically simple sub-
domains (zones or blocks) for which independent structured grids can be easily generated.
Chimera approach keeps the simplicity of structured grid, while having the ﬂexibility in

grid generation and on the control of grid qualities. This approach can also dramatically

reduce the degree of difﬁculty and man power required for the grid generation. Since its
introduction to the CFD, Chimera approach has been very attractive and extended to solve
both steady and unsteady Navier—Stokes equations for complex geometries (Eberhardt er
al., 1985, Parks et al., 1991, Rogers et al., 1992, 1994, Soetrisno et al., 1994). Several
widely used CFD codes in the aerospace industry such as OVERFLOW (Buning et al.,
1991) and CFL3D have include the Chimera approach.

With Chimera, each independent grid is generated in each zone and the grids between
the neighboring zones are overlapped. During the computation the values of ﬂow variables
between the overlapped grids must be exchanged. In the traditional approach, the inter-
grid communications are accomplished through data interpolation. As a results, the con-
servation form of the numerical discretization is lost even if a conservative scheme is used
in each grid of the system. Although there are debates on whether conservation is neces-
sary for smooth ﬂows, it is generally recognized that the non-conservative data interpola-
tion approach may cause problems across discontinuities. There was an evidence that
conservation for the internal ﬂow through a nozzle must be guaranteed to correctly predict
the thrust of a rocket engine (Davis et al., 1992). Furthermore, the Navier-Stokes equa-
tions are obtained from the conservations laws and Lax (1972) mathematically proved that
the solutions of a conservative scheme guarantee the correct propagation speeds of discon-
tinuities provided the solutions convergence.

Since it was realized that the data interpolation has the deﬁciency of losing conserva-
tion when the Chimera was ﬁrst proposed, extensive research efforts have been devoted to
improved the conservative properties of Chimera (Berger, 1987, Moon et al., 1989, Wang

and Yang, 1995 ). But the full conservation in Chimera is not successfully accomplished

until recently Wang and Yang (1995) developed a fully conservative Chimera approach.
Wang et a! ﬁrst implemented the fully conservative Chimera into an inviscid ﬂow solver
and later successfully extent to tackle three-dimensional ﬂow problems for the ﬁrst time.
The fully conservative Chimera was developed originally based on ﬁnite volume
approach (Wang et al., 1994, 1995). The OVERFLOW which is being widely used in the
industry and government agency was written based on ﬁnite different approach. The
OVERFLOW includes the ability of Chimera with traditional data interpolation. It is our
objective to extend the fully conservative Chimera to the OVERFLOW. The implementa-
tion of the fully conservative Chimera into the OVERFLOW allows the errors in the tradi-
tional Chimera approach due to non-conservative data interpolations to be evaluated
through direct comparison. The OVERFLOW is a huge software packages for the simula-
tion of compressible ﬂows. The OVERFLOW of version 1.6ap after the implementation of
the fully conservative Chimera contains about a hundred thousand lines of source codes in
FORTRAN. The work of extension and implementation of the fully conservative Chimera
in the OVERFLOW is part of SBIR Phase H project of “A Fully Conservative Chimera
Approach for Structured/Unstructured Grids in Computational Fluid Dynamics” sup-
ported by NASA Ames Research Center with Dr. 2.]. Wang from CFDRC as the principle

investigator and Dr. Pieter Buning as the technical monitor.

1.3 Outline of the Thesis

Chapters 2, 3 and 4 are devoted to the ﬁrst part of the thesis, the study on optimized
upwind schemes for computational aeroacoustics and applications of the schemes to

acoustic problems. Chapter 2 derives the acoustic ﬁled equations in different forms from

the Euler equations by ignoring second and higher order small perturbation terms. The
properties and Jacobian matrices associated to the acoustic ﬁeld equations in the general-
ized curvilinear coordinate system are also derived in this chapter. Chapter 3 is devoted to
the development and construction of the optimized upwind schemes. The properties of
model wave equations are introduced. The schemes in one dimension are developed ﬁrst
by minimizing the difference of numerical effective wavenumber and actual wavenumber.
A fourth order accurate, 7-stencil, optimized, upwind scheme is given. The generation of
the schemes to the acoustic ﬁeld equations in the generalized curvilinear coordinate sys-
tem is included. Detailed boundary conditions and optimized one-side biased schemes
used for grid points near the boundary are also given. In Chapter 4, the developed 7—stencil
optimized upwind scheme is ﬁrst validated with one-dimensional model wave equation.
Then applications of the optimized upwind schemes to multi-dimensional benchmark
problems are carried out. Comparisons of the numerical solutions with analytical solutions
are made to show the advantages of the optimized upwind schemes over the central DRP
method. Finally the optimized upwind schemes are applied to the two-dimensional and
three-dimensional acoustic problems with the presence of curved objects.

Chapter 5, 6, 7 and 8 are for the extension of the fully conservative Chimera to ﬁnite
difference schemes, implementation in the OVERFLOW, and validation and demonstra-
tion of the implementation. Chapter 5 reviews the fully conservative Chimera based on the
ﬁnite volume approach. The differences between ﬁnite volume and ﬁnite difference dis-
cretizations are illustrated. Chapter 6 describes the implementation of the fully conserva-
tive Chimera in the OVERFLOW. In this chapter the emphasizes are on the implicit

treatments of the patch interface which is very sensitive to the stability and convergency.

Chapter 7 is devoted to the validation of the implementation of the fully conservative Chi-
mera in the OVERFLOW. A sequence of inviscid, laminar and turbulent ﬂows are simu-
lated by the OVERFLOW with the conservative Chimera. The direct comparisons
between the fully conservative Chimera and the traditional Chimera are made. Chapter 8
demonstrates cases to show the capability of the conservative Chimera in tackling more
complicate ﬂow problems.

Chapter 9 summarizes the conclusions obtained in the thesis. The suggestions for the

future work is also included in this chapter.

10

2 ACOUSTIC FIELD EQUATIONS

The generation and propagation of sound in aerodynamics ﬂow ﬁelds can be mathe-
matically described by Navier-Stokes equations. As mentioned in the introduction, it is
prohibitively expensive to predict sound propagation and radiation by the direct numerical
simulation of the full N avier-Stokes equations. On the other hand, the disturbances associ-
ate with the sound wave are small and the effect of viscosity on the propagation and radia-
tion of the sound is negligible. Therefore the modeling of acoustic ﬁeld can be based on
the linearization of equations governing the motion of an inviscid, non-heat-conducting

perfect gas, i.e. the Euler equations.

2.1 Equations in Vector Form

The continuity equation which describes the conservation of mass can be written as:
%?+(V-V)p+p(V-V) =0 (2.1)

_)
where p , p and V are density, pressure and velocity respectively. The inviscid momentum

equations, i.e. the Euler equations, can be written as the following:

8" -> -> _
p(§V+(V-V)V) - —Vp (2.2)

For the inviscid ﬂow without the absorption of radiation originating outside the system or

the local emission of radiation by the ﬂuid itself, the energy equation derived from the

11

conservation of energy is usually written as:

2 2

3%... Y2.)+ (ix . V)(e+ Y2.) = —V - (pi?) (2.3)

_)
where e is the internal energy per unit mass and V is the magnitude of velocity V, i.e.

V2 = V - V. Using the notation of substantial derivative (material derivative)
D a ->
— = —+ V . V .
Dr 3: ( ) (2 4)
and noticing
D V2 D —> Dv’
e
-D—t(e+-2—)—Ft+V'-D—t- (2.5)
V-(pV)=pV-V+V-Vp (2.6)

and combining the Euler equation (2.2), we can write the energy equation (2.3) in the fol-

lowing form:
De 9
For a perfect gas, the equation of state, the speciﬁc heat at constant volume cv and the
internal energy e can be respectively written as the followings:
p = pRT c = -R—-

" r—l " Y-lp

where T is the temperature, 7 is the speciﬁc heat ratio and R is the speciﬁc gas constant.

The speciﬁc heat ratio and the gas constant have different values for different gases. For
air at the standard condition, for example, y = 1.4 and R = 287 J / (kg - K). With the
above relations, we have

2‘?— _1_ 2 E __1__(QE_.122)
th " y-lth(p) " y-l Dt th (2'8)

Using the continuity equation (2.1) and substituting (2.8) into (2.7), we get the energy
equation for the inviscid perfect gas written as the following:

%IP+(V-V)p+yp(V-V) = O (2.9)

We use the continuity equation (2.1), the Euler equations (2.2) and the energy equation
(2.9) to derive the acoustic ﬁeld equations for the perfect gas. Consider the small perturba-

_)
tions on the mean state of a ﬂow, p0 , p0 and V0 , so that

P = Po+P'
p=po+p'
-) -> -)
V = V0+V'

The resulting acoustic ﬁeld equations, after second order and higher order terms in the
small perturbations in the continuity equation, Euler equations and the energy equation are
ignored, are written as the following:

8p'

57 +(i’10. V)p'+(V' - V)p0 + p0(V - V') + p'(V - V0) = 0 (2-10)

13

_)
p'(-51-Vpo) + 90(83—3 + (130- VW' + (3' - Vile) = —Vp' (2.11)
0

«13'

at +(V0 V)p'+(V'- V)p0+yp0(V V')+yp'-(V V0) = O (2.12)

They are representing acoustic continuity, acoustic momentum and acoustic energy
respectively. Here we have used the fact that the mean states pO , p0 and V0 satisfy the
continuity equation, Euler equations and the energy equation. From now on, the prime
used for the perturbations will be dropped for simplicity. Putting all the derivative terms of
the perturbations to the left hand side and all the terms without derivatives of the perturba-

tions to the right hand side, we can write the acoustic ﬁeld equations as the following:

3.2.02) V)p+po(V v>= -<‘v’-V>po—p<vi’o> (113)

gag/+(Vo- V)V+plO Vp = —(V V)V0+—129VP0 (2'14)
Po

8.5—? + (I30 - V)p +7pO(V - V) = - (I; ' V)Po‘YP(V ' i>10) (2'15)

They are called acoustic continuity, acoustic momentum and acoustic energy respectively.
It can be easily shown that the above equations will remain the same in non-dimensional
form providing the scaled factors are given as the follows: the density p is scaled by p, (a
reference density), the velocity V) is scaled by c, (the reference speed of sound), pressure
p is scaled by p re? , time t is scaled by L/c, where L is a suitable reference length, and

the spatial coordinates are scale by L. If there is no mean ﬂow, i.e. V0 = 0 , there are only

14

acoustic disturbances. From equations (2.13), (2.14) and (2.15), we have

aa—fwowil’) = 0 (2.16)
_)
av 1 _
57+an _ o (2.17)
9
g—€+Yp0(V- V) = 0 (2.18)

Now the acoustic continuity equation (2.16) and energy equation (2.18) are identical
because p = y?p . Taking time derivative on equation (2. 18) and substituting (2.17) into
0

(2.18), we have

32p pOV V _ 2

which is the standard second order scalar wave equation with the speed c0 of sound as the

propagating speed, where C: = 7:319 and A = V - V is the Laplace operator.
0

2.2 Equations in Cartesian Coordinates

2.2.1 Three-Dimensional Equations

. . . . 9
In three-dimensronal Cartesran coordinate system (x, y, 2), we have V = (u, v, w) ,
where u, v and w are the velocity components in x-, y- and z-directions respectively.

Denote variable vector 3 = (p, u, v, w, p) which is referred as primitive variable vector,

15

we can write acoustic ﬁeld equations in three-dimensional Cartesian coordinate system as:

a) 39 39 39_"
Eq+AO$q+BOEq+Coa—Zq - S (220)

9
where the components of vector S are the right hand side of acoustic ﬁeld equations
->
(2.13), (2.14) and (2.15). S depends on the ﬂow ﬁeld variations. For a uniform mean ﬂow
->
and no noise sources are presented, S = 0. The Jacobian matrices A0, Bo and Co, which

depend on the mean states only, can be written respectively as the following:

up -

“09000

0 OOuOOO ()

0 0 0u00
_0YP00 O'uo

 

 

 

 

va 0 po 0 0T
0 v0 0 0
Bo = o 0 v0 0 i (2.22)
Po
0 v0 0
_ Ypo 0 v0.1
r 1

(2.23)

 

 

16

2.2.2 Two-Dimensional Equations

_)
In two-dimensional Cartesian coordinate system, say, (x, y), then velocity V = (u, v) and

ﬂow variable vector 2)] = (p, u, v, p). The acoustic ﬁeld equations in two-dimensional

Cartesian coordinate system become:

.)
3372, + A0572); + 30312, = s (2.24)

The matrices A0, and Bo can be respectively written as follows:

 

 

-u0 po 0 0-
0 u o i

A0 = 0 p0 (2.25)
0 0 “0 0
O 'ypo 0 “0_
Fvo 0 po 0-
0 v0 0 0

B = (2.26)
0 0 0 v0 i
Po
_ 0 0 'YPo VOJ

 

 

2.3 Equations in General Curvilinear Coordinates

2.3.1 Three-Dimensional Equations

For a general curvilinear coordinate transformation in three-dimension, we have

g: §(x’y’z)

17

n = n(x,y,z)

C = C(x1y,z)

If the general curvilinear coordinate transformation is introduced into the equation (2.20),

and the Cartesian variable vector 2)] = (p, u, v, w, p) is retained as the dependent variable

vector, we have

%2} + Zigg—g?) + 80%?) + (303:5 = is (2.27)
where
I10 = Aogx + Bogy + C05, (2.28)
1730 = A0"). + Bony + ConZ (2.29)
60 = Aogx + B013), + Cogz (2.30)

and the Jacobian matrices, A0, 80 and C o , are

kx

O 60 0 O —

Po

~ ~ ~ k
A0130.0r C0 = 0 0 00 0 1 (231)

PO

k2

Po

0 Ypka Ypoky YpOkz eO

 

 

18

where

90 = uokx+kay+W0kz,

with k = E for A0, k = n for B0 and k = l; for C0. The Jacobian matrices A0, B0 and

C 0 have a set of eigenvalues and a complete distinct set of eigenvectors which are related

to the propagation of vorticity, entropy and acoustic waves. Therefore, similarity transfor-

mations can be used to diagonalize A0, B0 and C0, so that

where

A0 = TaAngl
80 = “Ann“:
60 = TCAchl

. ~ ~ ~ ~ 2 2 2 ~ 2 2 2
- D1ag(u0, “0’ “0’ uO + COJEX + 5y + £2, u0 — cO ﬁx + g), + éz) (2.32)
. ~ ~ ~ ~ 2 2 2 ~ 2 2 2
D1ag(v0,v0,v0,v0+c0,/nx+ny+nz,vo—COA/nx+ny+nz) (2.33)

. ~ ~ .. .. 2 2 2.. 2 2 2
D1ag(w0, W0: W0, W0 + cO,/§x + Cy + CZ, W0 — c0 Q, + Cy + CZ) (2.34)

 

= (“ng + voéy + woﬁz

= “Onx '1‘ Von), + WOnZ

{1'20 = “091: + VOCy + W09;

DiagOtl, A], ll, 71.2, 2.3) is a diagonal matrix with 7‘1: A], A], 12 and 13 as its ﬁve diago-

19

nal elements respectively. Co is the speed of sound which can be evaluated from the

expression c3 = ypO/po. Tk and Ak with k = C, n, C are referred as eigenvector matri-

ces and eigenvalue matrices, because they consists of the set of eigenvectors and eigenval-

ues respectively. The eigenvector matrices and their inverse matrices can be written as the

following:

"I
a:-
I
k-
N
O

 

-1_..

 

 

 

 

k 90 po
3 £60 «Eco
y 12 22

-12}: 2 _ﬂ

J5 Ji

0 12. i2.

J5 J2

0 Poco 8359
[2 J5-

~ ~ k
kz —y ——:
Co

.. 12y

O kx ——2
Co

~ k
_x 0 __§
co

J5 [2 fzpoc.
y kz 1

 

 

Where

20

 

(2.35)

(2.36)

.. k - k .. k
k, = x k y k = Z (2.37)

a y: ’ Z
,/ki+k§+kf ,/k§+k§+kf ,/k§+k§+kf

The eigenvalues of the Jacobian matrices determine the direction of propagation of waves.

 

 

 

For positive value of {to in the matrix A0 , for example, we can roughly say that the eigen-
vectors corresponding to 170 are propagating in the direction along which Ej-axis is
increasing. In another word, the eigenvectors are propagating in positive § direction.
Upwind schemes use different stencil points for wave propagating in different directions.
In the numerical methods shown in next chapter, the Jacobian matrices, A0, B0 and Co,
are split into two matrices. One has only non-positive eigenvalues and the other has only

non-negative eigenvalues. Therefore the following matrix is needed in the matrix splitting:

 

 

 

 

 

 

 

 

 

 

 

 

I} I? I?
A] p0 xcx p0 ya p0 2a __l_i
-2 .. . ~ ~ ~
k, 18,1.y 1.ka k,r
0 25+)” 2 B 2 B 2p0c0
1 I} I? k2 I? I? I?
T AT‘ = x y .2 y z y (2.38)
" " O 2 B 25”] 2 B 2p0c0
~ ~ ~ ~ ~2 ~
1.ka kykz k, k,
0 2 B TB 2‘3“”l 2p0c0
Pocolfx 905078: p060~z l
O 2 2 or 2 or 50.2+7t3)

 

 

where

A = Diag(}.l,7tl,)t,,kz,7t3) ,or = 12—13 and B = 924‘13‘27‘1'

21

2.3.2 Tho-Dimensional Equations

With a general curvilinear coordinate transformation in two-dimension

E = §(x,y)
n = n(x,y)
we have
a 9 ~ 3 -> ~ 3 -> _ 9
W+AW+BOT4 — S (239)
.310 = A0§x+80§y (2.40)
80 = AOnX+BOny (2.41)

and the Jacobian matrices, A0 and Bo , are

kx

~ ~ 0 00 O _
A0 or B0 = ‘30 (2.42)

k

0 0 90 1

Po

 

 

with k = t; for A0 and k = n for BO. Where 00 = “0kx+ voky. The Jacobian matrices
A0 and BO also have a set of eigenvalues and a complete distinct set of eigenvectors
which are related to the propagation of entropy, vorticity and acoustic waves. A0 and Bo

can be diagonalized and written as:

22

11>
o
l

—1
" T§A§T§
%=TAT

“1171

where

>
.m
I

. ~ ~ ~ 2 2 ~ 2 2
- D1ag(u0,u0,u0+c0,’§x+§ ,uO—co §x+§y) (2.43)

. .. - .. (2 2~ 2 2
An = D1ag(v0,v0,v0+cO nx+ny,v0—-c0 nx+ny) (2.44)

{to = u0§x+ V053);

VO = “011x + Von),

' 7

Po P0
ﬁco ﬁco
" 12x Ex
0 ky — ——
T7. = J5 «[5 (2.45)
0 :kx £2. -5

ﬂ J5
9000 13000

_ f2 72—.

p—a

 

 

-l

"I
R-
II
an an

(2.46)

 

 

 

23

with

 

 

(2.47)

By the same argument as in the three-dimension, we need to obtain the matrix TkAT;1

with the diagonal matrix A = DiagOLl, 7.1, 71.2, 2.3) . The manipulation of the matrix mul-

tiplication shows that

 

 

 

21.1 pkaa pOkya
2c0 200
122 1212
-1 0 .2113”, "2413
12,12 12
0 ,4 a) M
0 p0C0kx pocoky
_ 2 2

 

where or and B are deﬁned as in three-dimension.

2.4 Wave Decomposition

1

 

 

 

(2.48)

To illustrate that the Jacobian matrices have a set of complete distinct set of eigenvec-

tors which are related to the propagation of waves, we consider acoustic ﬁeld equation in a

uniform ﬂow in one-dimension, which can be written as:

24

-)
—— + A0— = (2.49)

where wave vector 3 = (p, u, p). The Jacobian matrix is

 

 

“0 Po 0

l

A = 0 u —

0 0 Po

0 rp u0

and it can be diagonalized as

1 Po Po F1 0 _l
Jic ,/2c c2
0 0 “0 O O 0

- 1
A0=TAT1= 0 i —i 01.0wO 0

l
«[2- 4/2 0 «73 £13000

0 0 u—c
Po Po 0 O 1 1

0 — — 0
L ﬁco ﬁco J2 «59060

 

 

 

 

 

It is seen that the matrix T contains three independent eigenvectors of the Jacobian matrix
A0. Denote these three eigenvectors as x1, x2 and x3 respectively. We can then project the
disturbance vector 3 on these three independent eigenvectors and write the disturbance
vector as 3 = 31X1 +32X2 + 83x3, where s,, 52 and s3 are the three projections in the
three directions. Substituting the above expression into equation (2.49) and using the prop-

erty of eigenvectors, we obtain

&(slxl) + angx-(slxl) + (2.50)
a a

5(52X2) + (“o + Co)§;(32X2) 4‘

$0310) + (u0 — ope—3150323) = 0

25

Therefore the disturbance vector can be decomposed into three vectors which are propa-
gating with speeds of “0 , uO + c0 and “0 - c0 respectively. Note that xl is traveling along
with the ﬂow and has no contributions to pressure disturbance p, which is referred as
entropy wave. x2 and x3 are propagating with the relative speeds of sound to the ﬂow and
they have inﬂuences on the pressure disturbances, which are known as acoustic waves.
Similar analogy can be applied to the eigenvalues and eigenvectors of the Jacobian matri-
ces in multi-dimension. The disturbance waves in the multi-dimension can be decomposed
into different waves, namely entropy, vorticity and acoustic waves. The entropy and vor-
ticity waves are propagating with the ﬂow and have no effects on the pressure disturbance,
while the acoustic waves are propagating with the relative speeds of sound to the ﬂow. The
contributions to the pressure disturbances are only from the last two columns of Tk in
equation (2.35) in three-dimension and equation (2.45) in two-dimension, which implies

that pressure disturbances are only affected by the acoustic waves.

26

3 OPTIMIZED UPWIND SCHEMES

In this chapter, we start from one-dimensional scalar model wave equation with a con-
stant propagating speed. The basic properties of the model wave equation, the Fourier
transform, the traditional ﬁnite different approach and the newly developed Dispersion-
Relation-Preserving (DRP) (Tam and Webb, 1993) on the model wave equation are brieﬂy
reviewed. The optimized upwind spatial discretizations are introduced and developed on
one-dimensional uniform grid ﬁrst. The generalization and extension of the spatial discret-
izations to the multi-dimensional uniform grids are derived subsequently. Finally time

integration and boundary conditions are discussed.

3.1 Model Wave Equation

Consider the scalar model wave equation

Bu Bu _

37 +c$ — 0 (3-1)
with the initial distribution

u(x,0) = u0(x) (3.2)

where c is the propagating speed of the wave and u = u(x, t). Without the loss of gener-
ality, we assume that c > 0. The solution for the above initial value problem of the equa-

tions (3.1) and (3.2) can be easily obtained and written as the following

27

u(x, t) = u0(x — ct)
or

u(x + ct, t) = u0(x)
which simply implies that the solution at point (x, t) is the initial value at x — ct. In another
word, the solution is a right traveling wave. For example, consider a “snapshot” of the
solution taken at to. After time of interval At has passed, the wave solution has moved to
the right (in the positive x-axis direction) by cAt.

Consider the Fourier transform of the solution u(x, t) and its inverse transform:

&(a, t) = ﬁ£u(x,t)e“‘°‘xdx (3.3)

u(x, t) = J“ ii(or,t)eiaxd0t (3.4)

where i = ./——1 and or is the wavenumber. Application of Fourier transform requires
some limitation on the solution u(x, t), e.g. u(x, t) is absolutely integrable. We assume
that Fourier transform is always applicable on the solution. The Fourier inverse transform
indicates that the solution u(x, t) can be decomposed into the inﬁnite sum of simple har-
monic waves. Applying the Fourier transform to the initial value problems of equations

(3.1) and (3.2), we have
d ~ . ~
Em" t) + rcoru(or, t) = 0 (3-5)

am, 0) = 170(1):) (3.6)

28

where
170(0t) = 517p.“ u0(x)e-iaxdx.

Solving equation (3.5), we obtain:

&(om) = &0(a)e"c°”

u(x,!) = r 120(a)eia(x_md0t (3.7)

Therefore the solution of initial value problem of equations (3.1) and (3.2) is the sum of
sinusoidal-traveling-waves with various wavenumbers. The expression (3.7) also indicates
that all the sinusoidal-traveling-waves are propagating to the right with the speed c. As a
matter of fact, for scalar linear wave equation, the behavior of the solution can be under-
stood and studied through the fundamental solution of the sinusoidal-traveling-waves in

the following form:

«ax—wt)

(p(x, t) = e (3.8)

where a) is the angular frequency. Denote ¢ = ax — on which is referred as the phase.

Substituting equation (3.8) into equation (3.1), we obtain
(0 = cor (3.9)

This is the dispersion relation for the one-dimensional scalar wave equations (3.1). It
states that for a wave with a given wavenumber or , its frequency a) is determined by the
dispersion relation. On the other hand, we may say that for a wave with a given frequency

(0, its wavenumber is determined by the dispersion relation. In general, for an arbitrary

29

scale wave equation, to satisfy the wave equation, the wavenumber or and the angular fre-

quency a) are related through an equation which can be formally written as the following:
6(0), or) = O (3.10)

The relation between 0) and or through the above equation is called the dispersion relation
of the wave equation. The wave phase velocity is deﬁned as vi, = m/or which is the prop-
agating velocity of the surface at which the phase does not change, i.e. 4) is a constant. If
angular frequency can be further expressed explicitly as a function of the wavenumber, i.e.
(1) = (0(a) , the wave group velocity can be introduced and deﬁned as v8 = dco/dor
which is the propagating velocity of different wavenumbers. For the one-dimensional sca-
lar wave equations (3.1), VI, = v8 = c. The wave phase velocity and wave group velocity
are the same constant regardless of wavenumbers. Although an initial arbitrary wave may
consist of the sinusoidal-traveling-waves with different wavenumbers, its waveshape and
its magnitude remain unchanged after arbitrarily large propagating distance. Therefore the
model equation (3. 1) is dispersionless and dissipationless.

For the system of wave equations in multi-dimension, equation (2.20), for example,
the behaviors of the solutions of the system can also be understood and studied through

a
the fundamental solution of the sinusoidal-traveling-waves. Denote X = (x, y, z) and

-> . .
or = (or), 02, a3) . The the fundamental solution can be wrrtten as

-> - ."_
choir) = 138‘“ X w” (3.11)

where R is a constant vector. By the same argument as that for the one-dimensional scalar

. -> .
wave equation, the or and o.) are related through the equation

30

-> ->
G((i),0t) = O (3.12)

The relation between (i) and (it through the above equation is called the dispersion relation
of the system. The phase velocity is now deﬁned as

3

.)
lar“

_)
VP:

which is in the same direction as a. This is because the surface with speciﬁc phase 410 is

given by

the normal vector of the surface 1S just or and the phase veloc1ty IS in the same direction as

the normal.

3.2 Optimized Upwind Spatial Discretizations

Finite difference schemes discretize partial differential equations by approximating
and replacing partial differential derivatives with ﬁnite differences. The discretizations
result in algebraic equations which can then be solved by various well-developed numeri-
cal methods. The algebraic equations resulted from the discretizations should be consis-
tent with the original differential equations and their solutions should be stable and
converge to the solutions of the differential equations. In another word, the solutions of the
algebraic equations approach the solution of the differential equations when the grid sizes
and time steps are approaching to zero. There are basically two ways to discretize the dif-

ferential equations. One is to treat spatial and temporal derivatives in the differential equa-

31

tions separately. In this approach, the discretizations of spatial derivatives result in a
system of semi-discrete ordinary partial differential equations in which the temporal deriv-
atives can then be further discretized by either well-developed methods such as Rung-
Kutta methods and multi-stages methods or by optimized time discretizations. The second
approach is to approximate the spatial and temporal derivatives simultaneously. No inter-
mediate semi-discrete forms are produced. Classical Lax-Wendroff and MacCormack
ﬁnite difference schemes belong to the second category. In this thesis, the optimized
upwind schemes following the ﬁrst approach only are studied and developed.

Refer to the scalar model wave equation (3.1). Consider the approximation of the ﬁrst

order spatial derivative Bu/ax by a ﬁnite difference on a uniform grid of spacing Ax

M:

Bu 1
(5)125 “rum (3.13)
j=-N
. Bu Bu . .
where “1+j = u((l+1)Ax,t), B—x I = 5;(1Ax, t) and aj is coefﬁcrent needed to be

determined according to required order of accuracy and other properties. I is an integer
representing grid node. For a pure initial value problem, 1 can be any integer. In real com-
putation, however, it requires ﬁnite computational domain. Therefore I must have lower
and upper bound. The ﬁnite difference approximation (3.13) requires M stencils of u to the
right and N stencils of u to the left of the node 1. Substituting (3. 13) into the model wave
equation (3.1), we obtain a system of semi-discrete ordinary differential equations which

can be written as the following:

(I M
u C
(3181—. .8 = o (3.14)

32

Consider the fundamental solution with a given wavenumber a in the form of

iale—wt
e( )

(pS(le, t) = (3.15)

where subscript s stands for semi-discrete to distinguish the fundamental solution of the
one-dimensional scalar wave equation (3.1). Substituting equation (3.15) into the semi-

discrete equation (3.14), we have
ei(a1Ax-W)[— i01+ i 2 ajeiaij] = 0 (3.16)

Denote

1'11ij

6: = aje (3-17)

M:

i
Ax
j-N

then equation (3.16) gives the dispersion relation of the semi-discrete equation (3.14) as

the following:

0) = ca (3.18)

It is observed that equation (3.18) is identical to the dispersion relation in equation (3.9), if
(it is replaced by or. By the analysis similar to the derivation of the solution for the one-
dimensional scalar wave equation (3.1) with the Fourier transform, it is necessary that at
approximates or accurately for relative broad wavenumbers to assure that the solution of
the semi-discrete equation (3.14) can be a good approximation to the solution of the wave

equation (3.1). (7 is called the effective numerical wavenumber of the ﬁnite difference

33

scheme (3.13). For ﬁxed spacing Ax, ('i is a periodic function of or. The spacing Ax
depends on the range of wavenumber which needs to be resolved in the real computation.
For example, if the scheme is being used to resolve the waves with the wave length longer

than 4Ax, in another word, using four grid points per wavelength (PPW), 27:5 2 4Ax or

E
2

it is seen that (TLAx is now a periodic function of orAx with 211: as the period.

orAx S . Therefore, it is convenient to consider orAx together as a variable. From (3.17),
Since {it can be a complex number, bi = 6L, + £611., where 6t, and 6t,- are the real and
imaginary parts of (it respectively. From the dispersion relation (3.18), we can write the

fundamental solution with a given wavenumber or for the semi-discrete equation (3.14) as

i(orle—c6rt) _ eialle—c(ﬁ,/a)tlecﬁ.-t

(950/315,!) = e (3.19)

The wave phase velocity VP: C(("ir/a) is not the constant c, unless a, = or which is not
realistic in a computation. vp depends on the wavenumber or. Waves with different wave-
numbers are propagating with different speeds. Therefore for an initial arbitrary wave
which may consist of the sinusoidal-traveling-waves with different wavenumbers, its
waveshape is no longer able to remain unchanged. This phenomena is called dispersion.
Note that if (3,- ¢ 0 , in addition to travelling with different speeds, waves with different
wavenumbers are now growing or diminishing with different factors. From the expression
(3.19), if there exists a wavenumber so that cﬁi > 0, the magnitude of the wave corre-
sponding to the wavenumber or increases without an upper bound as time increases. The
numerical solution is unstable and the discretization is referred as unstable. If cﬁi < 0 for
all wavenumbers, the magnitude of the solution has a upper bound as time increases. The

discretization is, therefore, referred as stable. The phenomena that the magnitude of a

34

wave grows and diminishes as it travels is called dissipation. Therefore, contrast to the dis-
persionless and dissipationless one-dimensional scalar wave equation (3.1), the semi-dis-
crete equation (3.14) is dispersive and dissipative if 6i, at or and 6t,- ¢= 0. Comparing

fundamental solutions (3.19) and (3.8), we have

(PAIN, t) ic(a - 6,): ca):
_ e

W — e (3.20)

Therefore the phase and magnitude errors introduced by the approximation of the one-
dimensional scalar wave equation (3.1) with the semi-discrete equation (3.14) are respec-

tively as;

|¢,(t) -¢(t)| = |c(a- 61ml

G(t) = eca‘t
where (1 represents the phase and G(t) is known as the ampliﬁcation factor. If the ampli-
ﬁcation factor has upper bound for all t, the numerical discretization is stable. For the
equation (3.14), it is easy to see the sufﬁcient and necessary condition for the stability is
ct-xl- S 0.

There are many ways to construct the coefﬁcients in the ﬁnite difference approxima-
tion (3.13) based on different criterions with which the wavenumber or is approximated by
the effective numerical wavenumber 6:. Traditionally the coefﬁcients are constructed
solely based on the order of accuracy of the truncation errors in the right hand side of
(3.13) or equivalently the right hand side of (3.17). Using Taylor expansion, we have

(ionij)2 + (iaij)3

1'01ij
6 21 3!

=1+iaij+

 

+ ...... (3.21)

35

Substituting the above expansion into the expression of the effective numerical wavenum-

ber (3.17), we obtain

(ian)3

. 2
(1%“,ij 2! (ajj3)+ ...... ] (3,22)

2!

 

, M
a = —.Aix 2 [aj+(iorAx)(ajj)+
j=—N

If we let

M
2 aj = 0 (3.23)
j=—N
M
2 ajj = 1 (3.24)
j=-N
M
Z ajjk = 0 f (3.25)
j=-N

for k = 2, 3, ..., N + M , there are N + M + 1 linear algebraic equations given in equa-
tions (3.23), (3.24) and (3.25) with N + M + 1 unknowns aj, j = —N, —N + 1, ..., M.
The unique solution of the equations will determine the coefﬁcients a J. of the ﬁnite differ-
ence approximation (3.13). Using these equations and equation (3.22) we can write effec-

tive numerical wavenumber as

ii = a+0(orAx)N+M (3.26)

From another point of view without employing the concept of wavenumber, we can

expand u I + j in the ﬁnite difference (3.13) at node I, i.e. the point [Ax , using Taylor series

36

as the following

 

_ Bu au (ij) 81150483
“i+j"“'+(ax,j)(Ax)+[a:2]1_212+[axsl 3! + (3'27)

Substituting the above expansion into the ﬁnite difference (3.13) and using algebraic equa-

tions given in equations (3.23), (3.24) and (3.25), we have

 

where the power of N+M is called the order of the truncation error. The coefﬁcient before

(Ax)N+M

in the truncation error depends on the (N+M +1 )th order of derivative of u with
respect to x. The N+M +1 )th order derivative of the fundamental solution with wavenum-

. N M 1
bet or contains the factor of or + +

. Therefore, these two different view points, using
and without using the wavenumber, are in fact equivalent. The one with the concept of
wavenumber, however, is preferred in this study because it infers how big the wavenumber
can be resolved, which is quite important in dealing with the acoustic waves.

We have discussed the ﬁnite difference scheme for which the coefﬁcients are deter-
mined purely from the standard Taylor expansions. The above procedures show that order
M+N of the truncation error is the maximum order which can be achieved with N+M+I

stencils used in the ﬁnite difference scheme. Therefore if orAx is small, which is true

either for small or (the relative long wave), or for small Ax (using a lot of grid points per

37

wavelength), the ﬁnite difference determined solely with standard Taylor expansions gives
the highest order of accuracy. One problem with this type of differencing, however, is that
the scheme quickly becomes inaccurate for relative short wave (large or) propagations or
relative course-grid discretizations (large Ax). As it is mentioned in the introduction,
acoustic waves have high frequency and broad wavenumber which demand the schemes
with less PPW. It is reasonable for the scheme to use 6 PPW, in another word the scheme
should be able to resolve the waves with orAx 5 it/ 3. When orAx is not close to 0, for
example, orAx is close to n/ 3 for the scheme with 6 PPW, the truncation error of the
scheme obtained solely with standard Taylor expansions is no longer small. To overcome
this problem, Miranker in 1971 presented a variety of difference schemes derived from
constrained minimization of L2 norms of the local truncation errors for a class of hyper-
bolic differential equations. Due to its complexity in solving the minimization analytically,
the technique apparently went relatively unnoticed until is was recently rediscovered by
many authors (Liu, 1996). This technique is used for deriving so called optimized schemes
including Tam and Webb’s Dispersion-Relation-Preserving (DRP) schemes in which L2
norms of truncation errors for both wavenumber in space and the frequency in time were
minimized.

Instead all the coefﬁcients in the ﬁnite difference approximation (3.13) are determined
by the Taylor expansion to achieve the maximum order of truncation error, some of them
are left to be free and used to minimize L2 norms of the truncation errors of the approxi-
mation of effective numerical wavenumber ('11 to the actual wavenumber or. To demon-
strate the idea, let P be an integer and 0 < P < M + N . Utilizing equations (3.23), (3.24)

and (3.25) for only k = 2, 3, ..., P, i.e. the ﬁrst P+1 equations instead of all M+N+I

38

equations, is unable to uniquely determine all the coefﬁcients in the ﬁnite difference
schemes (3.13), leaving M+N-P of coeﬂicients free. Without the loss of generality, we
assume the coefﬁcients from a_N, a_ N + l to a M_ P_, are free. From equations (3.23),
(3.24) and (3.25) for k = 2, 3, ..., P , we can express other coefﬁcients as the functions of

the free coefﬁcients, which can be formally written as the following

aj=aj(a_N,a_N+l, ...,aM_P_l),j = M-P,M-P+l,...,M

The above functions give inﬁnite sets of coefﬁcients in the ﬁnite difference (3.13) so that
ﬁnite difference (3.13) approximates the spatial derivative art/8x with P order of accu-

racy in the truncation error. In term of effective numerical wavenumber, we have

a = or + 0(an)P.
P is chosen to be greater or equal to one to guarantee that the ﬁnite difference approxima-
tion has at least one order of accuracy. The free coefﬁcients can be used to minimize the
following L2 norms of the local truncation errors introduced by the approximation of the

effective numerical wavenumber ﬁt to the actual wavenumber or:

E = “Mam-1) — [112 + (1 — x)[1m(ﬁ) + exp(—1n2(”3l(; “)2)]2}d0 (3.28)
45.

 

where B = orAx and E = &Ax. Re and Irn represent the real and imaginary parts of a
complex number respectively. [30 is predetermined number which gives the optimized
range of wavenumbers. The parameter 7. is a weighting coefﬁcient whose value is
between 0 and I. As it is discussed before, we consider orAx instead of or as an indepen-
dent variable. The ﬁrst term in the expression of E in equation (3.28) minimizes the dis-

tance between (TirAx and orAx in the sense of L2 norm. The second term instead

39

minimizes the distance between bit-Ax and 0, it minimizes the distance in the sense of L2
norm between aim and a Gaussian function, which is almost 0 when orAx is far from it
and is decreasing to -1 when orAx approaches to it. Remember that the wave propagating
speed c in the discussion is a positive number. Therefore the second term in the expression
of E, which forces the imaginary part of the effective numerical wavenumber to be nega-
tive during the minimization, guarantees the stability of the scheme as it is discussed
before. The Gaussian term is chosen to allow the imaginary part of the effective numerical
wavenumber to be very close to zero for the waves with wavenumber within certain scope
by adjusting 80. The second term also allows controls over short wave or high-frequency
damping by adjusting parameter 0'. From the expression of the effective numerical wave-
number in equation (3.17), using 8101ij = cos(jorAx) + isin(j0tAx) , we can write the
real part and imaginary part of B in the following forms:

M
Re(l3) = 2 a ,sinijB)

j=-N
_ M
mm) = — 2 a jcos(jB)
j=-N
Note sin(x) is an odd function and cos(x) is an even function. Im([—3) has same sign for
both positive and negative wavenumbers. Therefore we consider E in the following form

instead of (3.28):

a. M 2
E = ”L2 ajsin(jB)—B:| dB (3.29)
0 '='N

+ (1 - Min; ENajcos(jB)—exp(—ln 2(P%1t)2):l2dﬂ
0 = _

40

The necessary conditions that E is a minimum for all the free coefﬁcients are

QJ
D1

= 0,1 = —N,—N+ l, ...,M—P—l (3.30)

O.)

at

Directly taking the derivative of E in equation (3.29) with respect to the free coefﬁcients
a I , we obtain

M

(23—51: 22.19 :2»! a sin(jB)- 13:“: 2N8 a—Esinﬁﬂﬂdﬁ

”3a

+2(1—2.)j9 =2- Najcos(jB)—exp(—ln2(Bgn)2)][ 2 a—aicos(i[3):ld[3
0

i= —N

 

M M Ba
11;,” up ZN—jlsin(jB)sin(iB)- Bsin(il3)ldl3}

—N i=-

 

+2(l—A.) 2a12{ M a—‘:f[cos(j[3)cos(i[5)-- exp(— ln2(B; n)z)cos(iﬁ)]d8}

j=—N i=—Na
In the above expression for the derivative of E with respect to al, the integrations of
sin(jB)sin(iB), cos(jB)cos(iB) and Bsin(iB) over [0, [30] can be analytically calcu-

lated. Numerical integration methods, however, are needed to evaluate the integration of

 

_ 2
exp (— In 203 0' 1‘) )cos(iB) over [0, Bo] . In the real practice, the trapezoidal algorithm is
used to evaluate the integration. The partial derivatives of ai over a, appeared in the
above expression for aE/aa, can be obtained by solving the following M+N+I equations

forl = —N,—N+1,...,M—P—1:
—" = 5,, for i = —N,-N+ 1, ...,M-P—I (3.31)
and

41

. j_ _
2 j 52—, — 0 for — O, 1,...,P (3.32)

where 81-, is the Kronecker delta function (5,, = 1 if i=1 and 51-, = 0 if i¢ l) and equa-
tion (3.32) is from equations (3.23), (3.24) and (3.25) by taking derivatives with respect to
a,. Combining the equations (3.23), (3.24), (3.25) for k = 2, 3, ..., P with equation
(3.30), we can uniquely determine the coefﬁcients of ﬁnite difference approximation
(3.13) for a given order of accuracy P and parameters B0 , 7t and o.

If M=N and a_j = —aj , the ﬁnite difference scheme (3.13) becomes a central differ-
ence scheme. This is because it uses same number of stencils and weights in the left and
right of a node. From the expression of imaginary part of B derived above, for a central
difference scheme we have

_ M M

aim: Im(B) = —- 2 ajcos(jB) = a0— 2 (ajcos(jB) + a_jcos(—jB)) = a0

j=—N i=1

Equation (3.23) and a_j = -—aj imply that a0 = 0. Therefore the imaginary part of effec-

tive numerical wavenumber 6!. for a central difference scheme is zero. For the central dif-

ference scheme, the L2 norms of the local truncation errors E in equation (3.28) should not

include the imaginary part of [3. 7L and o are no longer necessary. Although the general
procedure to derive the coefﬁcients in the ﬁnite difference scheme (3.13) still applies for
the central difference scheme, slightly modiﬁcations are needed. For example, since
M M
2 0ij = zajjk[1-(—-l)k] = O for even k
j=—N i=1

the number of free coefﬁcients are now M-P/2 and P is an even number. Tam and Webb

42

(1993) studied and developed a fourth order accurate, 7—stencil, central ﬁnite difference
scheme in space with M = N = 3 in their DRP method for CAA. It was obtained
through minimization procedures described above for a central scheme with P = 4 and
B0 = 1|:/ 2. As we discussed before, the imaginary part of effective numerical wavenum-
ber 6t,- introduces the error in the wave magnitude. It seems in the ﬁrst appearance that
central difference schemes are more accurate than non-central schemes. But in reality,
because of the presence of objects, discontinuities in the ﬂow ﬁelds, computational errors
etc., the effective wavenumber corresponding to very short wavenumber may have small
positive imaginary part instead of zero. As it is pointed out before, this results in instabil—
ity. In fact, both the real and imaginary parts in the effective numerical wavenumbers con-
tribute to the numerical errors. The error in magnitude should be comparable to the error
in phase. Therefore, in equation (3.20), if (ii is in the same order of magnitude as a — 6t, ,
or in another word, error introduced by 6:,- is not dominant, non-central difference
schemes with negative (it. enforce the stability while keeping the competitive accuracy to
the central difference schemes.

Based on the above analysis it is possible for a non-central difference scheme to make
6t, — 0t and ("ii be almost zero for a wide range of wavenumbers, while keeping cﬁti nega-
tive to ensure the stability. Therefore this kind of scheme is expected to automatically
damp out spurious short waves while keeping the advantage of optimized scheme of using
less PPW. We studied 7-stencil of ﬁnite difference scheme (3.13) with M=2, N=4 and
P=4. This is an upwind biased scheme, because the wave is propagating to the right
(c > O) and the scheme uses more stencils upwind (from the left) than downwind (from

the right). If c < O i.e. the wave is propagating to the left, the upwind biased scheme uses

43

more stencils from the right than from the left. For example, the 7-stencil of ﬁnite differ-
ence scheme (3.13) has M=4, N=2 in the case of c < O. The coefﬁcients of the ﬁnite dif-
ference scheme (3.13) for c < 0 can be directly obtained from its counterpart for c > 0 , if
the same number of P as well as the same number of stencils are used. To see this, denote
a?” for the coefﬁcient a j in the ﬁnite difference scheme (3.13) with index j from -N to M
and (Ii-"N for those with index j from -M to N. The a?!" can be interpreted as the coefﬁ-
cients of the ﬁnite difference scheme which uses N stencil of points from the left and M
stencil of points from the right. a?!” has the similar meaning in its convention. Using the

same convention for B, note that

N M M
2 02ml). = 2 ail-NH)" = (—1)" 2 af'j’vj" (3.33)
j = -M j = "N j = -N

N M M
Re(B””) = 2 ajf’”sin(j[3) = 2 af’j”sin(—j(3) = — 2 af’j”sin(j[3) (3.34)

j=-M j—-N ‘=-N
_MN N N N

Im(B ) = - 2 a7~cos(jl3) = — 2 achos(—jB) = — 2 dill-”cos(jﬂ) (3.35)

j=-M j=-M j=-M
Ifwe chose

MN _ NM
“-1 " — 1
we obtain
MN .k k+l NM .1:
X J " (‘1) 2 1'
J=- j=-

and

Re<BM”) = Raﬁ”)

Im<BMN> = —Im(ﬁ”’”)

Therefore the scheme with coefﬁcient a?!” has the same order of accuracy as the scheme
. . NM . . .

wrth coefﬁcrent aj . The effectrve numerical wave number of the scheme wrth coefﬁ-
. MN . . . . NM . .

crent a 1- IS the conjugate of that of the scheme With coefﬁcrent aj . That opposrte srgn

in the imaginary parts of the effective wavenumbers of these two schemes are desired to

ensure the stability for the waves propagating in two contrary directions, namely positive

and negative x-axis directions.

Table l Coefﬁcients for the 7-stencil optimized upwind schemes

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

N=4,M=2,c>0 N=2,M=4,c<0
a3 0.016140071346698814 a3; 0.041382855555706463
afﬁ 0.12265083451 1 12346 a3? 0.44077420643183318
03 0.45448643568845881 aﬁ“ 0.50020513450976445
a3? 0.12475721579099250 af“ 0.12475721579099250
8232 0.50020513450976445 a3“ 0.45448643568845881
a? 0.44077420643183318 a§4 0.12265083451112346
a? -0.04138285555570646:l| ai“ 0.016140071346698814

 

 

45

 

Table 1 gives the coefﬁcients for 7-stencil optimized upwind ﬁnite difference scheme

obtained by the above described procedures with B0 = 1t/ 2, A = 0.964 and

0' = 0.26751t. Figure 1 compares the effective numerical wavenumbers for the optimized
upwind scheme (N=4, M =2, P=4) with coefﬁcients given in Table l to other well-known
schemes. It indicates that the optimized upwind scheme has much less dissipation error
than those of the standard six-order accurate upwind scheme and Tam &Webb’s backward
DRP scheme used for grid points near boundaries. It also indicates that the optimized

upwind scheme is able to resolve the waves with as high wavenumbers as those resolved

by Tam &Webb’s central DRP scheme. The parameter BO = 7t/ 2, )0 = 0.964 and
o = 0.267511: used for determining the coefﬁcients of the optimized upwind scheme were
obtained through trial and error. During the trial, a lots of combinations of [30, 3. and 0'

were used. For each combination the curves as those shown in Figure l were plotted and
compared. The coefﬁcients with the minimum errors based on visual evaluation were then
chosen. Therefore the coefﬁcients given in the Table 1 might not be the best choice among

the schemes with N=4, M=2 and P=4 because of trial and error.

46

0.05

 

 

 

 

 

 

 

 

   

  

 

 

 

 

 

 

 

 

 

 

 

0.04 ~ — Optimized Upwind (a) 4
- — - Tam 81 Webb’s Central DRP
0,03 _ —- - Tam & Webb's Backward DRP
---- Standard Upwind
0.02 L
21 0.01 f ’x’
E?
I 0.00 — _‘- '
k ”K-‘-__ o
lg —o.01 ~ "
-0.02 L
-0.03 —
-0.04 ~
-0.05 ‘ 1 * * 1 1 1 . 1
0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0
OLAx
0.05 T T Y I T I 1 r v I y
0 04 l ._ Optimized U ' (b)
. pwnnd —
— - - Tarn & Webb's Central DRP
0.03 L , —-- Tam & Webb's Backward DRP _
;; ---- Standard Upwind
0.02 F 3
0.01 ~
S... 0.00
Id
-0.01 F
-0.02 b
-0.03 L
-0.04 —
_0.05 1 1 1 1 1 1 1 \ 1 n J
0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0
an

Figure 1 Real part (a) and imaginary part (b) of the effective numerical wavenumber
versus actual wavenumber for 7-stencil schemes: Optimized upwind scheme (N=4,
M=2, P=4); Tam &Webb’s central DRP (N=3, M=3, P=4); Tam &Webb’s
backward DRP used near boundary (N=4, M=2, P=4); Standard six-order accurate

upwind scheme (N=4, M =2, P=6)

47

3.3 Extension of Optimized Upwind Schemes to Acoustic Field
Equations in Multi-Dimension

The acoustic ﬁeld equations in general curvilinear coordinate system in three-dimen-
sion have been derived in section 2.3. Here we rewrite the acoustic ﬁeld equations which

are given in equation (2.27):

39 99 39 ~39_"
374+ Aw” 0‘64 C0524”

The Jacobian matrices Ag, 30 and 50 can be diagonalized as

A0 = TgAngl,
~ —1
30 = TnAnTn ’
6‘0 = TgAngl.

The diagonal eigenvalue matrices are given in equations (2.32), (2.33) and (2.34). Tk and
T: with k = i, T] and C are given in equations (2.35) and (2.36). The diagonal matrices
determine the directions of wave propagating. The values of the diagonal elements in the
eigenvalue matrices may have different signs, which implies that the waves corresponding
to different elements may travel in opposite directions. In order to apply the optimized
upwind scheme which is derived from the model wave equation with a single propagating
speed, we need to split the diagonal matrices into two parts. One part contains elements
with only non-negative values and the other part contains elements with only non-positive
values. Therefore the waves corresponding to each part can then be treated as the waves
propagating along the same direction, either positive axis direction or negative axis direc-

tion. To be more speciﬁc, denote

48

+ -
Ag '—" A§+A§

ANA?)

f

+ .-
AC = AC + AC
where

A; = 500140 1,: lug—1A0

. 1 - 1

An = 5(An+|An|) An = 5(An’lAnl)
+ l - 1

AC = §(Ac+lAcl) A; = amt-Md)

Here we ad0pt the convention that the absolute value of a diagonal matrix is the diagonal
matrix which is resulted after each diagonal element in the matrices has been taken its
absolute value. Therefore diagonal matrices A; and A}, with k = E, n and C, contain

only non-negative and non-positive elements respectively. Let

.4; TéAngl 21;, = TgAng'

BS = T1422": Eb = TM:

&3 = TgAETEI 6;, = TCA'CTEI
The expressions for matrices .43, A6, 33, Eb, 63 and 06 can be obtained by substituting

the corresponding diagonal matrices deﬁned above to equation (2.38). From the properties

of the Jacobian matrices, we have

.20 = A's-+132)

49

60:654-60
The acoustic ﬁeld equation becomes
39 ~+39 “-39 ~ +93 ~ 39 “+39 "-39 9
— +A A +B5 +B = 3.36
a" 9:34“ ”a?" 071" WWWWW S ‘ )

We apply only 7-stencil optimized upwind schemes derived in the last section to dis-
cretize the spatial derivatives in the acoustic ﬁeld equations. As it can be seen later that the

procedures described here can in fact be used for any upwind schemes. Consider a grid of
uniform spacing A}; An and AC in 7;, 11 and C coordinates respectively. Denote node (i,

j, k) as the grid point with curvilinear coordinates (ii, 11]., Ck) which deﬁnes E,- = mg,

‘1,- = jAn and Q, = kAC. The cartesian coordinates of this node can be obtained from

the curvilinear coordinate transform. Denote 3,; 13k as the acoustic variable vector 3 eval-
uate at node (i, j, k) and at time t. Use the same convention of the subscripts i, j, k for
matrices and other variables. We can then approximate the spatial derivatives in the acous-
tic ﬁeld equation (3.36) by the optimized upwind ﬁnite difference schemes according to
the wave propagating directions as the following:

~ 3
(143%),“ = Al—E(A0)i,j,kl_Z-4a12 qi-I-l, j, k

A1€(A0)i,j,k 2 a124qi+l,j, k
I=-2

A
>1
Q.)
01°"
9
>4
1.
R'-
Dl

50

2
~+ a l ~+ 42
(30—21) = —(BO)i,j,k z 01 6911,14.“
3 r' 'k All
I l=-4

4
~- 3 9 _ 1 ~- 249
(BOEKILM — 33(30):,1', 111—2201 qi,j+1,k

2
~+39 l ~+ 429
(COW): j k = E(C0)i,j, let-2:401 qi',j,k+l

4
1 ~- 249
- [EU/‘0); j, k1 2201 (It, j,k+l

/—‘\
Q1
65’ '
J‘CIQ’)
1%»
”\_/
K.
R"
I

In the two-dimension, after the split techniques are applied to the Jacobian matrices,

the acoustic ﬁeld equation can be written as

39 ~+39 "-39 ~+39 ~-39_‘>
-a—tq+Ao-azq+Ao-a-Eq+B05-ﬁq+Bo§—q—S (3-37)

The diagonal matrices A; and A}, with k = § and T] and the split Jacobian matrices

.313, :16, 33 and 36 are deﬁned similarly as those in three-dimension. Expressions can be

obtained from the corresponding two—dimensional formulae given in section 2.3.2. For

example, the expressions for 3;, A5, 3; and 36 can be obtained by substituting the diag-
onal matrices to equation (2.48). The same discretization procedures as those in three-
dimension are then applied to approximate the spatial derivatives in the two-dimensional
acoustic ﬁeld equation (3.37). Note that the above described procedures of extension of
the optimized upwind schemes are limited to the acoustic ﬁeld equation. They can be

applied to any wave equations in general.

51

3.4 Time Discretizations

The acoustic ﬁeld equation (2.27) in three-dimension reduces to a system of ordinary

differential equations in time after spatial discretization, which can be written in the form:

431) __+..
(dt 1);.” R.,,,k (3.38)
where
9 ~+39 ~-39 ~+39 ~-39 ~+39 ~-39 9
Ri,‘,k = (A— +A +B— +B +C +C ) -Si, ,1: (3.39)
’ °a§q “3:4 “am" ”5174 “92" “WW ’

with the spatial derivatives at grid node (i, j, k) replaced by the corresponding optimized
upwind ﬁnite differences as described in the previous section. Similarly the acoustic ﬁeld

equation (2.39) in two-dimension becomes

d_?1 __->..
(ml ,- _ R1,) (3.40)
Where
9 ~+39 "-39 "+39 ~-39 9
Ri, ' = (A +11 +3 +8 ) -Si,' 3.41
, 0x4 0x4 0864 (06".,- , ‘ ’

with the spatial derivatives at grid node (i, j) also replaced by the corresponding optimized
upwind ﬁnite differences. In order to treat the systems of ordinary differential equations
(3.38) in three-dimension and (3.40) and in two-dimension together, we will drop the sub-

scripts which denote the spatial grid nodes.

52

Since acoustic wave propagation problems are generally not stiff, explicit methods are
appropriate. Optimized Adams-Bashforth and Rung-Kutta families for the wave propaga-
tion have been investigated by many researchers (Hardin et al., 1995). Many low-storage
multi-stage methods were introduced to integrate ordinary differential equations (Zigg et
al., 1993). The main purpose in the ﬁrst part of the thesis is to develop and demonstrate the
optimized upwind schemes. We therefore simply apply Tam and Webb’s optimized time
discretization developed in their DRP method to integrate equations (3.38) and (3.40)
(Tam and Webb, 1993). Let us refer to the model wave equation (3.1). In order to advance
to the next time level, Tam and Webb considered a four-level ﬁnite difference approxima-

tion to the derivative of 3u/3t in the following form:
(n - 1)
u(n + 1) _ u(n) z At 2 Egg) (342)

where At is the time step and superscript n is an integer which denotes time level. There-
fore variable with the superscript n means that the variable is evaluated at time nAt. For
example, um = u(x, nAt). Similar to the construction of a ﬁnite difference approxima-
tion to a spatial derivative, the coefﬁcients b I can be obtained solely through Taylor
expansion. As it is argued before, the ﬁnite difference scheme obtained with solely Taylor
expansion is good only when the frequency is relative low. Therefore the optimized proce-

dure is needed. To be more speciﬁc, after the time discretization and without the spatial

discretization the model wave equation (3.1) becomes a semi-discrete equation:

3 (n .. I)
(n + 1) 0!) 3a
— .4
u - u — ~Atl .2 0b,c(——ax) (3 3)

53

Consider also the fundamental solution with a given frequency 00 instead of wavenumber

or in the following form:
<pr, nAt) = e“°‘“ “A” (3.44)

Substituting the above fundamental solution into equation (3.43), we have

3
e'HoAt) - 1 = {—iAt 2 bled—101A!) )ca (3.45)
I = 0

Denote

. i(-—(DAI)
(T): ‘l" ‘1] (3.46)

3
At 2 bled—031A!)
l = 0

 

From equation (3.45) the dispersion relation of the semi-discrete equation (3.43) is given

as:
(T) = Ca (3.47)

Therefore the goal is to construct the coefﬁcients b I such that (T) is close to (o for relative
large band of frequency. (T) is then called the effective numerical angular frequency. The
procedures described for the optimized spatial ﬁnite difference schemes can also be
applied for the temporal ﬁnite difference scheme. Tam and Webb used Taylor expansion to
ensure the (T) approximates a) to order of (At)3 , i.e. (T) = o) + 0(At3) and minimize the

error L2 norm E:

54

0.5
E: j {MRe((T)At) — (0At]2 + (1 — k)[Im((T)At)]2}d(0)At) (3.48)
—0.5

They chose 3. = 0.36 and obtained the coefﬁcients bO = 2.3025581 , bl = —2.4910076,

122 = 1.5743409 and b3 = —O.38589142.

Using Tarn and Webb’s temporal ﬁnite difference scheme, we can discretize the sys-
tem of semi-discrete ordinary equations (3.38) in three-dimension and (3.40) in two-
dimension. After dropping the subscripts for the grid node and discretizing the temporal
derivatives in equations (3.40) and (3.38), we obtain the full—discrete algebraic equations

in the following form:
1) < ) 3 90: I)
+ _
5‘" = 2’)" -At 2 b,R (3.49)
l = 0

I

with 23‘"- ) given by equation (3.39) in three-dimension and by equation (3.41) in two-

. . . (n -l) .
dimensron. In our computanons, we let 1% = O for n S 1, 1n another word we

assume that time derivatives are zero when the actual time is less than zero.

3.5 Boundary Conditions

Since the computation domain can not go to inﬁnity, outer boundaries which enclose a
ﬁnite interior domain are required. Through the outer boundary, acoustic waves are radiat-
ing from the interior domain to the outer region. There are also ﬂuids leaving and entering

the interior domain through the outer boundaries. Since the outer boundaries are chosen

55

far enough away from the source, ﬂows near outer boundaries can be treated to be uni-
form. Therefore there are no disturbances coming into the interior domain. As it is dis-
cussed before, the disturbances consist of entropy, vorticity and acoustic waves, each
having distinct wave propagation characteristics. The outer boundary conditions should be
transparent to allow the disturbances to leave the outer boundaries without signiﬁcant
reﬂections. Since acoustic waves are propagating with the relative speed of sound to the
ﬂuid velocity, and the entropy and vorticity waves are propagating along with the ﬂow,
there are basically two types of outer boundary conditions, namely radiation and outﬂow
boundary conditions. Radiation boundary condition is used only for acoustic waves radiat-
ing to the far region through the boundaries. The outﬂow boundary condition, however, is
used for the combination of outgoing acoustic, entropy and vorticity waves. Near the sur-
face of solid objects inside the interior computational domain wall boundary condition is
needed to enforce the non-permeating ﬂow condition. Signiﬁcant progress has been made
in the development of numerical boundary conditions for high order ﬁnite difference
schemes on uniform Cartesian grids (Hixon et al., 1995, Tam 1997, Tam et al., 1996).
Radiation and outﬂow boundary conditions proposed by Tarn & Webb based on asymp-
totic solutions of the governing equations have been numerically shown to be very effec-
tive. They derived the radiation and outﬂow boundary conditions by applying Fourier-
Laplace transform to the two-dimensional acoustic ﬁeld equations in a uniform ﬂow along
the positive x—axis and employing sophisticated mathematical techniques of integration
analysis. In this section, we derive the radiation boundary conditions through intuitive
observations in both two-dimension and three-dimension. We will also discuss 7-stencil

forth order one-side biased ﬁnite differences (M+N: 7, P=4) used to discretize the spatial

56

derivatives near the wall and outer boundaries.

3.5.1 Optimized One-Side Biased Spatial Discretizations Near Boundaries

The general curvilinear coordinate transformations given in Chapter 2 are employed to
map the physical domain into computational domain in which the body surface is usually
transformed to a boundary of computational domain permitting easy application of surface
boundary conditions. Figure 2 schematically shows mapping a two-dimensional physical
domain with a wall into a computational plane. The index of grid points in n direction
shown in Figure 2 starts from the minimum 1 which denotes the wall to the maximum JM
which represents the outer boundary. Since only the values of variables inside the compu-
tational domain are available, one-side biased spatial ﬁnite difference schemes are needed

to discretize spatial derivatives for the grid points near the boundaries.

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

JM 5 i §(x, y)
11 " “(x’ y) M
2
y TI 1
K L_’x / i g Ghost Point
Physical Plane Computational Plane

Figure 2 Schematically showing the mapping of a physical space to a computational
space

57

Table 2 gives the coefﬁcients a? and 060 of 7-stencil one-side biased spatial ﬁnite

difference schemes. Here the subscriptions have the same meaning as those described

before. a?! was obtained by the procedure described in section 3.2 with P = 4

[30 = 1t/ 2, k = 0.01 and o = 0271:. 060’ however, represents standard six order accu-

rate one-side scheme and was obtained solely with Taylor expansion. The reason for

choosing the standard scheme for 060 is that we found based on visual evaluation that

C O I O O O 15 5]
optimization did not make any improvement for this case. Note that a j = -a_j and
06 _ _ 60
a] - a_j 0

Table 2 Coefﬁcients for one-side biased schemes used in boundaries

 

 

 

 

 

 

 

 

 

N=5,M=1 I: N=6,M=O

51 60

a_5 0.0283927780285557 a_6 1/6
51 n 60

a_4 0.192107686530459 a_5 -6/5
51 l 60

a__3 0.617980095557289 a_4 15/4

a5; 1.28536574416235 afg -20/3

a3} -2.14340185401956 a3 15/2

agl 1.11244509330028 af‘,’ -6
51 60

a, 0.199856203612323 a0 49/20

 

 

 

 

 

 

 

58

0.05 I l ' l T T

 

 

' (a)
0.04 ~ Present 1' l 4

— - - Tam 8. Web 1' I
0.03 — - -- - Standard ' , a

 

 

 

 

0.02 ~
a
51 0.01 i
5
1 0.00

 

 

L
'3 -0.01 4

 

-o.02 » .
-0.03 ~ 4
—0.04 ~ —

 

 

_ .05 1 1 1 1 L 1 1 i 1 1
0 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0
OtAx

 

0.05 T I ' I I I '1

 

 

I
0.04 L Present '
— - — Tam & Web ’
0.03 r ---- Standard ’

 

 

 

 

0.02 ~ ,’ / ‘\
0.01 —

a
<.L 0.00
15

 

-0.01 L
-0.02 r
-0.03 F
-0.04 ~

 

 

 

_0.05 1 1 . 1 1 - .
0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0

01A):

 

Figure 3 Effective numerical wavenumber versus wavenumber for one-side biased
7-stencil schemes: Present optimized scheme used near boundary (N=5, M=l,
P=4); Tam &Webb’s backward DRP used near boundary (N=5, M: l, P=4);
Standard six-order accurate scheme (N=5, M = l, P=6)

59

Figure 3 compares the effective numerical wavenumbers for the present optimized
scheme (N=5, M = l, P=4) whose coefﬁcients are given in Table 2 with other schemeslt is
shown that the standard six order accurate scheme gives the best approximation of the real
part of the effective numerical wavenumber to the actual wavenumber for GA): 5 1 among
all schemes shown in the ﬁgure. Its imaginary part of the effective numerical wavenumber,
however, gives the poorest approximation to the zero. As it is discussion in section 3.2, the
numerical errors come from both the real and imaginary parts of the effective numerical
wavenumber. Good scheme should balance the errors from both parts. Figure 3 shows that

the optimized scheme is the best choice based on overall accuracy for OtAx S 1 .

3.5.2 Wall Boundary Condition

To illustrate the general ideal of imposing the wall boundary condition we consider
here only in two-dimension, because it can be easily extended to three-dimension. We also
assume that the curvilinear coordinates on the wall are orthogonal which is important to
ensure the accurate implementation of the wall boundary condition. For non-orthogonal
coordinates other approaches for implementing high order wall boundary condition should
be implored. It requires extensive study for this topic alone as it is pointed out by Tam
(1997). On a solid wall the normal velocities of both the mean ﬂow and the acoustic wave
are zero. The Jacobian matrix in the normal direction has only two non-zero eigenvalues.
One is positive and the other is negative. It implies that only one boundary condition
should be imposed on the wall based on local characteristic analysis. From the acoustic
ﬁeld equations, we can express the derivative of pressure disturbance in term of velocity,

density and curvature only. The derivative of pressure disturbance is used to provide the

60

boundary condition on the wall.

Let’s refer to Figure 2 in which the wall is represented by a constant 11 line in physical
plane. We derive the derivative 3p/3n on the wall ﬁrst. Denote 5t and y the unit vectors

in Cartesian x-direction and y-direction respectively. I: and f] are the unit vectors in the E

and 1] directions respectively. We have

W)
I

- (xi)? + ygyy x: + y:

x ,. ,. 2 2
'1 — (xnx + yny)/,/xn + y,I

Note that 1] axis on the wall is now in the same direction as the normal 3. On the wall we

_) A A
have V = V15. Using the orthogonality of E, and ﬁ , and {ther metric properties of the
transform, on the wall we have I

“0x11 + voyn = 0, uxn + vyTl = 0
“0nx+v0ny = O, “lei'Vle = 0
Multiply x“ to the second component equation and y11 to the third component equation in

acoustic ﬁeld equation (2.39). Adding them together and noticing
11.x“ + 11 yyn = 1
we obtain
1 3p
(1“ng + voﬁyxxnug + ynvg) + 93-0- = Jrns2 + yns3

where s2 and 33 are the second and the third components of the known vector 3‘ in the
right hand side of equation (2.39). Use the fact (xnu + ynv)§ = 0 and

x1)“: + ynvg = (xnu + ynv)g — “xﬁn — vygn = —(ux§n + Vyén)

61

Finally on the wall we have the expression for the pressure derivative,

3% = p[(u0§x + voﬁyXuxgn + vygn) -1-(xns2 + yns3)] (3.50)

This expression can also be derived by linearizing the well known expression of

3p: sz / R and employing other facts about the transform. Here R is the curvature
radius of the wall.

We are now implementing the wall condition. Again refer to Figure 2 which also
shows that a level of ghost points denoted by the index g is added to the computational
domain. The level of ghost points is used to impose the pressures derivative on the wall.

One the wall we let

3p
(—)1 = Tn[2::a15"pj+a“pg]

with 3% given in equation (3. 50). We solve the above equation for the only unknown p g
and store its value in the ghost point g. For all variables except for pressure one-side

biased ﬁnite difference schemes are used with index starting from I. For example, for u at

index I, 2, 3 and 4 we respectively use

7
3a 1 06
('37,), = 302 01-1“)
1=1
7
314 1 15
(a), = raga--216-
7
3u 1 24
(3)3 = 33211010“)

62

7
au 1 33
(55). = A_n 7- “1—4‘0‘

J=1

where the coefﬁcients a? are from Tam and Webb’s central DRP scheme which is used

for one level of grid points near boundaries only. For pressure, the one-side biased ﬁnite
difference schemes are similarly used. The index, however, includes the ghost point g
instead of starting from index I. For example for index 2, we use
3p 1 6 24 24
(5)2 = -A—n[j§=;1 ai" 2p} + a-Zpg]
We should note here it is not necessary to calculate dp/dn at ghost point 3 which is not
used in discretization of acoustic equation, and the derivative at index 1 which is already

known from its expression.

3.5.3 Radiation Boundary Condition

As it is pointed out that the outer boundary is chosen far away from the source, there-
fore the uniform ﬂow near outer boundary can be assumed. The behavior of the acoustic
disturbance approaching outer boundary from asymptotic analysis is that it is propagating
at the relative speed of sound to the ﬂow. Refer to Figure 4 which schematically shows
vorticity, entropy and acoustic waves near the outer boundary in a two-dimension uniform

ﬂow. The acoustic disturbance is approaching to the point (x, y) on the outer boundary

with the speed 9 as shown in Figure 4 in two-dimension. We can describe the behaviors

of the acoustic vector 2’] at a point on an outer boundary with distance r from the source as

63

-> ~ -> ~
2’1= f(r-Vt) and 2’1: f(r—Vt)

J; r

for two—dimension and three-dimension respectively (W hitham, 1974). Taking the deriva-

tive of 3 with respect to time t and radius r; we have the following equations

3* “ 3 9 l9 _ 3.51
Eq+v(§?q+2rq)'0 ( )
39 ~ 39 19 _

—atq+ V(—qar + -;q) — 0 (3.52)

for two-dimension and three-dimension respectively.

 

 

acoustic wave (x, Y)
//
13
O
VOI'thlty,
_ entropy
4>_B and acoustic
j . waves
acoust1c wave

 

 

 

Figure 4 Schematically showing acoustic, vorticity and entropy wave sources in a
uniform ﬂow

Equations (3.51) and (3.52) are transferred to the curvilinear space and discretized with
one-side biased ﬁnite difference schemes. Near the outer boundary where there are only
acoustic waves radiating to far ﬁeld, we solve the above equations instead of the acoustic
ﬁeld equations. But the speed {7, however, must be determined ﬁrst. From the triangle

shown in Figure 4, it can be easily checked that
l7 = cocos<l>1 + Vocos<l>2 (3.53)

Let i" be the unit vector pointing from the source center to the point (x, y) on the outer

boundary shown in Figure 4. Notice that

 

 

 

 

Therefore we have

 

.. 2
v = cﬁ—V§+(i’lo-r) +6304) (3.54)

In two-dimension we can write {I in term of components of mean ﬂow velocity and
the angle 0 between the vector i" and x-axis shown in Figure 4. Chose the center of source

as the origin of the Cartesian coordinates, then we have
? = (cosB, sine)

Noticing that

_)
V0 = (uo, v0)

65

2 2 2
V0=u0+VO
and
2 -> .2 2 2 . 2 . 2
VO—(Vo-r) = u0+v0-(u0cose+v0sm9) = (uOSInG—vocose)
weobtain

 

{1(9) = JCS - (uosinO - vocose)2 + uocose + vosinO (3.55)

which is the same as those obtained by Tam and Webb with Fourier-Laplace transform and

integration techniques.

3.5.4 Outﬂow Boundary Condition

As it is discussed in chapter 2, the outgoing vorticity and entropy waves are propagat-
ing with the ﬂow and do not have effect on the pressure disturbance. Therefore the equa-
tion (3.51) in two-dimension and equation (3.52) in three—dimension should also be
applied as outﬂow boundary condition for the pressure disturbance. For other variables,
the acoustic ﬁeld equations are used. Since uniform ﬂow can be assumed near outer
boundaries, all the derivatives of mean ﬂow variables are zero. The acoustic ﬁeld equa-

tions become

%§+(v’o-V)p+p0(V-v’) = 0 (3.56)
33 —> —> 1
—+(Vo-V)V+—Vp = O (3.57)
at Po

66

g—l:+(i30-V)p+yp0(V-t7) = 0 (3.58)

Combined with acoustic energy equation, the acoustic continuity equation can be written

as

‘12 " . _ __1_ 5113 " .
at+(Vo V)p _ C(z)[at+(V0 V)p] (3.59)

In the computation near the outﬂow boundaries we use the radiation boundary condition

for the pressure p, acoustic continuity equation (3.59) for the density p , and acoustic

-> . . . .
momentum equation (3.57) for the velocity V. These equatlons are also solved in curwlln-

ear coordinates and discretized with one-side biased schemes.

67

4 APPLICATIONS OF OPT IMIZED UPWIND SCHEMES

To test the effectiveness of the developed 7-stencil optimized upwind scheme and the
implementation of the scheme in the multi-dimension in the presence of solid objects, a
sequence of direct numerical simulations have been carried out. Exact solutions which are
available for some test cases are calculated to validate the numerical results. Furthermore,
the standard six-order 7-stencil upwind scheme and Tam and Webb’s 7-stencil central
DRP scheme are applied to some test cases, and results computed with these schemes are
compared with those obtained using the optimized upwind scheme. In all simulations the
temporal difference scheme given in section 3.4, is used. For the simulations in which the
comparisons of various schemes are made, boundary conditions, time step, and other
parameters etc. are set to be the same. Therefore the only differences between the solu-

tions with various methods should come from the spatial discretizations.

4.1 One-Dimensional Model Wave Equation

One-dimensional model wave equation (3.1) has been intensively used to validate
numerical schemes. It is necessary for numerical schemes to pass the tests on one-dimen-
sional model wave equation before they can be further applied to more complicated prob-
lems. Two pure initial condition problems are considered here. The propagating speed for
both problems are normalized and set to be unit, i.e. c=1. The initial conditions are given

as the following:

68

case 1 u0(x) = exp[—ln2(§)z]

1 if -SOSxSSO

case 2 u0(x) = {
0 Otherwise

A Gaussian pulse and a step wave are released initially in the ﬁrst and the second cases
respectively. For the both cases the spatial step Ax and time step At are set to I and 0.1
respectively. The ﬁrst case is used to test whether the schemes are able to resolve the wave
with the wave length longer than 6Ax and the second one is used to test the ability of han-
dling discontinuities. The computational boundaries are chosen far enough so that the
waves do not reach the boundaries within the considered time.

Figure 5 shows the results at time 400, i.e. after 4000 time steps, computed with the
optimized upwind scheme, Tam and Webb’s central DRP scheme without damping and
the standard six-order accurate scheme. Figure 6 shows the computational errors caused

by these three schemes. The error is deﬁned as error = u These ﬁgures

numerical ‘ uexact'
demonstrate that the optimized upwind scheme has better accuracy near the peak than the
standard scheme derived solely by Taylor expansion, although its order of accuracy is
lower than the standard scheme. Figure 6 also indicates that the result obtained with the
optimized upwind scheme is more accurate than that obtained with central DRP scheme.
Oscillations in the computational error for the central DRP scheme is propagating

upstream, instead of being conﬁned near the peak as it is for the optimized upwind

scheme.

69

 

 

 

 

 

 

 

 

 

 

1.5 r a I . v . .
L Exact ‘
1-3 IL --- Optimized Upwind ‘
i — - - Tarn & Webb's Central DRP 1
M t 1 ---- Standard Six—Oredr Upwind 5
t .
0.9 . ’ .
i i
a 0.7 e 5
* i
0.5 ~ -
‘ i
0.3 ge -
,
0.1 — 1’ -
- —* / ‘1' " .

-0

.1 I 1 r r r 1 .
380.0 390.0 400.0 410.0 420.0
X

Figure 5 Comparison of exact and computed solutions at t=400 with different
schemes for the one-dimensional modal wave equation with an initial Gaussian
pulse.

 

 

 

 

 

 

 

 

0.12 E 3 . 5 . , § , 1
; [ ~-- Optimized Upwind
g — - — Tam & Webb’s Central DRP
0.03 ........... ---- Standard Six-Oredr Upwind ................
3 1| 3
g ; 3 1 g
' : “g ;
0 04 , ........ .......... ‘. ..... I 11.. .................................
5 I 2 r ’i I? .
- i , ' i -' ; ;
LTJ lt : I. \' \ i 1': q “1
u : //\ l‘ii/t/ 4 ,§
0 00 W~%7A+*7ﬁ’f‘ ""r't--'.---‘ .1914. i >73...“
1 lieu .*
. ' \ Ii . 1.! n ; -
\ I ' '11!” i" i
; u \ ”1 t’ '
.004 i ............................................................................. .V‘: .......... I ....................................
L!
0.0 370.0 390.0 410.0 430.0
x

Figure 6 Comparison of computational errors by different schemes for the one-
dimensional modal wave equation with a initial Gaussian pulse.

70

 

16

 

f 1
’ l Exact (a) 5
i — - - Optimized Upwind ,
l.
t

 

h

 

'1

1.
1.2

i

i

 

1.0

0.8
3 r 1
0.6 ~ 5
0.4 ~ 4

0.2

'_T

 

 

 

 

AA
0.0 v» V'

 

 

 

0.0 A 50.0 A 100.0 ‘ 150.0 ‘ 200.0 ' 250.0 ‘ 300.0
X

 

1.6 [— ' T ' T ' T ' T ' l

Exact (b)
— - — Tam & Webb's Central DRP

 

 

1.4 Ir

 

 

 

1.2 r
1.0 _
0.8 r

0.6 ~

 

 

 

1|.

. . . . 5 . .‘1 . 1

100.0 150.0 200.0 250.0 300.c
X

 

 

4

 

2 1
0.0 50.0

Figure 7 Comparison of exact solution and computed results at t=200 with (a): the
optimized upwind schemes, (b): Tam and Webb’s central DRP without damping, for
the one-dimensional modal wave equation with an initial discontinuity

71

Figure 7 shows computed results at time t=200, i.e. after 2000 time steps, with the
optimized upwind scheme and Tam and Webb’s central DRP scheme without damping. As
it is shown in the ﬁgure, the solution quality of Tam and Webb’s central DRP scheme
degrades by the presence of numerous ﬁne scale oscillations. Fine tuned artiﬁcial dissipa-
tion term is needed to remove the oscillations. In contrast to Tam and Webb’s central DRP
scheme, the solution of the optimized upwind scheme conﬁnes the oscillations near the

discontinuity due to the built in dissipation in the scheme.

4.2 Acoustic Radiation from an Oscillating Piston

xA

 

Outer
Boundary

Axis of
Symmetry

Wall

1.... E\\\\\\\\\\\\\\\\\\\\\\\\\\\\ ”

Figure 8 Schematically showing the computational domain for acoustic radiation
from an oscillating piston

 

 

 

 

72

Acoustic radiation from an oscillating circular piston in a wall is considered here to
test the behaviors of the optimized upwind scheme in multi-dimension. The radius of the
piston is 10. A cylindrical coordinate system originated at the center of piston is used. The
computational domain is set as 0 S x .<_ 100 and O S r S 100 as shown in Figure 8. The
governing equations for this case can be easily obtained from equations (2.16), (2.17) and
(2.18), which can also be written as equation (2.24) with y replaced by r and the source
term i = (—v/ r, O, 0, —v/ r) . Here v is the velocity component in r direction. The piston

is moving up and down with the speed of

= -4- 1’)
u 10 sm(5

where u is the velocity component in x direction. The grid space Ax and Ar are both set to
1. In the implementation of boundary conditions on the wall and in axis of symmetry, one
layer of ghost cells is placed near the boundaries as they are illustrated for the solid wall
boundary condition in Chapter 3. The derivatives of pressure along the directions normal
to the wall and normal to axis of symmetry can be easily obtained respectively as the fol-

lowing: On the wall, we have

 

 

 

.4
10 mcoscit) O s r < 10
5 5
3'8 = * 10-4111 111
x — =
10 cos( 5) r 10
L 0 r) 10

Since there is discontinuity at r = 10 , the derivative of pressure at r = 10 is set to be the

average of those from the left and the right to r = 10. In the axis of symmetry, we have

73

319.
37-0

Since r = O in the axis of symmetry, the source term is singular. But v = 0 too in the
av . . . . .

axis of symmetry, therefore v/ r can be replaced by— 37. .S1nce it is an axrsymmetnc prob-

lem in which acoustic waves are actually propagating in three-dimension, the radiation

condition derived for three-dimension instead of two-dimension should be used for the

outer boundaries shown in Figure 8.

The exact solution for the pressure disturbance can be written as (Hardin et al., 1995):

p: Re[eR0)I— —J (RE,)JO(r§)exp(— Jéz —(02 x — iwt)d§] (4.1)
H
where Re means the real part of a complex number, 8 = 10'4 , R = 10 , 0) = 1%, Jo and
J, are Bessel functions of order 0 and 1 respectively. The evaluation of the above integra-
tion is not trivial. Note that the integrand is singular near the point z = a) , therefore accu-
rate analytical estimation for the integration over small interval [0) — A0), a) + Am] must
be made. Since a numerical computation cannot go to inﬁnity, a ﬁnite upper bound instead
of inﬁnity is also required and should be chosen so that the error introduced by truncating
the integration from the upper-bound to the inﬁnity is within an acceptable small value. A
suitable numerical integration algorithm is needed to evaluate the integration over the
remaining region. Tam and Webb’s DRP method with and without artiﬁcial damping are
also applied to this case to compare with the optimized upwind scheme. The results shown
here are obtained at 1:200, at which the initial disturbances has propagated out of the

computational domain and the pseudo steady state is reached. The pseudo steady state

74

means that the solution is now periodic sinusoidal wave with a cycle of the piston oscilla-
tion as its period. Figure 9 and Figure 10 show surface plot and contour plot for the pres-
sure distribution obtained with the optimized upwind scheme at the beginning of a cycle
respectively. Figure 11 show the pressure contours at the beginning of a cycle computed
with Tam and Webb’s central DRP without damping. For the surface plot, regions with the
same pressure value are shown in the same color. For the contour plot, the lines with the
same pressure value are drawn in the same color. The contours obtained by the optimized
upwind scheme would have been indistinguishable from those of exact solution to an
accuracy corresponding to the thickness of the contour lines, if we had plotted them in the
same graph. The pressure contours obtained with the central DRP scheme without damp-

ing, however, have seriously spurious oscillations.

 

Figure 9 Pressure distribution computed with the optimized upwind scheme at the
beginning of a cycle

75

 

 

 

 

Figure 10 Pressure contours at the beginning
optimized upwind scheme

 

 

 

 

Figure 11 Pressure contours at the beginning

of a cycle computed with Tam and
Webb’s central DRP method without damping

76

Figure 12 shows the comparison of pressure distributions along the axis of symmetry
at the beginning of a cycle. It is shown that the pressure obtained with the optimized
upwind scheme is right on top of the exact solution. This ﬁgure further validate that the
results obtained with the optimized upwind scheme are identical to the exact solutions to
an accuracy corresponding to the thickness of lines shown. The pressure distribution
obtained with the central DRP scheme without damping, which is far away from the exact
one, are also shown in Figure 12. Our computation with the central DRP without damping
experienced that the numerical solutions obtained at small time were not too far off the

exact solution. But the solutions are getting deteriorated as the computation time

 

 

 

  
  

 

 

increases.
2 5 'x 104 l l f T l
l i Exact Solution 1 j
l 1;;- - <1 Optimized Upwind l
15 r l -——- Central DRP without Damping 1‘
7 i
9 0.5L 1
g ,1. ’ h. ,
9
0- -0.5 —‘
l” 1‘
i
—1.5 T
_2'5 . i l 1 1
0.0 20.0 40.0 60.0 80.0 100.0
X

Figure 12 Comparison of pressure distributions along the axis of symmetry at the
beginning of a cycle

77

1.0X10' j ' T ' ' W T 7

 

 

i — Optimized upwind
l — - — Central DRP without Damping
1 - - - - Central DRP with Damping

 

 

Error of Pressure

 

 

 

 

_1_0 . 1 . . i . i .
0.0 20.0 40.0 60.0 80.0 100.0

X

Figure 13 Comparison of computational errors for pressure along the axis of
symmetry at the beginning of a circle.

The spurious oscillations caused by the central DRP method without damping can be
eliminated by adding an artiﬁcial damping term. The artiﬁcial coefﬁcients, however, need
to be determined through trial and error. Fine tuning is required to get good results. Here
we follow the approach of Tarn et al. and used the same damping coefﬁcients, i.e. an artiﬁ-
cial mesh Reynolds number of 5, as it is described in the literature (Hardin et al., 1995 ).
With the inclusion of the artiﬁcially selected damping terms, spurious waves are effec-
tively eliminated in the numerical solution. If we had plotted pressure contours and pres-
sure distribution computed using the central DRP scheme with artiﬁcial damping together
with those computed using the optimized upwind scheme, they would be indistinguishable

to an accuracy corresponding to the thickness of lines shown. Figure 13 shows the com-

78

parison of computational errors for the pressure obtained with the two schemes. It demon-
strates that computational error associated with the optimized upwind scheme is
signiﬁcantly less than that associated with the central DRP with ﬁnely tuned damping
term. It demonstrates the optimized upwind scheme minimizes not only the disserpation
errors but also the dissipation errors. Furthermore, the optimized upwind scheme requires

less CPU time than the central DRP with damping.

4.3 Acoustic Propagation in a Two-Dimensional Uniform Flow

An acoustic, an entropy and a vorticity disturbances are initially released at different
locations in a uniform mean ﬂow. The uniform mean ﬂow is in the direction along x=y
with free stream Mach number of 0.5. Figure 14 schematically shows the computational
domain, boundary conditions and the locations of initial disturbances. More speciﬁcally,

the initial conditions for the disturbances are given as the following:

00:. y) = CXP{-(ln2)[(x2 + y2)/9l} + 0.1expi—(1n2)1(<x - 67)2 + (y - 605/251}
u(x, y) = 0.04(y — 67)exp{—(ln2)[((x - o7)2 + (y — 67)2)/25]}

v(x, y) = —0.04(x- 67)exp{—(ln2)[((x — 67)2 + (y — 67)2)/25 ]}

ptx. y) = expt-(lnzxoc2 + y2)/9)}

It is seen that the acoustic disturbance is released initially at the origin of the coordinate
system. The acoustic disturbance is propagating at the speed of sound relative to the mean
ﬂow. The entropy and vorticity disturbances are released initially at (67, 67). They are

propagating at the speed of mean ﬂow.

79

 

 

11 Y
Outﬂow Boundary
/
Entropy, Outﬂow
/ VQmCltY Boundary
Radiation Dlsturbances
Boundary \ 5

 

V

\

P\\\\\

Acoustic Disturbance

¥

 

 

A 45 Radiation Boundary

 

Figure 14 Schematically showing the computational domain, boundary conditions
and initial disturbances for acoustic propagation in a two-dimensional uniform ﬂow.

   

Density

Figure 15 Density distribution computed with the optimized upwind scheme at 1:80

80

The exact solution for density can be written as (Hardin et al., 1995):

 

p = 5:17] cos(gt)JO§J(x — 0.5:)2 + y2§d§ + 0.1exp{—0t2[(x — 67 — 0.5:)2 + y2]} (4.2)
0

where a1 = (ln2)/ 9 , 0L2 = (ln2)/25 and J0 is Bessel function of order 0. It also
requires numerical algorithms to evaluate the integration for the exact solution given in
equation (4.2). The computational domain is set as —100 S x S 100 and —100 S y S 100.
The griding spaces Ax and Ay are both set to 1. Note that the radii of the disturbances in
the expressions of initial conditions are chosen to test the ability of schemes to resolve the
waves with 6 PPW. Figure 15 shows the surface plot for the density distribution computed
with the optimized upwind scheme at time 80, where the acoustic disturbance partly prop-
agates out computational domain, catches up and interacts with the entropy and the vortic-
ity disturbances. No visible reﬂections near the boundaries have been observed, which
implies the outﬂow boundary conditions work quite well for this case. The density wave-
front along the line of x = y computed using the optimized upwind scheme is compared
with that computed using the central DRP method without damping, and the exact solution
(see Figure 16). The solutions are almost identical except the density obtained by the cen-
tral DRP doesn’t agree well with the exact solution at center of acoustic wave. The opti-
mized upwind scheme requires about 10 percent more CPU time than the central DRP
without damping. Figure 17 compares the computational errors of above schemes. The
error associated to the optimized upwind scheme is signiﬁcantly less than that associated
to the central DRP scheme. Figure 17 clearly demonstrates the advantage of optimized

upwind scheme over the central DRP scheme.

81

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

0.12 Y i v I v T ' T ' T ' l

—— Exact Solution 1

0.10 e — — - Optimized Upwind '

i -- - —- Central DRP without Damping n 1

0.08 — i

+ 1

0.06 1 ﬂ 1

2‘ _ .

2 0.04 . ‘J

a ii

0.02 — d

h .l

0.00 ~ A 5

-0.02 - if ‘
-0.“ 1 i r 1 1 1 A 1 * 1 L

-150.0 -100.0 -50.0 0.0 50.0 100.0 150.0

8

Figure 16 The comparison of density distributions along the line y = x at t = 80.

4.00 X 10:

 

r . r ' j ' r ' l

 

—— Optimized Upwind
1 — - — Central DRP without Damping

 

 

 

2.00 ~ 3

.. uh
-—--
-t—IIO-I-va-

— -

c- -

 

 

0.00

Error of Densny
.2
7.?—
L
at.
, as
.25

 

 

 

—2.00 e i i! -
1 i!
'1 1
_4'00 . L 1 l . 1 . 1 r 1 .
-150.0 -1oo.0 -50.0 0.0 50.0 100.0 150.0

S

Figure 17 The comparison of computational errors for density along the line y = x at
1:80

82

4.4 Acoustic Reﬂections

4.4.1 Reﬂection of an Acoustic Pulse off a Straight Wall in a Two-Dimensional
Uniform Flow

Let’s refer to Figure 14. The Mach number of the mean ﬂow is 0.5. But the angle of
attack is zero degree. The computational domain is now set as —100 S x S 100 and
0 S y S 200. The lower boundary is set to be a straight wall instead of a radiation bound-

ary. The initial conditions for the disturbances are now set as u=v=0 and

P06 y) = P(x. y) = cxr>{-(1I12)[(x2 + (y — 25)2)/25]}
In another word, only acoustic disturbance is released initially at (0, 25). The exact solu—
tion can be obtained by superposing the two mirror imaged solutions of acoustic propagat-
ing in a uniform ﬂow without a wall with two acoustic disturbances initially released at (0,

25) and (0, -25).

PreSSure

   

Figure 18 Pressure distribution computed with the optimized upwind scheme at 1:80

83

Figure 18 shows the surface plot for the pressure distribution computed with the opti-
mized upwind scheme at time t = 80, where the acoustic disturbance hit the wall, incident
and reﬂected waves propagate downstream and move out of the computational domain.
Figure 19 and Figure 20 show the comparison of computational errors for the pressure at t
= 80, along the wall and the line with x20 respectively. These ﬁgures also clearly demon-
strate the advantage of the optimized upwind scheme over the central DRP scheme. Figure
19 indicates that large computational errors encountered downstream near the lower right
corner of the computational domain. These errors are caused by the interaction of the wall
and outﬂow boundary conditions. Further study on accurately imposing boundary condi-

tions at the comer of computational domain is needed to reduce the computational error.

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

2.00 X104. . - 1 . .
-— Optimized Upwind
— - — Central DRP without Damping
1.00 c
2
3 t
(I)
(I)
‘2
l 0.00
‘5
§
Lrl
-1.00 ~
-2.00 . E J * * ‘
-100.0 -50.0 0.0 50.0

X

Figure 19 The comparison of computational errors for pressure along the wall at t =
80

84

1.00 x 10-3 .

 

 

 

 

 

 

 

 

 

 

—— Optimized Upwind
- - — Central DRP without Damping
0.50 ~ /~ 1‘ 5
93
3
a)
U)
9
0- 0.00 ' ~
“5
e
m
-0.50 a
__1 00 . l . 1 1 1 .
0.0 50.0 100.0 150.0 200.0

Figure 20 The comparison of computational errors for pressure along the line x=0 at
r=80

4.4.2 Scattering of an Acoustic Pulse off a Cylinder

Unlike the previous applications which are limited to Cartesian coordinates system
with straight boundaries, this example is used to study the behaviors of the optimized
upwind scheme in the presence of a curved object. Figure 21 shows the diagram of noise
source, computational domain and boundaries. The noise source is located at 4 cylinder
diameters away from the origin in the x-axis. A uniform ﬂow is assumed in the z-direction
which is pointing out of the paper. Therefore a two-dimensional still ﬂow can be assumed

in the x and y plane. This example can be associated to the physical problem of ﬁnding the

85

sound ﬁeld generated by a propeller scattered off by the fuselage of an aircraft. The fuse-
lage is idealized as a circular cylinder, and the noise source is given at a speciﬁed location.
Since the problem is symmetric about x-axis, the computations are conducted only for the
upper half plane. Computational domain is chosen as a semi-circular region which is
bounded by an outer semi—circle, the cylinder wall, and symmetric lines. The meshes are
generated by the circumferential and radial lines.The polar coordinate systems is used to

transfer the acoustic ﬁeld equation in the Cartesian system as it is described in Chapter 2.

lly

  
 
 

  

/

outer boundary

 

 

axis of symmetry fuselage noise source

Figure 21 Schematic diagram showing the computational domain and boundaries of
acoustic scattering problem

TWO cases are considered. In the ﬁrst case a pressure pulse is initially released and

given by

We y) = expi-(Inzllor-4)2 +y2)/0.041}

86

In the second case a noise source is added to the acoustic energy equation. The right hand

side in equation (2.18) is replaced by the following periodic source term

p(x, y) = exp{—(ln2)[(x — 4)2 + y2)/0.04]}sin(81tt)
Wall and outer boundaries are treated by the procedures described in section 3.5. For grid
points near the symmetry lines, we have
u(x, y) = u(x, -y), V(x, y) = —V(x, —y), We y) = p(x. -y)

Since there are only acoustic waves in this case, the density and pressure are identical. A
grid system of 200 x 180 points with Ar = 0.0375D and A0 = 10 is used for the simu-
lation in the ﬁrst case, where D is the diameter of the cylinder. Time step is chosen as
0.001. Figure 22 shows the surface plot for the pressure distribution at t = 7 for the ﬁrst
case. Figure 23 shows the pressure as a function of time for three different locations A, B
and C. These ﬁgures clearly show the incident wave, reﬂection wave and the effect of wall
on the incident wave. A grid system of 400 X 360 points with Ar = 0.025D and
A0 = 0.50 is used for the simulation in the second case. Time step is chosen as 0.0005.
Figure 24 shows the scattered pressure wave patterns at t = 40, in which the effect of the
cylinder on the scattered sound is clearly demonstrated. Figure 25 gives the computed
directivity function D(0) = r1.)a at r= 9. The computed results are obtained by marching
time to t = 60 at which a pseudo steady state is achieved. Since we do not have the exact
solutions, the accuracy of the numerical results are unable to be evaluated. Our results
qualitively agree well with those reported in the literatures (Kim et al., 1997, Tam et al.,
1997). To obtain high quality results with Tam and Webb’s central DRP methods for this
problem, artiﬁcial damping with different damping coefﬁcients for different region is

required as it is pointed out by Tam et al (1997).

87

Pre55u re

   

Figure 22 Pressure distribution computed with the optimized upwind scheme at t = 7
for the ﬁrst case in which an acoustic pulse is released initially.

 

 

—— PointA
-—— Point B .
‘ —-- PointC """""""

 

 

 

 

 

 

 

Pressure

 

   

 

 

 

 

-0.02 L V j

-0.04 e . - . .
6.0 6.5 7.0 7.5 8.0 8.5 9.0 9.5 10.0
Time

 

Figure 23 Time history of pressure ﬂuctuation at points A (r = 5, 0 = 900), B
(r = 5,9 =135") and c (r = 5,0 =180") as indicated in Figure 21.

88

///////%///;/W ,,

lllllli

1W ‘

 

Figure 24 Pressure wave pattern at t = 40 simulated by the optimized upwind
scheme

 

3.0

iii/111 . A
ll

 

 

 

 

 

D(0)

 

 

 

\/
E31V

0.0

 

 

 

 

 

 

 

90 100 110 120 130 140 150 160 170 180
0(degree)

Figure 25 Directivity of radiated sound computed by the optimized upwind scheme

89

4.4.3 Reﬂection of an Acoustic Pulse off a Sphere

This problem is almost same as the ﬁrst case in the scattering of acoustic waves off a
cylinder except that the object is a sphere. The computation is carried out with the spheri-
cal coordinate transform to test the generalization of the optimized upwind scheme to
acoustic ﬁeld equations in three-dimensional curvilinear coordinate system. A pressure
pulse given below is initially released at 2 sphere diameters away from the center of the

sphere without the mean ﬂow

p(x, y) = —0.01exp{-—(ln2)[(x — 2)2 + y2 + z2)/0.252]}
Since the transform of the Cartesian coordinate system to the spherical coordinate system
has singularity at poles, special treatments for the grid points near the poles are needed for
the numerical stability. Let (r, 0, (p) be the spherical coordinated system, the transform

can be written

x = rsinGcostp
y = rsin0sintp
z = rcosO

The three-dimensional acoustic ﬁeld equation (2.27) without source term becomes:

&E’] + (Aosinecosrp + Bosinesintp + Cocos0)§a;-3 (4-3)

+ %(A0cosecostp + Bocosesintp—Cosin0)aa—O?1

 

. 3 -) _
rsin9(—Aosmtp+ Bocosrp)%q — 0

Therefore for grid points at the north pole where 0 = 0 , the above equation is not suitable
because the denominator in the last term in left hand side is 0. For the north pole, we go

back to use the equation (2.20) in Cartesian coordinate system. The derivatives are evalu-

90

ated by the following:

Qv

8421:1323 aj=1§§ and 9.3.:

ax ra_9(p=o’ 8y ;BO¢=n/2 az

QJIQJ
‘1

where the subscripts denote the directions or planes which are used to evaluate the deriva-
tives. Since the grids are screwed near the north pole, it impose very strict limit on the time
step. In the calculation, the values of acoustic variables in ﬁrst three levels of grid points
near the north pole are obtained by third-order accurate interpolations of the values at the
pole and the values at the fourth and ﬁfth level of grid points surrounding the pole. By
using the properties of symmetry, only the region bounded by 0 S 0 S tt/ 2 and 0 S (p S 1:
is needed for the simulation. A grid system of 75 x 45 x 60 points with Ar = 7D/30,
A0 = 20 and Atp = 30 is used. Time step is chosen as 0.001. To run this case in SGI
INDY R5000 up to t = 2.5, it requires about 11 CPU hours. To accurately evaluate some
properties like the pressure history at different points or the directivity for a periodic
source with the frequency given in the scattering of waves off a cylinder, grid system of
points at least double as many as the points presently used in each direction is needed.
Taking into account of the smaller time step required for ﬁner grid, it will need at least
about 176 CPU hours, or a week in SGI INDY R5000 for a decent result. Therefore simu-
lating the real life three-dimensional acoustic problems requires the use of super computer.
Figure 26 and Figure 27 show the surface plots of pressure patterns on the surface of the
sphere and various planes of x=0, y=0 and 2:0 at t=2.5 computed by the optimized
upwind scheme. They clearly demonstrate the radiation waves, reﬂected waves and the
interaction of waves and the sphere. There are not noticeable differences in the wave pat-

terns near the north pole, which validate the treatments for the grid points near the pole.

91

 

Figure 26 Pressure patterns on the surface of the sphere and the plane of z = 0 at t =
2.5.

 

Figure 27 Pressure patterns on the surface of the sphere and the planes of x=0, y=0
and z=0 at t=2.5

5 F INITE DIFFERENCE VERSUS F IN ITE VOLUME
DISCRETIZATION S

Two discretization approaches are widely used in CFD, namely ﬁnite difference and
ﬁnite volume schemes. In a ﬁnite difference approach the governing equations are ﬁrst
transformed from the physical space (x, y, z) to a computational space in the general curvi-
linear coordinate system (g, n, C). Then the derivatives in the governing equations are
approximated and discretized with ﬁnite differences in the computational domain. In a
ﬁnite volume approach, the grid cells are treated as control volumes, and the governing
equations are integrated directly in the control volumes to yield the discretized equations.
It can be shown that with proper interpretation of metric terms in a ﬁnite difference, the

ﬁnite difference scheme is identical to a ﬁnite volume scheme.

5.1 Navier-Stokes Equations

In this study the three-dimensional unsteady Navier-Stokes equations in Cartesian coordi-

nate system are employed, which can be written in conservation form as the following:
a a B F 9 8 9
até + —(dx 1% Ev) + By( Fv) + az(G 3 ) 0 ( )

where a is the conservative ﬂow variable vector, 2, F and G are the inviscid ﬂux vectors
1n x, y and 2 directions respectrvely. v, v and CV are v1scous ﬂux vectors correspondlng

to x, y and 2 directions respectively. In detail, we have

93

Di

and

 

pu
pv
pw

e
PrJ

 

(11¢

 

T f
P“ P"
pu2 -l- p puv
2
puv gr? =1 pv +p
puw pvw
pu(e, + 2) pv(e, + B)
p 1 k p
0
1:xx
gv = TX),
sz
MT '1' VT '1' WT -Ka—1:
xx xy xz ax
0
Ty):
2v = Tyy
Tyz
0T
“Tyx + v1:yy + WTyz—K-é;
0
1:2):
—)
Gv = TZ)’
1:22
GT
u’tzx + vrzy + wtzz—x-a;

 

 

 

 

 

94

 

 

l

 

J

 

J

Q1

 

pw
puw
pvw

2
9w +p

pw(e,+g)

l (5.2)

 

(5.3)

(5.4)

(5.5)

where the symbols used here for ﬂow variables have the same meanings as those used in
acoustic ﬁeld equations. The new notations e, and K represent the total energy per unit
mass and the thermal conductivity. The total energy can be obtained from the expression
of e, = e + 0.5(u2 + v2 + wz). The relations of viscous stresses to the rates of strain, the

velocity gradients, are given as the following:

In = — £11041, + vy + W2) + Zuux
tyy = - gum, + V), + Wz) + Zuvy
I“ = - gmux + V), + W2) + 211.11)z
I” = Tyx = u(uy+vx)
Ta = In = “(“z‘l'wr)
Tyz = IO = u(vz+wy)

Where it is the dynamic viscosity.

 

S
-—-.I

 

 

 

 

 

 

 

: i ’1 i :

I . I, ' ' l
A i ;_:.—-E.--;:_,._-_I-_-1

l---:-L- '/--: ---i
§=§(xIsz) ’1 I ' ... :- ,I

' i 1'.-.LI:-L.]’!‘., .I’--L--

n=n(x.y.z) "T'T"1';?'“C :

L'----.1---L----_--'
§=§(x,y.z) ; g
T: t n

 

Figure 28 Coordinate transformation and the grid cell around grid point (i, j, k)

The Navier-Stokes equations can be transformed from the physical space (x, y, z, t) to

curvilinear space (§, 11, C, ‘6) shown in Figure 28 and written as the following:

95

2.

at“; - 0v): 0 (5.6)

a-- a - ~ 3 ~ ~
5;Q+-a—§(E—Ev)+ﬁ(F-Fv)+

~ -1 -> . . . . . . .
where Q = J Q. Inv1sc1d and VlSCOUS ﬂux vectors 1n the curv111near coordinates are

is: 14095:, + 3% + 3&2) (5.7)
17‘: J-l(i)3nx + 2n, + 311,) (5.8)
5;: 14(ng + Iggy + GE.) (5.9)
135,: J"1 (13:, + ﬂay + avgz) (5.10)
15v: 1'1 (fimx + Any + 3.11,) (5.11)
6,: J“ (qux + icy + 3.4,) (5.12)

Note that the conservation form of Navier-Stokes in the Cartesian coordinates system is
preserved after the coordinate transformation. In another word, equation (5.6) is in the
conservation form in the genera] curvilinear coordinate system. J is the determinant of
Jacobian matrix of the coordinate transform. In this study only stationary grids are consid-
ered. Therefore {3, n and C are functions of x, y and 2 only. They are not function of time
t. The reciprocal of the determinant of the J acobian matrix of the coordinate transform can

be explicitly written as:
—l
J = xtynzt + x§yézn + xnyczz ‘ xéyczn“xnyézc ‘ xcynzé (5-13)

96

Let A , f3 and E? be the three Jacobian matrices of [:3 , F and E} with respective to Q, i.e.,

24:35 i=8: 53:95.3 (5.14)
8Q 8Q 8Q
Then they can be diagonalized as
A: igxgig‘ I9: inﬁnigl 6: Tciigigl (5.15)

where 11g, 1111 and ﬁg are the diagonal eigenvalue matrices which are the same as those
for the acoustic ﬁeld equations. T5,, in and TC are the corresponding eigenvector matri-
ces. We should note that .21 , E , 6‘ and T5, Tn, TC are different from those for the acoustic
ﬁeld equations. In the acoustic ﬁeld equations, we do not consider the conservation forms
which are necessary for the computation of the mean ﬂow when the discontinuities like
shock wave appear. The detailed formulae of the matrices A , l9 , 6 and T5, Tn, T; for the
conservative Navier-Stokes equations in the general curvilinear coordinate system can be

found in the literatures (Pulliam et al., 1980, Warming et al., 1975, 1978).

5.2 Finite Difference Versus Finite Volume Discretizations

Finite difference approach starts from the differential equation (5.6) and seeks the
numerical approximations to the derivatives appearing in the equation. The optimized
upwind schemes developed in the study on CAA belong to this approach. For example, the
partial derivative in the § direction at grid point (i, j, k) can be discretized by a ﬁnite dif-

ference with up to second order accuracy and written as the following:

97

9—(3 — 13".)

BE; (ff—Fr)”1/2,j,k-(l:5—l:§v)i-1/2.j,k (5.16)

Ill

where E and Ev are the ﬁnite difference approximations to if and 53,, at surface (i+I/2, j,
k) respectively. Without loss of the generality we have assumed that the grid sizes Aﬁ , A11
and AC in the computational domain are all set to l.

A cell-centered ﬁnite volume approach starts from the following integral equation

I &QdV + “that! —Z'v) + ny(z‘— 122v) + rig-(I)? — 3v)}ds = 0 (5.17)
S

i,j.k

where v, j, k is volume of grid cell with index (i, j, k) and 3 = (nx, n nz) is the outward

y,
normal vector at the cell surface S. The surface S of cell (i, j, k) shown in Figure 28 is com-
posed of six faces in i, j and k directions, with two faces in each direction. Equation (5.17)
can be derived by integrating equation (5.1) over grid cell (i, j, k) or directly from the
Navier-Stokes equations in the integral form. Finite volume approach seeks the numerical

approximations to the integrations over the cell surfaces. For example, for cell (i, j, k)

shown in Figure 28, the integration on surface (i+ ”2, j, k) is approximated by

[MAE—23:.) + nydi-fv) + nz(?}— 3.)}ds (5.18)
S
_=_ {S[nx(l73—E‘v) + ny(F—Fv) + nz(Z';— 6,)1},+,/2,j,k

where E and E“, are the numerical approximations to E and g1; at the surface respec-
tively. Similarly Ff}, and 3,5,, are the approximations to F, F v and 3, av at the same
surface respectively. S in the right hand side of the above equation is the area of the sur-

face.

98

Refer to Figure 28 and consider the cell surface between the grid points (i, j, k) and
(i+ I, j, k). In the ﬁnite difference approach, the metric terms can be evaluated in the fol-

lowing way,

—1

Ji,j,k = Vi,j,k (5.19)
(J-IE-ix)i+l/2,j,k = (snx),-+1/2, M (5-20)
(1"i,>.../2,,-,. = (Sm-Writ (521)
(149,“), M = (5312),.“ a“ (5.22)

Substituting above metrics into equation (5.16), we get

(ﬁ—év)i+l/2,j,k = {.<i[n,(1‘§-—E,)+ny(F—F,)+nz(f‘;—'c';,)]},.+1 ,2, j, k (5.23)

which is the same as the ﬁnite volume discretization at the cell surface given in the equa-
tion (5.18). For the ﬁnite difference approximation at surfaces in other directions, we can
similarly show that the ﬁnite difference approximation is equivalent to the ﬁnite volume
approximation provided the grid metrics in the surfaces are evaluated in a way similar to
equations (5. 19) to (5.22). If the cell shown in Figure 28 is deﬁned as a control volume for
grid points (i, j, k) and metric terms are evaluated by formulae given in equations (5.19) to
(5.22), then the ﬁnite difference approach is identical to a ﬁnite volume approach.
Although OVERFLOW is a ﬁnite difference based code, the metric terms are calculated in

a pseudo-ﬁnite volume manner. If the control volume for the metric terms around each

99

ﬁnite difference grid point are unambiguously deﬁned in OVERFLOW, the ﬁnite differ-
ence method in OVERFLOW can be viewed as the same as a cell-centered ﬁnite volume
method on the control volume mesh. The fully conservative Chimera, which was devel-
oped originally based on ﬁnite volume scheme, can then be naturally applied to the ﬁnite

difference scheme and implemented in the OVERFLOW.

5.3 Finite Difference Grid and Finite Volume Grid

As it is mentioned before, the OVERFLOW is a ﬁnite difference based code. There-
fore the computational grid is interpreted as a ﬁnite difference grid. To implement conser-
vative Chimera in the OVERFLOW, the control volume for each ﬁnite difference grid
point needs to be deﬁned unambiguously and relevant metric terms should be calculated
based on the deﬁned volume. Since the information on the patch interfaces and other rele-
vant information between the overlapped grids are needed, an independent module called
the Zonal Interface Generator (ZIG) is developed by Wang in CFDRC. The ZIG is also
based on a ﬁnite volume grid. Therefore another grid, namely the ﬁnite volume grid corre-
sponding to the ﬁnite difference grid, is required.

In the implementation the given input grid system which is generated by some other
codes, for example, CFD-GEOM from CFDRC, is treated as a ﬁnite volume grid. The
enhanced ﬁnite difference grid is then computed as the mutual grid. It is formed by simply
connecting all the centroids of control volumes, boundary face centers and the comers.
Figure 29 schematically shows the concept in the two-dimensional sense. In such a
approach, all the geometric quantities are computed based on the given original ﬁnite vol-

ume grid which guarantees the accurate representation of geometry conﬁguration. Thus

100

the given original ﬁnite volume grid is input to the conservative interface generator (ZIG)
to generate patch interfaces and calculates the relevant information for the conservative
Chimera. This grid is also used by OVERFLOW with the conservative Chimera when
determining cell volumes, cell face areas and other relevant metric terms. The enhanced
ﬁnite difference grid is directly input to the OVERFLOW to allocate dynamic memory and

to apply boundary conditions.

 

 

 

 

 

 

I I I I I I
I I I I I I
I I I I I I
a... ooooooooooooo I ..- . nun-.unocb uuuuuuuuuuuuuuu
I I I I
I I I I I
I I I I I I
I I I I I I
I I I I I I
I I I I I
..., ooooooooooooooooo . ........ ' ........ 'IIIJ-IC‘CIC
I I I I I I
I I I I
I I I I I
I I I I I I
I I I I I I
I I I I I I
I I I I I
--.:-..L-.-'--uuduou.--- III-.--.-----:--.dIon-...-
I I I I I
I I I I I
I I I I I I
I I I I I I
I I I I I
I I I I I I
..., ....... . ..... O ..000. I 0.....-.-.' ....... ‘Ioo
I I I I I
I I I I I I
L I I l I
I I I I I I
I I I I I
I I I I I I
I I I I I I
on... ....... . --------- . ........ .-.-.--.-:--------.---
I I I I I I
I I I I I I
I I I I I I
I I I I I
I I I I I I
I I I I I I
cant ------- 'n-nodonuu. ........ . ........ I ....... J...
I I I I I I
I I I I I I

 

 

 

 

 

 

 

 

Figure 29 Schematically showing ﬁnite volume grid (solid line) and mutual ﬁnite
difference grid (dash line)

101

6 IMPLEMENTATION OF CONSERVATIVE CHIMERA IN
OVERF LOW

6.1 Review of the Chimera Approach

As it is described in section 1.2 on the background of the Chimera, there are many
blocks in the computational domain. The grids in the neighboring zones are overlapped
and the information for variable values in the neighboring zones is needed to exchange.
Figure 30 schematically shows overlapped grids in two neighboring zones in two-dimen-
sional sense. As it is shown in the ﬁgure, the boundary of one grid intersects with the inte-
rior cells of the other grid. They are usually referred as the minor and major grids as it is
indicated in Figure 30. A grid can be a minor grid of several major grids. Reversely, a grid
can be a major grid of several minor grids. The boundaries of the minor grids which cut
into the interior cells of major grids usually form the patch interfaces.

In the traditional Chimera approach, for each cell near the patch interfaces in the minor
grid, the interior cells in the major grid surrounding the cell in the minor grid near the
patch interfaces are ﬁrst located. The values of ﬂow variables in the cell in the minor grid
are obtained by the interpolation of the values of ﬂow variables in the surrounding cells in
the major grid. The values in the cell R in the minor grid in Figure 30, for example, may be
obtained from the values in cells enclosed by the rectangle indicated by the thick lines in
the major grid through the interpolation. The required number of surrounding cells
depends on the required accuracy of the interpolation. The ﬂuxes ﬂowing into the minor

grid through the patch interfaces may not be in the same amount as those ﬂowing out of

102

the major grid through the interfaces. Therefore globally it results in the non-conservation.

 

 

 

 

4 34
L. '9’ . .....
I b.'\2 R .3 hhhh
a ' if .

 

' ~ . ’ Minor Grid

.' .
I~Q

. .‘
I

 

 

 

 

 

Major Grid

 

 

 

Figure 30 Schematically showing cut cells, smallest faces and patch interface

In the fully conservative Chimera approach developed by Wang et al (1994), the patch
interface is set to be a common interface of the minor and major grids. Cells of the major
grid which completely inside the minor grid are marked as hole cells. The ﬂow solver
skips all hole cells through a marker. Cells in the major grids which are cut by the patch
interface are named as cut cells, for example, the cell L shown in Figure 30. Cells in the
major grid that are outside the region covered by the minor grid are called normal cells.
Now the cell surfaces in the cut cells are reformed. Refer to the cut cell L showed in Figure

30, cell surfaces consist of faces I-b, b-c, c-3, 3-4 and 4-1. Faces like I-b and c-3 which

103

are cut from regular cell faces 1-2 and 2-3 by the patch interfaces are referred as cut faces.
The cell surfaces for the cell in the minor grid near the patch interface are also cut into
smaller faces by the cut cells in the major grid. For cell R in the minor grid in Figure 30,
for example, its cell surface labeled with a-d is divided into smaller faces a-b, b-c and c-d
which are referred as the smallest faces. With the ﬁnite volume approach the net ﬂuxes
through a cell are the sum of numerical ﬂuxes through cell faces. Therefore in order to pre-
serve the conservation globally, the amount of the ﬂuxes ﬂowing into the minor grid and
ﬂowing out of the major grid through the patch interfaces must be the same. Wang et al
(1994) proposed to let the contributions of the ﬂuxes on the smallest faces to the major
grid and minor grid be exactly the same. Refer to the smallest face b-c, we add the amount
of ﬂuxes on the smallest face to the net ﬂuxes of its left cell L and substrate the same
amount the of ﬂuxes from the net ﬂuxes of its right cell R. The remaining problems is how
to numerically evaluate the ﬂuxes on the smallest faces and cut faces with the desired

accuracy which is usually of second order in CFD.

6.2 Numerical Fluxes through the Patch Interface

Since the enhanced ﬁnite difference grid is constructed based on the given ﬁnite vol-
ume grid, a control volume for each ﬁnite difference grid point is precisely deﬁned. The
metric terms for each ﬁnite difference grid point are calculated based on the control vol-
ume for that point, which are calculated exactly in the way shown in section 5.2. In the
major grid, the cut face areas in each cut cell and each hole cell are assigned to zero ﬁrst.
Then the actual cut face areas and the actual volume in the each cut cell are updated using

the information generated by the module of ZIG. From equation (5.18) for the ﬁnite vol-

104

ume approach or equation (5.23) for ﬁnite difference approach with the grid metrics com-
puted as those illustrated in section 5.2, the amount of ﬂuxes in cell surfaces are
pr0portional to the areas of the surfaces. Therefore the correct amount of ﬂuxes through
the cut faces in each cut cell is automatically calculated in the OVERFLOW, and the ﬂuxes
through hole cell faces are zero because of zero face area assigned at the beginning. For
the cut cell L in the major grid shown in Figure 30, for example, the ﬂuxes through the cell
face 1-2 is equivalent to the ﬂuxes through the cut face I-b. Thus the total ﬂuxes out of this
cell calculated in the OVERFLOW, the sum of the ﬂuxes through four cell faces 1-2, 2-3,
3-4 and 4-1, are just the sum of ﬂuxes through out faces 1 -b, c-3 and normal faces 3-4, 4-1.
In the minor grid, the areas of the cell faces corresponding to the patch boundary, the face
a-d in the cell R shown in Figure 30, for example, are also set to zero at beginning. The
total net ﬂuxes out of the cell R obtained in the OVERFLOW are the total ﬂuxes through
a-f, f-e and e-d.

Therefore the remaining work to correctly compute the net ﬂuxes out the cut cell L is
to evaluate the numerical ﬂuxes out of the cut cell through the smallest face b-c. The
remaining work to correctly compute the net ﬂuxes out the cell R is to evaluate the numer-
ical ﬂuxes out of the cut cell through cell face a-d which consist of the ﬂuxes out the cell R
through the smallest faces a-b, b-c and c-d. To achieve ﬂux conservation for the overlap-
ping Chimera grids, the same amount of ﬂuxes across each smallest face with different
signs are added to the total net ﬂuxes for the cells sharing this face in both the major and
the minor grids. Refer to the smallest face b-c, let jibe be the numerical ﬂuxes through the
smallest face b-c with the face normal pointing to the cell R as it is shown in Figure 30.

ibc is added to the net ﬂuxes of cut cell L computed from the OVERFLOW and —1~Tbc is

105

added to the net ﬂuxes of cell R computed from the OVERFLOW. Therefore to implement
the fully conservative chimera, the remaining work is to evaluate the numerical ﬂuxes on
all smallest faces like b-c. The computation of the ﬂuxes through each smallest face is
done separately in the implementation of the fully conservative Chimera in OVERFLOW.
Consider the smallest face b-c shown in Figure 30. Let 61. be the variable vector at the
center of cut cell L and QR be the variable vector at the center of the right cell R. Then the
inviscid ﬂuxes through the smallest face b-c is determined by Roe’s Riemann solver (Roe,

1983). To be more speciﬁc, denote

iti’zr, 3) = Sinjstét) + 22,3151.) + nﬁtfzm (6.1)
26312.3) = S[an*(§R) + afﬁx) + nziidim (6.2)
m5... 3) = S[nxA(Qa) ”ya—£50) + nzcéan (6.3)

where S is the area of the smallest face b-c (the arc length in two-dimension), at! is the

_)
Roe’s average of Q1. and 5R , Jacobian matrices of the ﬂuxes A, B and C are

m
tlm

-> ->
14:95 B=a—F C:
36 35 do

The inviscid ﬂuxes through the smallest face b-c are approximated by the following

-> ~ 1 ->
Jin,E + a} + n.31ds = $1.. = 51 ﬂér. 3) + rr'ér. 3) - Iméa. inkér — on] (6.4)
s

where IAI = TIAIT-l with T the eigenvector matrix, A the diagonal eigenvalue matrix.

The convention for the absolute value of a diagonal matrix is the same as it is deﬁned for

106

the acoustic ﬁeld equations.
Refer to Figure 30 again. To calculate the viscous ﬂuxes through the smallest face b-c,

the following integral is evaluated

-> —> _ — —
j {11ng + nva + nva}dSE S[an,, +11va + nva] (6.5)
S

Viscous ﬂuxes Ev , RV and 3. consist of laminar and turbulent eddy viscosities, ﬂow vari-
ables and velocity and temperature gradients Vu , Vv , Vw and VT. All quantities on the
smallest face b-c in the viscous ﬂuxes except for the velocity and temperature gradients
are calculated as the average of the quantities in the cut cell L and the cell R. The key to
evaluate the viscous ﬂuxes through face b-c is to ﬁnd the gradients of u, v, w and Ton that
face. The thin layer approximation is used in calculating the ﬂow gradients on face b-c.

More speciﬁcally, under such an approximation, we have

(Vu),,c = 7 it (6.6a)
(Vv),,,. = gv—n in (6.6b)
(Vw),,,, = 9L: 3 (66¢)
(VT),,,. = a_: Z (6.6d)

In another word, we assume that the ﬂow variables do not change along the surfaces. This

approximation is usually applied in evaluating the gradients on surface when the schemes

107

up to second order are employed. Let us only consider velocity u in x-direction to illustrate
the procedure of evaluating the gradients. Neglect second order and higher order terms in

the Taylor expansion, we have

(uR— age—(V141,. - (ii-i1.) (6.7)

Where it, = (xL, yL, 2L) and irg = (xR, yR, zR) are the position vectors in the cell cen-
ters of the cut cell L and the cell R respectively. u L and u R are values of u evaluated at the
cell center of the cut cell L and the cell R respectively, which are both available at the cur-

rent time stage. Substituting equation (6.6a) into equation (6.7), we obtain

Bu = A1. (6.8)

87 lithe-h)

Substituting equation (6.8) back into equation (6.6a), we get the gradient of u on face b-c

as the following

(6.9)

 

= 3‘(’)’R-;L)

Vv , Vw and VT can be evaluated in the same way.

Since the cell cutting in the 216 is arbitrary, the direction of iR - :1, may be far away
form the direction it. The angle between i): and iR — i‘L may even be greater than 90
degree, i.e. i): - (iR — fl) is negative in some extreme cases. Then the velocity gradients
derived in equation (6.9) can produce large errors. In these cases we assume that the veloc-

ity gradients are in the direction of (ig— h) , i.e. the assumption in equation (6.6a) is

108

replaced by

 

(i —i )
(Vuh. = (Vania-Le"— <6-10>
Irr - rtl
Repeating the procedures of deriving the gradients, we now have
> 9
u — u r — r
(Vu)bc = ( L R)( R L) if ﬁ-(iR—hmo (6.11)

lie-i112

This special treatment is found to be necessary for boosting the stability and work very

satisfactorily in the tests which are shown in next chapter.

6.3 Implicit Treatment of the Patch Interface

Since a ﬁnite difference scheme can be identical to a ﬁnite volume scheme, in this sec-
tion we are using the ﬁnite difference approach to illustrate the implicit treatment of the
patch interface for convenience. After the spatial derivatives are approximated by the ﬁnite

differences, equation (5.6) in grid point (i, j, k) becomes

3 " _ 9. .
(E )1,“ _ —R.,,,r (6.12)

where the right hand side of the above equation can be formally written as

Rt, j,k = (1:3 - Ev)i+ 1/2, j, k — (1:? — 1?):- 1/2,j,k (6-13)
+(ff-Ev)”1/2,j,k—(E-Ev)i—1/2,j,k
+ (CT—Ev)”1/2,j,k—(5-5v)r—1/2,j,k

109

The notations associated to the ﬂuxes in the § direction are the same as they are deﬁned in

section 5.2. The notations associated to the ﬂuxes in the other directions are deﬁned simi-

larly. F and Fr are the ﬁnite difference approximation to F and F, respectively. 6 and

Or are the ﬁnite difference approximation to 6} and 6}, respectively. From now on, when-

ever the meanings are clear, we will omit the subscripts and superscripts for simplicity.
Since only the steady ﬂow is interesting in this study, considering ﬁrst order accurate

backward ﬁnite difference for the time derivative in equation (6.12). we obtain

 

~n+l ~n -)n+1 ->n -)n+1 -)n

Q —Q = —AtR = —AtR —At(R — R) (6.14)
Now utilizing

-)n+1 -)n -) ~n+1 9 ~n a} ~n+1 ~n

R -R =R(Q )-R(Q )=—. (Q -Q) (6.15)

8Q Q"
we have
(magazré = —AtRn (6.16)
3Q

where I is the identity matrix and AQ = Q" +1 — Q". Since only steady ﬂow is consid-
ered, the left hand side of above equation is not necessary to be high order accurate. As it
is seen that AQ approaches to zero as the ﬂow approaches to the steady state. Therefore
the left hand side (LHS) vanishes when the steady state reaches. The important thing is to

ﬁnd an approximation matrix to 33/36 so that the resulted linear algebraic equation

110

(6.16) can be easily solved while the matrix greatly boosts the stability and convergence
rate. One method implemented in the OVERFLOW is to use the central difference for R
in the left hand side of equation (6.16) to get the approximation matrix. Let us consider

inviscid ﬂows here for easy illustration, notice that

~ ~" ~

 

 

Eté"*‘)—E(é")e§ AQ = A AQ (617a)
do 6"

F(é"“)—F(é">=i A6 = FAQ (6171:)
do 6"

é<é"“)—6(é")z3_9 Ar”? = &"AQ (me)
at). d"

We have
[I+At(8§An+5nR"+5§Cn)]AQ = 45:13" (6.18)

where 8g, 5“ and 8C are the central difference operators in §, 1] and C directions respec-
tively. More speciﬁcally, the results of the operation of the central difference operator 5g

in § direction on AAQ at grid point (i, j, k) can be written as
.. .. I ~ ~ ~ ~
SEMAQ) = 5M” 1/2.j.kAQi+1,j.k -Ai—1/2, j, kAQt+ 1,j,k] (6.19)

Where A” ”2, M is A evaluated at Q,- + In, L k which is some kind of average of Q), 13k
and Q” 1, j, k at the surface (i+I/2,j,k). The linear system of algebraic equation (6.18),

however, is not easy to solve. It involves 3 grid points in each of L5,, 11 and Q directions.

111

Combining all grid points together, we need to solve linear system of 7-band diagonal
block matrices. Each block is a 5 x 5 matrix in three-dimension and a 4 x 4 matrix in
two-dimension. It will require huge amount of computational time. The splitting technique
is used to reduce to computational work in solving the linear system. Equation (6.18) is

split as

~n ~n ~n ~ -)n
(I + AtﬁgA )(I + AtBnB )(I + AISCC )AQ = —AtR (6.20)

The above equation requires solving the linear algebraic system of tri-diagonal block
matrix in each direction, which greatly reduced the computational work. Using the proper-
ties of Jacobian matrices in equation (5.15), we can approximately write equation (6.20)

as

igu + meg/Kai? inu + AtﬁnAn)T,—11T§(I + At8§A§)TEIAQ = —At it" (6.21)

Now the each block is only a diagonal matrix instead of a 5 x 5 or a 4 x 4 matrix. As a
matter a fact, it requires to solve tri-diagonal system of equations in each direction which
can be easily and effectively solved by the special algorithm developed for tri-diagonal
system of equations. If the viscosity terms are considered and smoothing terms for stabil-

ity are added, then equation (6.21) is replaced by (Pulliam et al., 1980)

ha + AtﬁgAg + Vlsg + SMogﬁgl in (6.22)
(I + AtﬁnAn + vrsn + SMOn)T,-]1T§n
(I + AtSCAg + VIS; + SMOQTEIAQ: 4:13"

This is so called the diagonal approximate factorization algorithm implemented in the

112

OVERFLOW as one method to march the solution to steady state. The viscous terms and

the smoothing terms can be respectively written as
2 2 2
VIS§= AthAgﬂlL + uT)(§x + Cy + CZ )I (6.23)

SMog = At[—V§A§ez+(V§A§)ze4]I (6.24)

Where 11L and pr are the laminar viscosity and turbulence eddy viscosity, A: and V: are

forward and backward-difference operators, e. g.
Agémﬂ: = éi+l,j,k-éi,j,k

Vgéi,j,k = Qi,j,k- Qi-1,j,k
£2 and 84 are the second order and the fourth order dissipation coefﬁcients. The second
order dissipation coefﬁcient is proportional to the maximum of absolute eigenvalues. We
should note here that AC, A11 and AC have been set to 1 for convenience.

Without any special treatment in OVERFLOW, ﬂux Jacobians on the smallest faces
are not included in the LHS implicit operator, although the ﬂux Jacobians on cut faces are
actually calculated correctly. Since the cutting in the zonal interface generator (ZIG) is
arbitrary, the LHS contributions from the patch interface may be dominant compared to
those on the cut faces for certain cut cells. The neglect of these dominant contributions can
cause severe convergence problems in some extreme cases. Since the central difference of
the inviscid ﬂux J acobian 5A does not contribute to the diagonal matrix, only the Jacobian
contributions of the smallest face due to the viscous terms and the smoothing terms need
to consider.

Refer to the Figure 30 again. Consider cell L and the smallest face b-c. The three com-

113

ponents ng, nn and "C of face normal 3 in the curvilinear coordinates C, 11 and C are cal-

culated, i.e ii = (n g, "71’ "§) . Then by analogy to equations (6.23) and (6.24) for viscosity

contribution and smoothing terms for a normal cell, and noticing that grid metrics in equa-

tions (5.19) to (5.22), the contributions of the ﬂux J acobians of smallest face b-c due to the

viscous term and smoothing term to cell L in the direction C are calculated respectively as

 

 

the following:
Arm +11 ) $2
(vrs§)bc = L 27 L bc-ngl (6.25)
PLVL
Atun+vn +wn +c S
(SM0§)bc = (I L x L y L 2' L) bc'ngl (626)

2VL

where c is the speed of sound, V is the cell volume, S b c is the area of the smallest face b-c.
The subscript L means that the values of the variable are evaluated at the center of cut cell
L. (VIS§)bC and (SMOch are added respectively to V185 and SMO: in the LHS of
the diagonal approximate-factorization algorithm (6.22) in the OVERFLOW for cell L.
Note that there are no ﬁnite difference operators in the expressions (VIS§)bc and
(SMO§)bc , which means that these two expressions are only added to the diagonal
matrix in the LHS of the implicit operator. The contributions of the smallest face b-c in

other directions to cell L and the contributions to cell R are similarly evaluated and added.

114

7 VALIDATION OF THE FULLY CONSERVATIVE
CHIMERA IN OVERF LOW

The purpose of the validation cases is to verify the basic implementation of the conser-
vative Chimera in OVERFLOW including mutual grid systems, inviscid and viscous
numerical ﬂuxes through the patch interface, special treatment in evaluating velocity gra-
dients at the patch interface and implicit treatment of the patch interface. When the data
generated by the Z16 is read in by the OVERFLOW, the test of free stream preservation is

used for error checking. All cases presented here passed this test.

7.1 Transonic Flow Through a Channel with a 10% Bump

This case has been widely used in the literature to test compressible ﬂow solver in cap-
turing shock waves. It is selected in this study to validate the implementation of the con-
servative Chimera in the OVERFLOW. The width of the channel is equal to the length of
the bump and the channel length is three times as long as the length of the bump. For the
inlet Mach number 0.675, a shock wave is generated over the bump. Two-zonal over-
lapped grids as shown in the Figure 31 are used in the testing. This two dimensional case
is actually run three-dimensionally on purpose to test the inviscid implementation of con-
servative Chimera in OVERFLOW. In the OVERFLOW two—dimensional ﬂow is run
three-dimensionally with no ﬂow in z~direction and grid points in that direction are set to

be minimum required by OVERFLOW.

115

 

Figure 31 Two-zonal overlapped grid system for the transonic ﬂow over a 10%
bump

Figure 32 shows pressure contours in one section computed by OVERFLOW with the
conservative Chimera. The contours across the patch interface appear smooth and the
shock wave are captured well. OVERFLOW with the original Chimera is also run for this
case. All the conditions such as the CFL number, smoothing coefﬁcients, etc., are set to be
the same as those in the run with the conservative Chimera. Figure 33 shows the pressure
contours computed by OVEFLOW with the original Chimera. The pressure contours in
the Figure 32 and Figure 33 agree well with each other. The comparison of the pressure
distributions on the lower wall of the channel computed with the conservative Chimera
and the original Chimera is displayed in Figure 34. The pressure distributions also agree
well with each other in general, although the pressure with the conservative Chimera has
better minimum values. Figure 35 shows the convergence histories for the both runs. Due

to the appearance of the irregular cut cells near the patch boundaries and the resultant

116

smaller local time steps, the convergence rate with the conservative Chimera is a little

slower than that with the original Chimera.

 

 

 

 

 

Figure 32 Pressure contours computed using OVERFLOW with the conservative
Chimera

 

 

 

 

 

Figure 33 Pressure contours computed using OVERFLOW with the original
Chimera

117

 

-2.0 a . r . r

 

—— Conservative

 

 

 

I
it - - - - Original
l

CD

 

 

 

1.0 i - . - 1 .
-1.5 -o.5 0.5 1.5

X

 

Figure 34 Comparison of pressure coefﬁcients proﬁles on lower wall of the channel

     
  

 

V I f I

 

[ — Original, Grid1

‘ -- — - — Conservative, Grid1
Original, Grid2 4

L—- — Conservative, Grid2

 

 

 

 

 

E -8.5 r a

:9

U)

Q)

EE

5 -9.5 - «
-10.5 l' q
_115 . 1 . 1 i "

0 500 1000 150(

Iterations

Figure 35 Comparison of convergence histories in using OVERFLOW with the
conservative Chimera and the original Chimera

118

Local grid reﬁnement is used to increase the resolution of the shock wave. An extra
zone and grid associated to the zone which reﬁned the major grid are automatically gener-
ated in the ZIG based on user speciﬁcations for the location of grid reﬁnement and the
reﬁnement factors. Figure 36 shows the computational grids after the local reﬁnement.
The ZIG also computes the interface patch information required for the conservative Chi-
mera in the OVERFLOW. Figure 37 shows the pressure contours in one section obtained
by OVERFLOW with the local reﬁnement. Compared with the pressure contours in Figure
32, the contours shown in Figure 37 display an improved shock wave in the reﬁnement

ZOHC.

 

Figure 36 Grids after local reﬁnement for the simulation of ﬂow over a bump by
OVERFLOW with the conservative Chimera

 

 

 

 

 

 

 

 

Figure 37 Pressure contours computed by OVERFLOW with the conservative
Chimera and with local grid reﬁnement

119

7.2 Laminar and 'Ihrbulent Boundary Layers over a Flat Plane

The ﬂat plane boundary layer problem is used to validate the implementation of vis-
cous ﬂuxes on the patch interface of the conservative Chimera in OVERFLOW. Figure 38
shows the two—grid system used in the testing, namely the major grid and the minor grid.
The major grid is the one ﬁlling the whole ﬂow ﬁeld over the ﬂat plane, which is plot in
black. The minor grid is the one contained in the major grid shown in red lines in the ﬁg-
ure. The lower boundary of the minor grid is constructed such that it extends well into the
boundary layer of the ﬂat plate. Hence the patch interface will be within the boundary
layer. Since the velocity changes rapidly inside the boundary layer, the correct evaluation
of velocity gradients, e.g 3%, is very important. Therefore this case is a good test for the
method illustrated in the sections 6.2 and 6.3. Our objective is to observe the smooth vari-

ation of streamwise velocity u in boundary layer across the patch interfaces.

 

Figure 38 Grids for the study of boundary layer growth over a ﬂat plane

120

Free stream Mach number is 0.2 and free stream Reynolds number is 106. OVER-
FLOW with conservative Chimera is run for both laminar ﬂow and turbulent ﬂow with
Baldwin—Lomax turbulence model. The boundary conditions for the ﬁrst several levels of
the grid points on the plane from the inlet are set to inviscid adiabatic walls; Figure 39 and
Figure 40 show the velocity component in the free stream direction near the region where
the boundary layer enter the minor grid for turbulent and the laminar ﬂows respectively.
There are gaps between the major and the minor grids shown in the ﬁgures because cut
cells in the major grid are not drawn in the post processing. As it can be seen, the variation
is smooth and the solution is identical to the single grid solution before the patch interface.
Figure 41 and Figure 42 shows similar contours near the region where the boundary layer
crosses the patch interface for the turbulent and laminar ﬂows, from the minor to the major

grid this time. Again, the solution variation across the interface is predicted correctly.

 

Figure 39 Velocity distributions cross interfaces between major and minor grids
inside the turbulent boundary layer over a ﬂat plane

121

 

Figure 40 Velocity contours cross interfaces between major and minor grids inside
the laminar boundary layer over a ﬂat plane

 

 

Figure 41 Velocity contours cross interfaces between major and minor grids inside
the turbulent boundary layer over a ﬂat plane

 

 

 

Figure 42 Velocity contours cross interfaces between major and minor grids inside
the laminar boundary layer over a ﬂat plane

122

7 .3 Grid Reﬁnement Study with Thrbulent Flow over a Two-Element
Airfoil

It is well known that non-conservative numerical discretizations will introduce errors,
especially in the presence of discontinuity. Wang (1995) compared the computed solutions
using conservative and non-conservative Chimera and clearly demonstrated the superiority
of the conservative scheme in the presence of discontinuities. Here the case of subsonic
turbulent ﬂow over a two-element airfoil is used to demonstrate that the conservative inter-
face scheme also has advantages over the non-conservative counterpart even without the
presence of a discontinuity. In this demonstration case three progressively reﬁned grids,
namely coarse grid, medium grid and ﬁne grid were used to study the variation of loading
on the surfaces of the airfoils. The two—dimensional grids around the main airfoil and the
rear airfoil were generated separately. To use the OVERFLOW for the simulation of two-
dimensional problems, three planes in the spanwise direction are needed for the original
Chimera, which is constructed by stacking two-dimensional grids in the spanwise direc-
tion. This three-dimensional grid is also the ﬁnite volume grid for the conservative Chi-
mera. Then the enhanced ﬁnite different grid which has one more plane in the spanwise
direction are needed and generated by a preprocessor for the conservative Chimera. The
coarse grid used in this study is shown in Figure 43. The Chimera holes for the conserva-
tive and the original Chimera are displayed in Figure 44 and Figure 45 respectively. It is
seen the interfaces are generated such that the shear layer of the main airfoil stays within
the same grid. This is done in the ZIG by specifying the ﬁrst grid line around the airfoil as
a solid wall instead of C-cut interface. The the ZIG automatically backs off the patch inter-

face from the solid wall.

123

 

Figure 43 The coarse grid for the turbulent ﬂow over the two-element airfoil

..u l ’
.- 'mlu11um1ummummmnmu
==E- llﬂlllllllll'll
=3."

l r’ 1 I, I
H/ / /
Ill

 

"llIMMMMIIIIIIII / // / f
lllllllll 111mm H
mm

 

I

/ /
IlImlmu1mmmuumumlm1H / / / l / /
mmmm mnmmmmm I l

mum"

 

1m
11;:

 

111
11,1:
unit

11

ill"

 

 

 

 

 

 

O

‘ ~ ~:
‘ “““. ev....~.- .
‘ ‘ .. O I O. 0. -
‘wxgegézizzrr:

.

o ‘1 '
.zii-‘i
.
..i.

 

Figure 44 Patch interface for the conservative Chimera in the coarse grid

124

 

Figure 45 Holes cut for the original Chimera in the coarse grid

Flow conditions are set as the following: free stream Mach number 0.4, angle of attack
2.260 and free stream Reynolds number based on free stream velocity and main airfoil
chord length 0.9x107 . Baldwin-Lomax boundary layer model is used in the direction nor-
mal to the surfaces of the airfoils and the Baldwin-Lomax shear layer model is used in the
direction normal to the C-cuts. In the other two directions, inviscid ﬂow is assumed.
OVERFLOW with both the conservative Chimera and the original Chimera are run for this
case. All the conditions, such as CFL number, time steps, numerical methods and smooth-
ing coefﬁcients etc., for both runs of OVERFLOW with the conservative and the original
Chimera are set to be exactly same. The CPU time for each time step with the conservative
Chimera is about twice as much as that with the original Chimera. This is because four

planes in the spanwise direction are supplied for the conservative Chimera, while only

125

three planes are needed for the original Chimera. In two-dimensional simulation, ﬂows on
the ﬁrst and the last planes in spanwise direction are not computed. Therefore for the orig-
inal Chimera in two-dimensional simulation, only ﬂows on one plane are actually com-
puted. For the conservative Chimera, ﬂows on two planes are computed, and additional
calculations and memory units are required for the ﬂuxes on the patch interfaces.

Figure 46 and Figure 47 show the Mach contours around the two-element airfoil simu-
lated by OVEFFLOW with the conservative Chimera and the original Chimera in the ﬁne
grid respectively. Figure 48 and Figure 49 show the pressure contours with the conserva-
tive and the original Chimera in the ﬁne grid. The ﬁne grid around the main airfoil is of
dimension 482 x 162 and the ﬁne grid for the rear airfoil has 242 x 62 grid points, which
is about twice in each direction as ﬁne as the coarse grid shown in Figure 43. On this ﬁne
grid Mach and pressure contours obtained with different Chimera treatment are shown to
agree well with each other. The contour lines are seen to smoothly cross the patch inter-
faces. The boundary layers growing along the airfoils and the interaction of the shear lay-
ers in the wake of the front airfoil and the ﬂow along the rear airfoil are clearly shown.
Figure 50 shows the convergence histories. It is shown the convergence rates with differ-
ent Chimera approaches are about same. Figure 51 shows the comparison of pressure dis-
tributions along the two-element airfoil obtained with the conservative and original
Chimera on three different grids. Figure 52 shows the comparison on just the rear airfoil. It
is shown that the Cp proﬁle obtained with the conservative Chimera on the ﬁne grid agrees
well with that obtained with original Chimera on the same grid. On the coarser grids
(coarse grid and medium grid) the conservative Chimera, however, yielded signiﬁcantly

better results than the original Chimera.

126

 

 

 

Figure 46 Mach contours around twoeelement airfoil computed by OVERFLOW
with the conservative Chimera

 

 

 

 

Figure 47 Mach contours around two-element airfoil computed by OVERFLOW
with the original Chimera

127

 

Figure 48 Pressure contours around two-element airfoil computed by OVERFLOW
with the conservative Chimera

 

 

 

Figure 49 Pressure contours around two-element airfoil computed by OVERFLOW
with the original Chimera

128

 

 

 

 

 

 

 

 

'8.0 _ T l r
if «1. — ~ - - Conservative, Grid1
’ "wt, -— Original, Grid1
‘3‘. — -— - Conservative. Grid2
‘H - ., - Original Grid2
. I 9
-9 0 )— AX —t
To " ,2x
3 -\ _
% «Ry: “~
o i ..Y ( x
e; x I
8 \K‘ 15‘5””. \
.1 ‘ ‘ , s...” .
‘10-0 *— M- \ g‘ ‘ \“' a. \Iur‘:"—*-'—ul:.;.:‘:u\ d
\\\\ ‘L\‘\. I mq\
...... \d ,
2 ‘ ‘s l
l- «‘r?" '01: ‘1‘“..MW‘
“Wad“; . 3,. f
-11.0 ' 1 - 1 . r . "illllrl'
0 400 800 1200 1600
Iteration

Figure 50 Convergence histories for OVERFLOW with conservative Chimera and
original Chimera in the ﬁne grids

 

  
  

 

   

 

 

 

 

 

'1.0 ' l V l v 1
R — -— Conservative, Coarse
FA) --- Original, Coarse
. Mil) ------ Conservative, Medium 1
,' ‘13., «~— Original. Medium
l ‘- - — Conservative, Fine
'0-5 * l Original, Fine
l
Q l
O
0.0 e
§
l-
o.5 - 1 - 1 . 1 .
-0.25 0.25 0.75 1.25 1.75
X/C

Figure 51 Pressure distributions along the airfoil computed by OVERFLOW with
the original and the conservative Chimera in different levels of grids

129

Cp

0.0

--— Original, Coarse
0.2 --— Conservative,

---- Original, Medium

— Conservative, Fine
0,4 —— Fine

 

1.15 1.25 1.35 1.45 1.55 1.65 1.75
X/C

Figure 52 Pressure distributions along the rear airfoil computed by OVERFLOW
with the original and the conservative Chimera in different levels of grids

7.4 Turbulent Flow over a Three-Element Airfoil

The case of ﬂow over three-element airfoil has been under extensive studied both
numerically (Soetrisno et al., 1994) and experimentally (Valarezo et al., 1991). Three
spanwise planes for the original Chimera and four spanwise planes for the enhanced ﬁnite
difference grids for the conservative Chimera need to be supplied to OVERFLOW. Four
different grids, namely coarse grid, medium grid. ﬁne grid and the ﬁnest grid, were used
for the study of the ﬂow over the three-element airfoil with OVERFLOW. The outer

boundary of the main element grid is about ten chord length away from the airfoil. Figure

130

53 shows part of the medium grid around the airfoil and the patch interfaces and holes
generated in the ZIG. Note that the interface between the main element grid and the front
element grid backs off to avoid the solid wall of the main element. The ﬁnest grids for the
main, front and rear elements have dimensions 475 by 160, 481 by 77 and 457 by 79
respectively, which are about twice as ﬁne in each grid direction as those shown in Figure
53. The ﬁne grids are of dimension 360 by 65, 269 by 57 and 297 by 49 for the three ele-
ments respectively. Flow conditions were set as the following: free stream Mach number
M“ = 0.2, free stream angle of attack a: 8.20 and Reynolds number based on free
stream velocity and main airfoil chord length Re: 0.9x107. The Baldwin-Lomax turbu-
lence model is used in the same way as in the two-element airfoil case. The run is stopped

after the pressure coefﬁcients around airfoil do not change for 500 time steps.

1111""
m
I'llllmum

   

      
  
     
  

      
   
  
  

\ \\\\ "'f l:
\ “ triennial {jgjzlllli'l'l'l'll'..'....

      
   
        
 
 

a)

%
//
//
Z/i

   
   
     
  
 
 

(a,

  
  
    
   

 

42/

  
   

    

  

//
%///
,/
é?)

/
// . g3"

  

//
//
//

//
/////
////
/////

 
 

\_

4%
///

 
  
     

   

 

{if}:
. . .»/II II/// .

”I ’III/// ,
0W

”///////
/////

   
 
 
 
 
 
   

 

\l
l‘l‘"
§g¢ ._
ii

     
 
   
  
 
 

_

  
   
   

l
1

  
 

I
ll 1
‘1

111
\\\\“(11\

  
    
  
  
    
 
  
  
  
   
  
  
 

\\

g;

   
 
  

\\

  
 
  
 

  
 
 
 
  
 

1 ll

  
 
 
 

 
 
 
 

 
   

Figure 53 Medium grids and patch interfaces for the thee-element airfoil

131

 

 

 

Figure 54 Pressure contours around the three-element airfoil computed using
OVERFLOW with the conservative Chimera on the ﬁnest grid

 

 

 

Figure 55 Pressure contours around the three-element airfoil computed using
OVERFLOW with the original Chimera on the ﬁnest grid

132

Figure 54 and Figure 55 show the pressure contours around the airfoil simulated with
OVERFLOW on the ﬁnest grid. It is observed that the contour patterns shown in Figure 54
and Figure 55 obtained with different methods are similar. Figure 56 and Figure 57 show
the distribution of the pressure coefﬁcients along the airfoil surface computed with the dif—
ferent grids by different methods compared to experimental data. It is shown that the
results are signiﬁcantly improved when the grid resolution increases from the coarse grid
to the medium grid. The results for the ﬁnest grid and the ﬁne grid are, however, almost
same. This shows that the ﬁnest grid is ﬁne enough for this demonstration case, at least for
the methods and turbulent model used in this case. The comparison of the pressure distri-
butions along the airfoil surfaces computed with the conservative and the original Chimera
is displayed in Figure 58. Note that the peak values in the pressure proﬁles computed with
the original Chimera are better than those with the conservative Chimera. On the other
hand the conservative Chimera yields better predictions near the stagnation points, espe-
cially for the front element. In general, the Conservative Chimera gave similar solutions as
the Original Chimera. OVERFLOW with the conservative Chimera and the original Chi-
mera is also run for the inviscid ﬂow over this three—element airfoil with same free Mach
number on the ﬁne grid. Figure 59 shows the comparison of pressure coefﬁcients along the
airfoils for the inviscid ﬂow with the experimental data. It is shown that the pressure pro-
ﬁles obtained with both Chimera treatments are almost same. Comparing Figure 59 with
Figure 56 and Figure 57, it is obvious that the results of the inviscid ﬂow agree better with
experimental data. The discrepancy between the experimental data and the numerical
solutions on the ﬁnest grid may imply that the Baldwin-Lomax turbulent model is not

good enough for turbulent ﬂows over the three-element airfoil.

133

'6.0 ' T ' i

 

 

 

 

 

 

 

 

Q3
, 1% 0 Experiment
‘31; D - - - - Coarse Grid
an, ———- Medium Grid
4-0 " 0 ~ - Fine Grid —
. ‘1‘"? — —- Finest Grid
1 l \b j
1 <10 3,
\\D
g :1
8 -2 0 ~53 ~ 5.61.. 3% -
0 d3 0 \ ‘\ D a 1! \
53‘1“. \ ‘ ~E‘Ef;\&293 E ‘13
l l I Q5121
', . l l V
0.0 atria; '1 D .
ii” 1 2
U @M 1
2.0 1 1 - 1 .
0.0 0.5 1.0 1.5
X/C

Figure 56 Pressure distributions along airfoil surfaces computed using OVERFLOW
with the original Chimera

 

 

 

 

 

 

-6.0 . r . ,
C’é
‘3 :33 0 Experiment
011 D -—-- Coarse Grid
llkgnb -——-— Medium Grid
-4.0 - CR1 — - - Fine Grid 3
' \‘D —«— Finest Grid
V C]
O. &l Xi? 1:1 1%
0 '2 0 F “660 I N l, d

0.0 F319;

 

 

2.0 . 1 . 1 1
0.0 0.5 1.0 1.5

X/C

 

Figure 57 Pressure distributions along airfoil surfaces computed using OVERFLOW
with the conservative Chimera

134

 

 

 
   
   

 

 

 

 

 

so . . . . l
Eh . I
r3 U Experiment ’1
. --— Conservative, Coarsel l
l ---— Original. Coarse ‘
4.0 _ 1 ‘ — Conservative, Finest
\ :— Original, Finest
3 -2.0 ﬁlm.
,,
0.0 — a
2.0 1 1 . 1 1 4
0.01 0.26 0.51 0.76 1.01 1.26

X/C

Figure 58 Comparison of pressure distributions along the airfoil surfaces computed
using OVERFLOW with the conservative and the original Chimera

 

 

 

 

 

 

 

 

'6-0 ' 1 l
l rt:
1 D a Experiment
1 E1"; ~ — - , Conservative
l l C * Original
-4.o - 11 Eb ~
1 i?
‘ K $1
0. _2 0 ‘Cl l D‘ _
O l Ell-5E; ‘41.; g 1:; %3\\
l ‘ D
=9 “ E1:810;
0'0 I??? l “in _
I WM
2.0 L 4* 1 1 1
0.0 0.5 1.0 1.5

X/C

Figure 59 Comparison of pressure distributions for the inviscid ﬂow along the airfoil
surfaces computed with OVERFLOW.

135

7 .5 Summary

Tests with inviscid, laminar and turbulent ﬂows with different geometries have vali-
dated the implementation of the fully conservative Chimera in OVERFLOW. The results
obtained with the conservative Chimera and the original Chimera approaches for most
cases are similar. For the turbulent ﬂow over the two-element airfoil the pressure distribu-
tions along the rear element obtained with the two Chimera approaches in the ﬁne grid are
almost identical. If the results in the ﬁne grid is assumed to be close enough to the actual
solution, the pressure distributions along the rear element obtained with the conservative
Chimera approach in the coarser grids (coarse and medium grids) are much closer to the
actual solution than those with the original Chimera approach. But for the turbulent ﬂow
over then three-element airfoil the peak values in the pressure proﬁles computed with the
original Chimera are better than those with the conservative Chimera. Therefore for the
smooth ﬂow we may conclude the following: In general the conservative Chimera and the
original Chimera yield similar results. The conservative Chimera may make a difference
as compared to the original Chimera for some cases in relatively coarse grids.

In the testing of the conservative Chimera with the turbulent ﬂow over the three-ele-
ment airfoil we found that the computation can not process without the special treatment
in evaluating velocity gradients at the patch interfaces. The implicit treatment of the patch
interface are found to be necessary to boost the stability and accelerate the convergence
rates for all the validation cases. Due to the appearance of the irregular cut cells near the
patch boundaries, without the implicit treatment of the patch interface the convergence
rate with the conservative Chimera approach is slow and not competitive to that with the

original Chimera approach.

136

8 DEMONSTRATION OF THE FULLY CONSERVATIVE
CHIMERA IN OVERFLOW

The purpose of the demonstration cases is to show the capability of the conservative
Chimera in tackling more complicate ﬂow problems. Since there are not experimental data
available for the demonstration cases, the results obtained with the conservative Chimera

can only be qualitively evaluated based on the visual images of surface and contour plots.

8.1 Hypersonic Flow Around a Launch Rocket

This case contains seven grid components: one main rocket and six smaller boosters.
Only half of the ﬂow ﬁeld was calculated because of symmetry. The four zone mesh in the
half ﬁeld is shown in Figure 60. For the main rocket 81 grids is distributed in the stream
wise direction, 39 grids in the direction normal to the body surface and 61 in the circum-
ferential direction. Each booster’s grid is of dimension 31 x 9 x 21 . The outer boundaries
of the boosters’ grids are used as the patch interfaces. The free stream Mach number is 6
and the angle of attack is set to 2.50. Figure 61 shows the Mach contours on the surfaces
of the main rocket and the boosters as well as the Mach contours on the last across stream-
wise plane in the main grid domains. Figure 62 displays a close-up view of Figure 61 on
the region around the boosters. The effects of the presence of the boosters can be clearly
seen on the surface of the rocket. These ﬁgures also show the intricate shock wave patterns
resulted as the ﬂow moves over the conﬁguration. In particular the shock patterns from the

smaller boosters, and their reﬂection off the surface of the main booster are well captured.

137

Q
- \‘~
-.1.:..

~ - .‘
.~ .*..~

 

Figure 60 Four zone Chimera grid for hypersonic ﬂow around a launch rocket

 

Figure 61 Mach distributions on the surfaces of the rocket and the boosters and on
the last cross-streamwise plane computed using OVERFLOW with the conservative

Chimera
138

 

Figure 62 Mach contours on the surfaces of the rocket and the boosters

8.2 Demonstration of Local Reﬁnement with Hypersonic Flow Around a
Launch Rocket

The case of hypersonic inviscid ﬂow over a rocket with six boosters has been shown in
the last subsection. This case is used again to demonstrate the capability of local reﬁne-
ment implemented in the ZIG. The grids without local reﬁnement is shown in Figure 60. A
region in the main grid, which encloses the grid of the middle booster as shown in the Fig-
ure 63, is chosen to be reﬁned. The factors of the reﬁnement in all three grid directions are
set to 2. Then the grid in this region is automatically generated in the ZIG by locally reﬁn-
ing the main rocket grid by the given factor in each direction. Let the grid of the rocket and
grids of three boosters be zone 1, 2, 3 and 4 respectively, and the grid in the reﬁned region

be zone 5. Then there are patch interfaces between zone 1 and zone 5, interfaces between

139

zone 1 and zone 2, interfaces between zone 1 and zone 4, and interfaces between zone 5
and zone 3. The grids of zone 2 and zone 4 share the major grid of zone 1. The grid in zone
5 is the minor grid for the interfaces between zone 1 and zone 5, as well as the major grid

for the interfaces between zone 3 and zone 5.

.~:.

.. - <
o.‘.~:.~ 1
..... .9

i
W?»
I, I.

l I
Hi

II“

’n
1

s
-
y
:6
.o

 

Figure 63 Grids around the rocket and the boosters with local reﬁnement

The ﬂow conditions are set to the same as those in the case without local reﬁnement.
Both OVERFLOW and CFD-FASTRAN from CFDRC with conservative Chimera are run
for this case. Figure 64 shows the Mach contours on the surfaces of the main rocket and
the boosters, as well as the Mach contours on the last cross-streamwise plane of the main
grid. Figure 65 shows the density distributions on the surfaces in the region around the

boosters. These results are simulated by OVERFLOW with the conservative Chimera. It

140

can be seen that the shock waves inside the local reﬁnement zone are signiﬁcantly sharper
than those outside. The reﬂection of the shock waves of boosters off the surface of the
main rocket are well captured. The ﬁgures also show that ﬂow variables are smooth across
the patch interfaces and intricate shock wave patterns that result as the ﬂow moves over
the boosters are resolved well. Comparison of the results with local reﬁnement shown in
Figure 64 and Figure 65 and the results without local reﬁnement shown in the Figure 61
and Figure 62 conﬁrms the improvement with local reﬁnement. The Mach contours on the
surfaces of the rocket and a cross-section plane simulated with CFD-FASTRAN are
shown in the Figure 66 and Figure 67. The improvements of the contour resolutions by the

local reﬁnement are also clearly seen in these ﬁgures.

 

/

Figure 64 Mach contours on the body surfaces and on the last cross stream plane
simulated with OVERFLOW with the conservative Chimera

 

141

 

Figure 65 Density contours on the body surfaces computed with OVERFLOW with
the conservative Chimera

 

Figure 66 Mach distributions on the body surfaces and some sections simulated with
CFD-FASTRAN with local grid reﬁnement

142

 

Figure 67 Mach contours on the body surfaces simulated with CFD-FASTRAN with
local grid reﬁnement

8.3 Flow over a Wing and Store Combination

The geometry is a combination of a symmetric wing and a store under the wing. The
physical domain around the store is chosen to be big enough so that the outer boundary of
the store grid intersects the wing. This is done on purpose to test the automatic interface
backoff algorithm implemented in the ZIG. Figure 68 shows the grids on the surfaces of
the wing and the store, and grids on two across spanwise sections. Figure 69 shows the
patch interface, the smallest faces and the holes cut automatically by the ZIG with differ-
ent viewing angles. It is seen that the smallest faces have arbitrary shapes, which are
shown in the different colors in Figure 69. The patch interfaces consist of all the smallest

faces, which are pushed back from the wing to avoid the wing as shown in the ﬁgure.

143

    
  
 

    
   
 

        

I
\\ \
1‘ \\\\“ \‘

\\\\\\\

Figure 68 Grids and the conﬁguration of the wing-store combination

       

  

\

 

 

Figure 69 Patch interfaces, smallest faces and holes cut by the ZIG

144

The free stream Mach number is 0.5 and the angle of attack is set to zero. A symmetric
Boundary condition was applied to the end of the wing close to the store. The other end of
the wing is extended far away from the store, simple extrapolation was applied to the
boundary at that end. OVERFLOW with the conservative Chimera was used to simulate
the inviscid steady ﬂow over the geometry. Figure 70 and Figure 71 show the pressure and
Mach contours on the surfaces of the wing and the store as well as the pressure contours
on one across spanwise section respectively. The effects of the presence of the store can be
seen from the contours on that section. Contours across the patch interfaces are smooth.
Figure 72 shows the convergence histories for both grids. This case once again validates
the inviscid implementation of conservative Chimera in the OVEFFLOW and demonstra-

tion the ability of the Chimera handling complex geometries.

 

Figure 70 Pressure contours on the surfaces of the wing and the store and on one
spanwise section computed using OVERFLOW with the conservative Chimera

145

 

Figure 71 Mach contours on the surfaces of the wing and the store and on a
spanwise section computed using OVERFLOW with the conservative Chimera

.__J__

— Wing Grid
— — — — Store Grid

 

 

 

1
73 -7.0 — i
:9
U}
a

-9.0 — i

1
-11.0 . 1
o 500 1000

1 500
Iteration

Figure 72 Residual history for the simulation of ﬂow over wing and store
combination by OVEFLOW with conservative chimera

146

8.4 Transonic Flow Around a Combination of Wing and Missiles

The geometry in this demonstration case is a combination of an extended delta wing
and two missiles under the wing. There are four ﬁns on the surface of each missile. Figure
73 shows the geometry including the shadows of the ﬁns. Figure 74 shows surface grids,
one section of grids around the wing and the outer boundaries of the two missile grids.
Note that the physical domains around the missiles are set to be big enough to intersect the
wing. The grid around the wing has 137 points in the streamwise direction, 39 points in
the spanwise direction and 31 points in the cross streamwise direction. The grid of each
missile is of the dimension 50 x 19 x 69 , with 50 points in the streamwise direction, 19
points in the direction normal to the body surface and 69 points in the circumferential
direction. Figure 75 shows the outer boundaries of the missiles after the cutting in the ZIG,
which are parts of the patch interfaces between the missile grids and the grid around the
wing. It is seen the patch interfaces are backed off from the surface of the wing to avoid

the wing.

.4. . 1.5rn11n-EIFVW 2"..7‘ “‘4’2'5-2‘147‘.“

  

 

  

Figure 73 Conﬁguration of the combination of wing and two missiles

147

. \\\“‘“
» \\\\\“
_ \\\\‘\
. 1 ’\\\\\
‘ \\\

ill,
1

'1'
Mill”

1"”l

I
l

   

I3:11,?" n ll
Will

ll

- u

, . ;\--.=_=‘.=_
g ... . . _ \\\\\\

\\
\\\\\\\
\\

 

      

.\\
\

          
 

    

Figure 74 Grids for the combination of wing and two missiles

 

Figure 75 Holes automatically cut in the ZIG

148

The free stream Mach number is 0.8 and the angle of attack is set to zero. Both CFD-
FASTRAN and OVERFLOW with conservative Chimera are run to simulate the ﬂow over
the combination of the delta wing and two missiles. The boundary condition at both ends
of the delta wing set to be symmetric. The ﬁns on each missile are set to be blockages for
CFD-FASTRAN, and CFD-FASTRAN can automatically treat the boundaries of the ﬁns
as walls. In OVERFLOW, however, solid wall boundary conditions must be speciﬁed in
the input ﬁle for the ﬁve surfaces of each ﬁn on the missiles. For each missile, twenty
boundary conditions for four ﬁns, and ﬁve more boundary conditions for the missile grid
are needed to set up the problem. Figure 76 and Figure 77 show pressure contours on two
cross sections and the pressure distributions on the surfaces simulated by OVERFLOW
and CFD-FASTRAN with the conservative Chimera respectively. Figure 78 and Figure 79
display the pressure distributions on the surfaces and two cutting planes simulated by
OVERFLOW and CFD-FASTRAN with the conservative Chimera respectively. The
effects of the missiles ﬁns on the ﬂows are clearly seen from these ﬁgures. It is observed
that the solutions computed with OVERFLOW and CFD-FASTRAN are quite similar to

each other as they are expected.

149

 

Figure 76 Pressure contours on the surfaces and two streamwise sections simulated
using OVEFLOW with the conservative Chimera

.4" 1-

k

' :1'4

4,
1i
ii

O

 

Figure 77 Pressure contours on the surfaces and two streamwise sections simulated
using CFD-FASTRAN with the conservative Chimera

150

 

Figure 78 Pressure contours on the surfaces and two spanwise sections simulated
using OVEFLOW with the conservative Chimera

 

Figure 79 Pressure contours on the surfaces and two spanwise sections simulated
using CFD-FASTRAN with the conservative Chimera

151

9 CONCLUSIONS

9.1 Summary

A computational methodology has been developed in the ﬁrst part of the thesis for the
simulations of acoustic radiation, prOpagation and reﬂection. The developed methodology
has high order of accuracy, uses less grid points per wave length comparing to standard
high order numerical methods, and automatically damps out spurious short waves. Fur-
thermore, the methodology can be applied to acoustic problems in the presence of objects
with curved geometries. To achieve these results, high order accurate optimized upwind
schemes, which are applied to discretize spatial derivatives on interior grid points, have
been developed. High order accurate optimized one-side biased schemes, which are only
applied to discretize the spatial derivatives on grid points near computational boundaries,
have also been constructed. The developed schemes are combined with a time difference
scheme to fully discretize acoustic ﬁeld equations in multi-dimension in arbitrary curvilin-
ear coordinates. Numerical boundary conditions were investigated and intuitively illus-
trated. Applications of the developed methodology to a sequence of one-dimensional and
multi-dimensional acoustic problems are performed. The numerical results have validated
the developed methodology, demonstrated advantages of the methodology, and showed
that the objectives of the methodology has been achieved. More speciﬁcally, the conclu-
sions and work done in the ﬁrst part of the thesis are summarized as the following:

A method of constructing a class of high order, optimized, upwind, ﬁnite difference

schemes for discretizing spatial derivatives based on the analysis on the model wave equa-

152

tion has been developed. An fourth order accurate, 7-stencil, optimized, upwind scheme
has been actually constructed based on the developed method. The comparison of the
effective numerical wavenumbers of different schemes shows that the 7—stencil optimized
upwind scheme has built in dissipation so that it can automatically damps out spurious
short wave. It has much less dissipation error than those of the standard six-order accurate
upwind scheme and backward spatial schemes developed in Tam &Webb’s DRP method.
It is able to resolve the waves with as high wavenumbers as those resolved by Tam
&Webb’s central DRP. The tests of a initial Gaussian pulse and a initial step wave on the
model wave equation have validated the advantages of the optimized upwind scheme.

Acoustic ﬁeld equations in different forms and in generalized curvilinear coordinate
systems are derived for integrity. The Jacobian matrices and their corresponding eigenval-
ues and eigenvectors associated to the ﬁeld equations are also derived. It is illustrated that
disturbance waves can be decomposed as entropy waves, vorticity waves and acoustic
waves. Acoustic waves are treated differently in the matrix splitting techniques according
to their propagating directions. Matrices similar to arbitrary diagonal matrices through
similarity transformations of eigenvectors matrices are derived. The correct forms of these
matrices are important in the matrix splitting techniques and convenient for the future
applications. The numerical tests of the multi-dimensional problems have validated the
correct derivations of the matrices.

The optimized upwind schemes have been extended to acoustic ﬁeld equations in
multi-dimension in generalized curvilinear coordinates with the matrix splitting tech-
niques. High order one-side biased optimized schemes used for grid points near the com-

putational boundaries have also been developed. The boundary conditions for curved wall

153

are developed. Radiation and outﬂow boundary conditions in the far ﬁeld are intuitively
illustrated. The boundary conditions are implemented with one-side biased optimized
schemes. The optimized upwind scheme together with the one-side biased optimized
schemes and the DRP temporal difference scheme gives an numerical methodology for the
computational aeroacoustics. With the extension of the optimized upwind schemes to
multi-dimension and the derived boundary conditions, the methodology can then be
applied to accurately solve the acoustic problems in the presence of objects with curved
geometrizes.

A sequence of direct numerical simulations on one-dimensional and multi-dimen-
sional acoustic problems have been carried out to test the effectiveness of the developed
methodology for CAA. Differences between numerical results obtained with the opti-
mized upwind scheme, central DRP methods, standard six-order accurate upwind
schemes, and analytical results are compared for many cases. In comparisons, the best
effort was made to ensure that all other factors remains the same. These factors include the
computational grid, the time step, etc. The tests show that the developed optimized upwind
scheme automatically damps out the spurious short waves, without the need of ﬁlters or
explicit dissipation terms. The tests show the optimized upwind scheme provides more
accurate solutions than the central DRP methods for all the simulations with the exact
solutions available. Central DRP methods without adding damping term cause spurious
oscillations in some cases. Adding dissipation can eliminate the seriously spurious oscilla-
tions of numerical results. The test of the acoustic radiation from a oscillating piston
clearly demonstrates that computational error associated with the optimized upwind

scheme is signiﬁcantly less than that associated with the central DRP with ﬁnely tuned

154

 

damping term. The optimized upwind scheme minimizes not only the disserpation errors
but also the dissipation errors, while retaining the numerical stability. The developed
methodology have been applied to the scatterings of acoustic pulses off a cylinder and a
sphere. The curvilinear transforms are employed to test the extension and implementation
of the methodology for acoustic ﬁeld equations in multi-dimension in generalized curvi-
linear coordinates. Although there are no direct comparisons available for these simula-
tions, the obtained wave patterns have validated the extension and implementation of the
methodology. The numerical results demonstrated that the developed methodology has the
ability to handle the curved walls and can be applied to acoustic problems when more
complex geometries involved. This will be discussed more in detail in the suggestions for
future work. The objective of investigating, developing and applying optimized upwind
schemes for CAA has been achieved in the ﬁrst part of the thesis.

The second part of the thesis deals with a fully conservative Chimera methodology.
The fully conservative Chimera was developed originally based on ﬁnite volume
approach. A ﬁnite difference scheme, however, is shown to be identical to a ﬁnite volume
scheme with proper deﬁnition of control volumes and matrices. The fully conservative
Chimera has been successfully extended to ﬁnite difference numerical scheme for viscous
ﬂows including turbulent models and successfully implemented into NASA’s widely used
code OVERFLOW. In the implementation the Roe’s numerical ﬂuxes are employed to the
inviscid ﬂuxes across the patch interfaces. The thin layer approximation is used in the cal-
culating of viscous ﬂuxes in the patch interfaces. The implementation also includes
implicit treatment of the patch interfaces.

Tests on several inviscid, viscous laminar and turbulent, exterior ﬂows over objects

155

 

with complex geometries were performed to validate the implementations and demon-
strate the conservative Chimera methodology. Numerical tests indicate that implicit treat-
ment of the patch interface is extremely important to the stability of the overall time-
marching procedure. Due to the appearance of very irregular cut cells with potentially
diminishing volumes and time steps, the time-marching procedure is found to be quite
unstable without the implicit treatment. With the implicit treatment the stability of the
numerical scheme with the conservative Chimera is Competitive to that of the original Chi-
mera.

The implementation of the fully conservative Chimera into the OVERFLOW allows
the errors in the traditional Chimera approach due to non-conservative data interpolations
to be evaluated through direct comparisons. A grid reﬁnement study with a subsonic two-
element airfoil ﬂow case shows that the conservative Chimera makes a difference as com-
pared to the original Chimera even for ﬂows without discontinuities. The difference, how-

ever, diminishes with global grid reﬁnement for smooth ﬂows.

9.2 Suggestions to Future Work

For the numerical consideration, the time marching scheme which may be more suit-
able than the DRP time difference scheme for the developed optimized upwind scheme
should be investigated and developed. The effective numerical wavenumber or frequency
of full numerical discretizations by both the spatial and temporal schemes should be con-
sidered together, because the numerical errors attributes to both schemes. Treating the spa-
tial and temporal schemes separately can not accurately reﬂect global effects of the

numerical errors caused by the approximations of the effective numerical wavenumber or

156

frequency to the actual wavenumber or frequency. The errors of numerical schemes on the
acoustic ﬁeld equations in multi-dimension should be analyzed to better understand the
effects of the matrix splitting techniques on numerical accuracy and stability. Furthermore
the numerical errors caused by the boundary conditions are still not clear and need to be
analyzed. The test case on the reﬂection of an acoustic pulse off a straight wall in a two-
dimensional uniform ﬂow clearly show the errors caused by the boundary conditions.

For the application consideration, grid generation for complicate geometries requires
the most of man power. Furthermore there is not easy way to accurate evaluated the matri-
ces of a given grid. Body-ﬁtted grid systems near arbitrary objects and uniform Cartesian
grid system in the farﬁeld with Chimera grid technology can be an effective way to handle
the complex geometries. The optimized upwind scheme combined with a ﬁnite volume
discretization approach should be investigated and developed. The optimized upwind
scheme with the ﬁnite volume discretization approach combining with the chimera tech-
nology will provide an effective way for simulating acoustic problems when the objects
with very complicate geometries are involved. The mean ﬂows can be obtained by recent
commercial CFD software. The proposed methodology can then be used to solve the
acoustic ﬁeld equations with the given mean ﬂows. For example, we may use CFDRC’s
geometrical modeling/grid generation package, CFD-GEOM for grid generation,
CFDRC’s ﬂow solver, CFD-FASTRAN for the simulation of mean ﬂows, a CAA module
with proposed methodology to predict aeroacoustic waves, and CFDRC’s post-processing
package CFD-VIEW for ﬁeld visualization. This will result in a computational environ-
ment which would present the opportunities to investigate sound generating mechanicsms,

directivity, spectral features and the propagation of the sound to the far ﬁeld.

157

REFERENCE

158

[l].

[2].

[3].

[4].

[5].

[6].

[7].

[8].

[9]-

[10].

[11].

[12].

[13].

REFERENCE

Ahmad, J., and Duque, E.P.N., “Helicopter Rotor Blade Computation in Unsteady
Flows Using Moving Embedded Grids,” AIAA 94-1922, 1994.

Anderson, D. A., Tannehill J. C. and Pletcher, R. H., “Computational Fluid
Mechanics and Heat Transfer,” Taylor & Francis, 1984.

Barth, T.J., “Higher Order Solution of the Euler Equations on Unstructured Grids
Using Quadratic Reconstruction,” AIAA Paper 90-0013, 1990.

Benek, J .A., Buning, RG. and Steger, J .L. “A 3D Chimera Grid Embedding Tech-
nique,” AIAA Paper 85-1523, 1995.

Benek, J .A., Steger, J.L. and Dougherty, F.C., “A Flexible Grid Embedding Tech-
nique with application to the Euler Equations,” AIAA Paper 83-1944, 1983.

Berger, M.J., “On Conservation at Grid Interfaces,” SIAM Journal of Numerical
Analysis, vol. 24, 1987.

Blanco, M. and Zingg D.W., “Fast Newton-Krylov Method for Unstructured
Grids,” AIAA Journal, vol. 36, No.4, 1998.

Buning, PG and Chan, W.M., “OVERFLOW/F3D User’s Manual,” NASA Ames
Research Center, Moffett Fiueld, California, February, 1991.

CFD Research Corporation, “CFD-FASTRAN User’s manual,” CFDRC Report
GR-93-2, 1993.

CFD Research Corporation, “CFD-GEOM User’s manual, Interactive Geometric
Modeling and Grid Generation Software,” Feb., 1998.

CFD Research Corporation, “CFD-VIEW User’s manual, 3-D Computer Graphics
and Animation Software,” Feb., 1998.

Chakravarthy, SR. and Osher, S., “A New Class of High accuracy TVD Schemes
and Their Applications,” Journal of Computational Physics, vol. 68, pp. 151-179,
1987.

Chen, RE, and Zhuang, M., “Application of Dispersion-Relation-Preserving
Scheme to the Computation of Acoustic Scattering in Benchmark Problems.” in
ICASE/LaRC 2nd workshop on benchmark problems in Computational Aeroacous-
tics, Tallahassee, FL., Nov., 1996,

159

[14].

[15].

[16].

[17].

[18].

[19].

[20].

[21].

[22].

[23].

[24].

[25].

[26].

[27].

Davis, D.L. and Thompson, H.D.L., “The Impact of Computational Zonal Interfac-
ing on Calculated Scramjet Performance,” AIAA Paper 92-0390, 1992.

Duque, E.P.N., Biswas, R., and Strawn, R.C., “A Solution Adaptive Structured/
Unstructured Overset Grid Flow Solver with application to Helicopter Rotor
Flows,” AIAA Paper 95-1766, 1995.

Duque, E.P.N., and Strawn, R.C., “An Overset Grid Navier-Stokes Kirchhoff-Sur-
face Method for Rotorcraft aeroacoustic Predictions,” AIAA 96-0152, Jan. 1996.

Eberhardt, S. and Baganoff, D., “Overset Grids in Compressible Flow,” AIAA
Paper 85-1542, 1985.

Eversman, W., “Theoretical Models for Duct Acoustic Propagation and Radiation,”
in Aeroacoustics of Flight Vehicles: Theory and Practice, ed. H.H. Hubbard, NASA
RP-1258, Chapter 13, 1991.

Farassat, F, “The Kirchhoff Formulas for Moving Surfaces in Aeroacoustics-The
subsonic and Supersonic Cases,” NASA Technical Memorandum 110285, 1996.

Hardin, J.C., Ristorcelli, J .R., and Tam, C.K.W., ICASE/LaRC workshop on
Benchmark Problems in Computational Aeroacoustics. NASA CP 3300, May 1995.

Harten, A., Engquist, A., Osher, S. and Chakravarthy, S., “Uniformly High-Order
Accuracy Essentially Non-Oscillatory Schemes 111,” Journal Computational Phys-
ics, vol. 71, pp. 231-323, 1987.

Hixon, R., “Evaluation of a High-Accuracy MacCormack-Type Scheme Using
Benchmark Problems,” NASA Contractor Report 202324, ICOMP-97-03, 1997.

Hixon, D.R., Shih, S., and Mankbadi, R.R., “Evaluation of Boundary Conditions
for Computational Aeroacoustics,” AIAA Journal, vol. 33, pp. 2006-2012, 1995.

Kim, C., Roe, PL. and Thomas, J .P., “Accurate Schemes for Advection and Aeroa-
coustics,” AIAA Paper 97-2091, 1997.

Lauder, BE. and Spalding, D.B., “The Numerical Computation of Tlurbulent
Flows,” Computer Methods in Applied Mechanics and Engineering, vol.3, pp. 269-
289, 1974.

Lax, P.L., “Hyperbolic Systems of Conservation Laws and the Mathematical The-
ory of Shock Waves,” SIAM, Philadelphia, 1972.

Lele, S.K., “Compact Finite Difference Schemes with Spectral-Like Resolution,”
Journal Computational Physics, vol. 103, pp. 16-42, 1992.

160

[28].

[29].

[30].

[31].

[32].

[33].

[34].

[35].

[36].

[37].

[38].

[39].

[40].

[41].

[42].

Li, Y., “Wavenumber-Extended High-Order Upwind-Biased Finite Difference
Schemes for Convective Scalar Transport,” Journal of Computational Physics, vol.
13. PP. 235-255, 1997.

Lighthill, M.J., “On Sound Generated Aerodynamically I. General Theory,” Proc.
Roy. Soc. A, vol. 211, pp. 564, 1952.

Lim, T.B., Sankar, L.N., and Hariharan, N ., “A Technique for the Prediction of pro-
peller Induced Acoustic Loads on Aircraft Structures,” AIAA Paper 93—4340, 1993.

Lin, S.Y., and Chin, Y. S., “Comparison of Higher Resolution Euler Schemes for
Aeroacoustics Computations,” AIAA Journal, vol. 33, No. 2, pp 237-245, 1995.

Lin, S.Y., and Chen, Y.F., “Numerical Study of MU SCL Schemes for Computa-
tional Aeroacoustics,” AIAA Paper 97-0023, 1997.

Liu, Y, “Fourier Analysis of Numerical Algorithms for the Maxwell Equations,”
Journal of Computational Physics, vol. 124, pp.396-416, 1996.

Lockard, D.P., Brentner, KS, and Atkins, H.L., “High-Accuracy Algorithms for
Computational Aeroacoustics,” AIAA Journal, vol. 33, No. 2, pp. 246-251, 1995.

Mavriplis, D.J., “Unstructured Mesh generation and adaptivity,” ICASE Report NO.
95-26, 1995.

Meakin, R.L., “A New Method for Establishing Intergrid Communication Among
Systems of Overset Grids,” AIAA paper 91-1586-CP, 1991.

Mecham, M. and Mckenna, J.T., “Cost, Not Size, To Drive Success of Super-
jumbo,” Aviation Week and Space Technology, Nov. 21, 1994.

Moon, Y.J., and Lion, M.S., “Conservative treatment of Boundary Interfaces for
Overlaid Multi-Level Grid Adaptions,” AIAA Paper 89-1980-CP, 1989.

Motsinger, RE. and Kraft, R.E., “Design and Performance of Duct acoustic Treat-
ment,” in Aeroacoustics of Flight Vehicles: Theory and Practice, ed. H.H. Hubbard,
NASA RP-1258, Chapter 14, Aug. 1991.

Parks, S.J., Buning, P.G., Steger, J .L. and Chan, W.M., “Collar Grids for Intersect-
ing geometric Components within the Chimera Overlapped Grid Scheme,” AIAA
Paper 91-1587-CP, 1991.

Pulliam, TH. and Steger, J .L., “Implicit Finite-Difference Simulations of Three-
Dirnensional Compressible Flow,” AIAA Journal, vol. 18, No. 2, pp159, 1980.

Roe, P.L., “Approximate Riemann Solver, Parameter Vectors and Difference

161

 

[43].

[44]-

[45].

[46].

[47].

[48].

[49].

[50].

[51].

[52].

[53].

[54].

[55].

Schemes,” Journal of Computational Physics, vol. 43, pp. 357, 1983.

Rogers, SE. and Pulliam, T.H., “Accuracy Enhancements for Overset Grids Using
a Defect Correction approach,” AIAA Paper 94-0523, 1994.

Rogers, S.B., Wiltberger, N.L., and Kwak, D., “Efﬁcient Simulation of Incompress-
ible Viscous Flow over Single and Multi-Element Airfoils,” AIAA Paper 92-0405,
1992.

Sankar, L.N., Reddy, N.N., and Hariharan, N., “A Third Order Upwind Scheme for
Aeroacoustic Applications,” AIAA Paper 93-0149, 1993.

Soetrisno, M., Imlay, S.T. and Roberts, D.W., “A Zonal Implicit Procedure for
Hybrid Structured-Unstructured Grids,” AIAA Paper 94-0645, 1994.

Steger, J .L., Dougherty, EC. and Benek, J .A., “A Chimera Grid Scheme,” ASME
Mini-Symposium on Advances in Grid Generation, Houston, Texas, June, 1983.

Steger, J .L., Dougherty, RC. and Benek, J .A., “A Chimera Grid Scheme,” Advances
in Grid Generation, K.N. Ghia eds., ASME FED, vol. 5, 1985.

Strawn, R.C., Oliker, L., and Biswas, R., “New Computational Methods for the Pre-
diction and Analysis of Helicopter Noise,” AIAA Paper 96-1696,1996.

Tam, C.K.W., “Advances in Numerical Boundary Conditions for Computational
Aeroacoustics,” AIAA Paper 97-1774, 1997.

Tam, C.K.W., “Computational Aeroacoustics: Issues and Methods,” AIAA Journal,
vol. 33, No. 10, pp. 1788-1796, 1995.

Tam, C.K.W., and Dong, 2., “Wall Boundary Conditions for High-Order Finite-
Difference Schemes in Computational Aeroacoustics,” Theoretical and Computa-
tional Fluid Dynamics, vol. 6, pp. 303-322, 1996.

Tarn, C.K.W., Kurbatskii, K.A. and Fang, J ., “Numerical Boundary Conditions for
Compositional Aeroacoustics Benchmark Problems,” in Proceedings of the Second

Computational Aeroacoustics Workshop on Benchmark Problems, ed. Tam,
C.K.W. and Hardin, J.C., 1997.

Tam, C.K.W., and Webb, J .C., “Dispersion-Relation-Preserving Finite difference
Schemes for Computational Acoustics,” Journal Computational Physics, vol. 107,
No. 2, PP. 262-281, 1993.

Tarn, C.K.W., Webb, J .C., and Dong, Z., “A Study of the Short Wave Components

in Computational Acoustics,” Journal Computational Acoustics, vol. 1, No. 1, pp 1-
30, 1993.

162

[56].

[57].

[58].

[59].

[60].

[61].

[62].

[63].

[64].

[65].

[66].

[67].

[68].

[69].

Tatsumi, S., Martinell, L., and Jameson, A., “Flux-Limited Schemes for the Com-
pressible Navier-Stokes Equations,” AIAA Journal, vol. 33, No. 2, 1995.

Tidriri, M.D., “Krylov Methods for Compressible Flows,” ICASE report No. 95-48,
1995.

Venkatakrishnan, V. “Implicit Schemes and Parallel Computing in Unstructured
Grid CFD,” ICASE report No. 95-28, 1995.

Valarezo, W.O., Momonic, OJ. and McGhee, R.J., “Multi-Element Airfoil Optimi-
zation for Maximum Lift at Higher Reynolds Number,” AIAA Paper 91-3332,
1991.

Walters, R.W., Thomas, J .L. and Switzer, G.F., “Aspects and Applications of
Patched Grid Calculations,” AIAA Paper 86-1063, 1986.

Wang, 2.]. “A Fully Conservative Interface Algorithm for Overlapped Grids,” Jour-
nal of Computational Physics, vol. 122, pp.96-106, 1995.

Wang, Z.J., Buning, P.G and Benek, J ., “Critical Evaluation of Conservative and
Non-Conservative Interface Treatment for Chimera Grids,” AIAA Paper 95-0077,
1995.

Wang, Z.J., Hariharan, N., and Chen RF, “A Fully Conservative Chimera
Approach for Structure/Unstructured Grids in Computational Fluid Dynamics,”
SBIR Phase II Final Report for NASA ARC, Contract Number: NA82-14226,
1997.

Wang, Z.J., Hariharan, N., and Chen R.F., “Recent Developments on the Conserva-
tion Property of Chimera,” AIAA Paper 98-0216, 1998.

Wang, 2.]. and Yang, H.Q., “A uniﬁed Conservative Zonal Interface Treatment for
Arbitrarily Patched and Overlapped Grids,” AIAA Paper 94-0320, 1994.

Wang, Z.J., Yang, HQ. and Przekwas, A.J., “Implicit Conservative Interfacing for
3D Overlapped Chimera Grids,” AIAA Paper 95-1683-CP, 1995.

Warming, RF. and Beam, R.M. “On the Construction and Application of Implicit
Factorized Schemes for Conservation Laws,” SIAM-AMS Proceedings, vol. 11,
pp.85-129, 1978.

Warming, RF. and Beam, R.M. “Upwind Second-Order Difference Schemes and
Applications in Unsteady Aerodynamic Flows,” Proc. AIAA 2nd Computational
Fluid Dynamics Conference, Hartford, Connecticut, pp. 17-28, 1975.

Whitharn, G. B., “Linear and Nonlinear Waves,” Wiley-Interscience, New York,

163

[70].

[71].

[72].

[73].

[74].

1974.

Yee, H. C., “A Class of High-Resolution Explicit and Implicit Shock-Capturing
Methods,” NASA TM 101088, 1989.

Zhuang, M., and Chen, R.F., “Optimized Upwind Dispersion-Relation-Preserving
Finite Difference Schemes for Computational Aeroacoustics.” AIAA Journal, Vol.
36, No. 11, 1998.

Zhuang, M., and Chen, R.F., “Applications of the Optimized Upwind Dispersion
Relation Preserving Schemes for Multi-Dimensional Acoustic Problems.” AIAA
98-2367, 1998.

Zigg, D.W., “A Review of High-Order and Optimized Finite-Difference Methods
for Simulating Linear Wave Phenomena,” AIAA Paper 97-2088, 1997.

Zigg, D.W., Lomax, H., and Jurgens, H., “An Optimized Finite-Difference Scheme
for Wave Propagation Problems,” AIAA paper 93-0459, 1993.

164