NASA/TM--2001-210985

Technique for Very High Order Nonlinear

Simulation and Validation

Rodger W. Dyson

Glenn Research Center, Cleveland, Ohio

July 2001
The NASA STI Program Office... in Profile

Since its founding, NASA has been dedicated to

the advancement of aeronautics and space

science. The NASA Scientific and Technical

Information (STI) Program Office plays a key part

in helping NASA maintain this important role.

The NASA STI Program Office is operated by

Langley Research Center, the Lead Center for

NASA's scientific and technical information. The

NASA STI Program Office provides access to the

NASA STI Database, the largest collection of

aeronautical and space science STI in the world.

The Program Office is also NASA's institutional

mechanism for disseminating the results of its

research and development activities. These results

are published by NASA in the NASA STI Report

Series, which includes the following report types:

TECHNICAL PUBLICATION. Reports of

completed research or a major significant

phase of research that present the results of

NASA programs and include extensive data

or theoretical analysis. Includes compilations

of significant scientific and technical data and

information deemed to be of continuing

reference value. NASA's counterpart of peerreviewed formal professional papers but

has less stringent limitations on manuscript

length and extent of graphic presentations.

TECHNICAL MEMORANDUM. Scientific

and technical findings that are preliminary or

of specialized interest, e.g., quick release

reports, working papers, and bibliographies

that contain minimal annotation. Does not

contain extensive analysis.

CONTRACTOR REPORT. Scientific and

technical findings by NASA-sponsored

contractors and grantees.

CONFERENCE PUBLICATION. Collected

papers from scientific and technical

conferences, symposia, seminars, or other

meetings sponsored or cosponsored by

NASA.

SPECIAL PUBLICATION. Scientific,

technical, or historical information from

NASA programs, projects, and missions,

often concerned with subjects having

substantial public interest.

TECHNICAL TRANSLATION. Englishlanguage translations of foreign scientific

and technical material pertinent to NASA's

mission.

Specialized services that complement the STI

Program Office's diverse offerings include

creating custom thesauri, building customized

data bases, organizing and publishing research

results.., even providing videos.

For more information about the NASA STI

Program Office, see the following:

• Access the NASA STI Program Home Page

at http:l/www.sti.nasa.gov

• E-mail your question via the Internet to

help@sti.nasa.gov

• Fax your question to the NASA Access

Help Desk at 301--621-0134

• Telephone the NASA Access Help Desk at

301-621-0390

Write to:

NASA Access Help Desk

NASA Center for AeroSpace Information

7121 Standard Drive

Hanover, MD 21076
NASA/TM--2001-210985

Technique for Very High Order Nonlinear

Simulation and Validation

Rodger Wo Dyson

Glenn Research Center, Cleveland, Ohio

National Aeronautics and

Space Administration

Glenn Research Center

July 2001
Acknowledgments

I would like to thank Dr. John Goodrich for introducing me to high order schemes. In addition, I thank

Dr. Ray Hixon, Dr. James Scott, and Dr. Edmane Envia for offering suggestions and encouragement.

NASA Center for Aerospace Information

7121 Standard Drive

Hanover, MD 21076

Available from

National Technical Information Service

5285 Port Royal Road

Springfield, VA 22100

Available electronically at ht_://gltrs._c.nasa.gov/GLTRS
Technique for Very High Order Nonlinear Simulation and Validation

Roger W. Dyson

National Aeronautics and Space Administration

Glenn Research Center

Cleveland, Ohio 44135

Finding the sources of sound in large nonlinear fields via direct simulation currently requires

excessive computational cost. This paper describes a simple technique for efficiently solving the

multidimensional nonlincar Euler equations that significantly reduces this cost and demonstrates

a useful approach for validating high order nonlinear methods. Up to 15 th order accuracy in space

and time methods were compared and it is shown that an algorithm with a fixed design accuracy

approaches its maximal utility and then its usefulness exponcntially decays unless higher accuracy

is used. It is conchtded that at least a 7 th order method is required to citiciently propagate a

harmonic wave using the nonlinear Euter equations to a distance of 5 wavelengths while maintaining

an ovcrall error tolerance that is low enough to capture both the mean flow and the acoustics.

1. Introduction

It is well known that aeroacoustic computations require high order numerical schemes to

emciently capture the nmltiscale l)henomena occurring over a relatively large time interval.

But the question of how high has not been answered because of a number of reasons, such

as the following:

• Acoustics may be mathematically solved in a number of ways such as via acoustic

analogies _, direct simulations 2, and perturbation techniques 3

• Commonly used numerical approaches such as Compact 4 Essentially Non-Oscillatory 5

Dispersion Relation Preserving 6 and Gaterkin 7 schemes are difficult to extend to

very high order.

• The numerical requirements vary considerably due to the wide variety of acoustical

phenomena which include linear propagation and scattering, inviscid gust interactions

with objects, receptivity effects, and turbulent noise generation.

Since most Computational Aeroacoustics (CAA) techniques have a fixed design accuracy

and resolution, they offer linfited opportunities for exploring the effects of design accuracy on

the simulation. Instead, an arbitrary accuracy Modified Expansion Solution Approximation

NASAfYM--2001-210985 1
(MESA) technique 8,9,10 is used here. This method was recently inq)lemented 1] in a way

that greatly simplified its use and demonstrated the possibility of subgrid scale resolution

by including tile spatiM derivatives of tile primitive variables on tile grid and by using

very high design accuracies 12. For these reas0ns_ t}ie MESA scheme is extended here to

the nonlinear Euler equations so that the design accuracy required for an efficient direct

numerical acoustical simulation can be quantified.

In addition, as tile accuracy requirements increase so does the complexity of the nonlinear simulator. Many subtle errors can be introduced, particularly for nonlinear systems

which have few analytical solutions available for testing with. Therefore, the Method of

Manufactured Solutions 13 is applied here and found to be a useful tool for developing and

testing very high order methods on complex nonlinear systems.

2. Problem Description

The following form of the two-dimensional nonlinear Euler equations are used 1.1,_5.

Op cOp Op ( Ou Ov _

0--/+_,_+v_+P,/_+_ =0_ (2.1)

OpOpOp (O_Ov)

CO(pro,) O(pu.___)+ O(pu 2 + p) + - Q3 (2.3)

Ox Oy

O(p_ 2 + p) 0(32) + + - (2.4)

Ot Ox COy

This form does not explicitly use the energy primitive variable since the equation of state

was used to eliminate it, resulting in only four unknown variables: p, p, u, v, which are the

density, static pressure, x-velocity, and y-velocity respectively. Tile ratio of specific heats is

defined as 3' = £a = 1.4.

CI_

Ordinarily, tile source terms (Q1, Q2, Qa, Q4) are zero. But, for testing purposes, these

terms will be modified so that the following analytical solution is defined:

p(.% y, t) = al cos(kx/rX) cos(kylry) cos(ktrrt) + c I (2.5)

u(x, y, t) = a2 cos(kxTrx) cos(kyrry) cos(ktTrt) + c2 (2.6)

v(x, y, t) = as cos(k_Trx) cos(kyrry) cos(kt_rt) + ca (2.7)

p(x, y, t) = a4 cos(kxTrx) cos(kyTry) cos(ktTrt) + c4 (2.8)

This simple harmonic solution is used to test the resolution of each method as their

design accuracy is modified and it enables comparisons to be made with previous linearized

results which also used this test 12. Tile source terms, Qi, are then found using computer

algebra and the resulting code is verified using the Method of Manufactured Solutions 16,17

Tile source terms corresponding to this analytical solution are shown in Appendix A.

It is important to insure tile density and pressure variables, p(x, y, t) and p(x, y, t), are

greater than zero when choosing the free parameters, (ai and e4), because tile numerical

NASA/TM--2001-210985 2
approachwill otherwisebenumericallyunstabledueto excessiveroundofferrorfl'omdivision

andnmltiplicationbysmallnumbers.ThelinearizedEulerequationswerepreviouslysolw_d

with a linearMESAscheme12whichdid not havethissourceof instability.However,since

p and p are not perturbation quantities in Eqs. (2.t-2.4), they will never be negative in

practical applications and therefore this instability is only of concern while imposing the

artificial solution, Eq. (2.5).

3. Numerical Approach

Tile numerical approach described here will produce single step, explicit, O = N(s + 1) - 1

(Vinteger s >_ O andN >_ 1) order of accuracy in space and time methods for the solution

of the two-dimensional nonlinear Euler equations using a N x N stencil containing the

prinfitive variables and their mixed spatial derivatives, o

OxO..._OyO...s •

The primitive variables and their spatial derivatives are advanced in time on a staggered

grid 12 such that time advancement is always performed in the center of each stencil by:

o"+bp(x'y'At) = _ 1C_b_At_' (3.9)

Ox"y b k! ' '"

k=O

Oxay b k! ' ' _

k=O

= Zc ; zxtk (3.11)

Oxay b k! ' '" k=O

oa+bp(x,Y, At) _ l Ir_Pk A÷k (3.12)

_" a_bY;_X 7- _ a,b,k,_, _

k=O "

where

y oa+b+kf(x, y, t = 0) (3.13)

C_,b, _.= Ox_ybtk

The mixed space-time (C[,b,k) derivatives are found by differentiating the governing

equations in space and time using the Cauchy-Kowalewski technique is. The difficulty in

doing this with nonlinear systems is in differentiating products and this was cited in the

past as a reason for not, using it 19. Fortunately, it can in fact be performed systematically

using a nmltidimensional extension of the Leibniz rule 2o.

oa+b+k f (x, y, t)g(x, y, t)

Ox_ybt k

b k b O(a-i)+(b-J)+(k-c) f (x, y, t) oi+j+Cg( x, y, t) ]

Ox a-i Oyb- j Otk- c Ox iOyJ Ot c J

(3.14)

By repeatedly applying the Leibniz rule to the governing equations, the explicit forms of

NASA/TM--2001-. 10985 3
the coefficients,C f

a,b,k '

DO k= I,0; DO b=O,O; DO a=O,O

are found. These coemcients are evaluated with the following loop:

_=0 i j c j=0 c=0

C_,O_,v,t) pp(x,v,_) , ,, . . ,..

C,'(x,y,O ,_,',(:r,y,O C(,(_'_J) .. . C,,(_-,y.O1 (a-i,b-j,(k-1)-c'Ji,j+l, c "-_ , - , -j,(k--1)--c i,j+l,c J

QI

-_-C¢_,b,k_ 1

_--o i j c

j=0 c=O

1)--C [+[,j,t "31-I (a_i,b_j,(_,_l)_t1,./i__l,j,c m

C,_,(z,_,,_)" C,(Z¢,_) ,.,,CpOr._.,) r,_.,(x,y.t)]

(a-i,b-j,(k-1)-c i,j+l,c -t- r (a_i,b_j,(k_l)_c'Ji,j+l,c ]

Q._

-[- Ca,b,k- l

C,p(x,y,O'_(:_,y,t)

,b,k

b _ £ £ £ a b k-1 i I jl c1

il =0 j_. _-0 ci =0 i2=0 j._=O c._=O [(,,)(j1)(Cl c2 )

(cZ._:_,___,,__,__o: c'_:_2_j_,_:__c:._,;?,2 +2c_(x.'_':

a--_ - ,o--?. _.,,l_r-- ....1)--el c_(_2__.,,__c:(:_;.L,+

• '2,.r, ,c:_ r'r (x'_"t) C'_(_:'Y't) CV(:_¢'O

_a-i! ,b-j_ ,(k- 1 )- c_ _-., -j2 ,c: -- C_+l,b,k_ 1

q_C Q3

a,b,k- 1

C,O(:_,.u.Ov(z,y,O ,b,k

a

- EEEEE _,

i_ ---0 j_ =0 c_ =0 i.>----0 j_ :-:0 c_=O

b k - 1 il j_

C"(x¢,_) C_,(x,y,O ,.-,,_,(z,_,,t)

a-i=+l,b-j,,(k-I)-c,_ i,-i_,j,-j._,c:-c2"-'i,_,j._,c._ + C:(-x'Y_)_j:,(k_l)_cl C_iY_)l,j__j_,c:_c_ C_'_(_;Yc' _) +

C:(-_':Y_)-yl,(k_l)_c] Ci:(x--}Y:_l_j._,cz_c2 C_'_'YlJ_,c.2 ) ] -- cPa,b+ i,k_l

+c_,__,

END DO; END DO

Perform Conversion Eq. (3.16)

END DO (3.15)

NASA/TM--2001-210985 4
DO a,b (a + b = 0, 1,...,20) and(O <_ a,b < O)

Co_.o_o _,_,_--

, , o,o,c_ \i-O j-O

a-1 k /

F_E i-O c-O t "i

b-lk(b ZZ 5

j-0 c-0

a b k CO t_±

c-o i j c a-i'b-j'k-c_Ji'j'cl

C Ca-i,O,k-cCi'b,c +

C O'b-j'k-cWa'j'c +

C[_' b,k =

,,lb1( i-o j-o i )(j )Cta-i,b-J 'OCi'j'k

c--O C O'O'k--c_a"b'e]

o,0,u ' ' o,o.o \i--0 j-0

a--I k (

"-- c--0 i

b-lk ZZ (b j

j--O c--O

+

)()() Z a b k (To f,v a_

c-o i j c "Ja-i,b-j,k-c"Ji,j'c[

C _a-i'O'k-c_i'b'c -}-

C C;'b-j'k-cCa'j'c +

i--O j-O i j Wa-i'b-j,OWi'j'k

_-_o c C°,°, k-_ C_''b'_

+

END DO (3.16)

After each loop of the k index in Eq. (3.15), the following momentum coefficients are

known:

CP,_b,k with (a,b): (a,b=0,1,2,...,O)

cPa_b,#with (a,b):(a,b=0,1,2,...,O) (3.17)

But the C_,b,k and C_',b,k coefficients are required to evaluate the coefficients, C_..o,o...,o,+1

at the next k + 1 loop. By repeatedly applying the Leibniz' rule to Eq. (3.17) an explicit

expression for the C,_,b,k and CV_,b,k coefficients is found, Eq. (3.16), in terms of already

known variables. Note, however, that these coefficients must be evaluated in the sequence

NASA/TM_2001-210985 5 "
of increasingderivativeorder(i.e.,Find C_,,b,_,and Ci, b,_. in the order (a,b):(a+b=0,1,..., 2

O). By apl)tying this conversion step once for each k-loop increment in Eq. (3.15), we can

efficiently continue increasing the index k until all the coefficients, Eq (3.13), are found.

4. Determining the Spatial Derivatives

The procedure described in the last section will advance the nonlinear Euler equations in

time at (x,y) with arbitrary accuracy if all the spatial derivatives:

cf, b,0V(,,, b) : a, b = 0, 1,..., O

are known at (x, y). However, only the spatial derivatives:

C fb,oV(a,b) : a,b = O, 1,...,s

are available on each grid point, (xi,y)). Therefore, the higher order spatial derivatives,

(a, b _> s + 1), need to be interpolated. In addition, since time advancement is performed

using a staggered grid, both the known derivatives on the grid and the higher order unknown

derivatives need to be interpolated to the stencil center, (x = 0, y = 0). A simple method

for accomplishing this on a N = 2 size stencil was shown previously 12. It is completely

generalized here to include any size stencil (N > 1) using any number of spatial derivatives,

s _> 0, and at any location such as on solid boundaries which require one-sided stencils.

A one-dlmensional osculating spatial interpolant using the primitive variables and their

spatial derivatives at the grid points, (x0, xl,..., XN-,), is:

N-1 k(s+l)+s

CL_ _(_) o,o,o = _ _--_

k=O i=k(s+l) Qi,i(j__o(Xk-1 ) s+i -- xj) (X -- Xk),_i-k(s+l)"

where the Qi,i are computed as follows:

(4.18)

DO k = O,N- 1; DO i = k(s + 1),k(s + 1) + s;

DO j = O,i - k(s + 1)

Qi,j = (1/j!)_

END DO; END DO; END DO

(4.19)

DOm=I,N-1; DOk=I,N-m

DO i = (s + l)(m + k - 1),(k + m)_ + m + k- 1

DO j = i+ 1 - (s + 1)k,i + s + 1 - k(s + 1) (4.20)

(Qi,.i- i -Q_-:,j-1 )

Qi,j = max

END DO; END DO; END DO; END DO

The real value of this algebraic form is the simplicity of interpolating the derivatives;

Ordinarily, by the product rule of differentiation, each differentiation of Eq (4.18), would

result in an equation that grew in size by a factor of N, up to N ° terms, and would render

NASAITM--2001-210985 6
very high ordermethodsintractable.Instead,the followingexplicit algebraicform avoids

this issue:

N-1 k(s+l)+s

cf,,(x) ,0,0= E E (4.2L)

k.-0 i=k(s+l)

where tile function Z_,i,_. is independent of space and time and can be computed simply as

follows:

DO a=O,(N- l)(_'+l)+s

DO i = O, s

a--1

z,,,_,o= (:,,- *o)(;-a) 1-[(i - _)

c--O

END DO

END DO

DO (, = O,(N - 1)(, + 1)+,

DOk=I,N-1

DO i= /c(s ÷ l),k(s ÷ l) ÷ s

Z,_,i,_o= (4.22)

zz... z

tO--0 r! =0 rk_:=0 \_'=0 L\ /A

{ 1

(a-E___-o I rp)-I (rk_i_l) k-1 (r/_-1-1) )

[*+ 1- _01 1-[ [i- k(._+ l) - _k]I] II [._+ l - _e?

el¢=0 fl=l e_=0

• II

eo=O

END DO

END DO

END DO

This technique is an extension of Hermitian divided-difference interpolation 21 and it

was derived by repeatedly applying the generalized Leibniz' rule in Eq. (3.14). For a fixed

stencil size on a uniform grid, Eq. (4.22) needs to be computed only once since it is constant

for all space and time. This results in significant computational savings since only the Qi,i

variables require further evaluation as the stencil location changes, but even these stay

constant as higher order derivatives are evaluated using Eq. (4.21) on the same stencil.

In addition, higher order interpolants may reuse the information from lower order interpolants to efficiently and adaptively increase local accuracy:

(X-1)(s+l)+s

CfN (_) ,.-,IN-_ (x) (4.23) ,0,0 = Vla,O,O ---- E Qi,iZa,i, k

i=(N-1)(,+1)

Multidimensional spatial interpolations are efficiently accomplished using the tensor

product of a series of one-dimensional interpolants 11a2

NASAfFM--2001-210985 7
5. Explicit forms for the Analytic Solution and its Source Terms

Coml)uter algebra can calculate tile derivatives of tile analytical solution, Eq. (2.5), and

Q,

the source terms, Ca,b,k_ l in Eq. (3.15), and then automatically write tile necessary FORTRAN code 23,13. However, for very high order algorMtms the code becomes too large and

complicated for most compilers and therefore it is necessary to develop an explicit form of

those expressions.

First define the Lagrangian operator as:

(l-IT_-o1 .rv[od(a, 4) - i) (I_a_,+l glod(a, 4) - i)

L,,,, (5.24)

And then define:

CCO._(¢z)

d =

Cdin(¢ z) =

(Ld,o eos(¢z) -- Ld,1 sin(¢z) -- Ld.2 cos(¢z) + Ld,a sin(¢z))¢ d (5.25)

(nd,0 sin(¢z) + nd,1 cos(¢z) -- nd,2 sin(¢z) -- Ld,a cos(¢z))¢ d (5.26)

where ¢ is any constant, d = a, b, or, k, and z = kxrrx, kyrry, or, ktTrt.

The derivatives of the assumed analytical solution, Eq. (2.5), may then be written as:

c::;:z'

'aanal

,b,k

C (L'Lra r_ d t

,b,k

Ca pana!

,b,k

,- _a ,- _b_- _k--,cos(k=*rz)_cos(k utry)_cos(k_ rrt)

= 5,_05b05_0cl + alttq-rr) t_cyrr) (_,rr) wa t-.b wk (5.27)

,- ,a,- ,b,- ,k.-,cos(k=_x).-,co._(ku_y).-,cos(_:_=0 (5.28) --_ (_aOSbOSkOC2 -'[- a2itCxTr) (I,CyTr) (I, gtTr ) (Ja (_b (Jk

= 5aOSbOSkOCa+ aa(rcxTr) ,- ,a,-_yrr) ,b,- ,k--,cos(k_ _z) .-,cos(k_ry) ..,cos(kt _t) (5.29) t_trr) w_ _b %.

,- ,a,- ,b,- ,k_cos(k_rrz),-,cos(k_rrU),-,cos(k,._O (5.30) = 5aOabOSkOe 4 + a4(tcxrr) (tcyrr) ttctTr) t_a t_ b t_ k

where ai0 is the Kronecker delta.

Tile source terms involve products of transcendental functions which may be differentiated explicitly using the Leibniz' rule, Eq. (3.14), and are given by:

cdtrigl (¢I z)trig2 (¢2 z)

c dtrigI (¢1 Z)trig2 (¢._z )trig3 ( ¢a z )

= d()d C'd_trig=(¢'_z) i C_trig.2(¢.2z)

i=0 i

= Z Z d i i_trig1(¢Iz)ctrig2(¢2z)ctriga(¢az )

i--0 i2=0 i i2 "_d-i i-J2 i2

where trigi is either the cos or sin function and ¢i is some constant.

Tile explicit forms of the derivatives of the source terms using these definitions are listed

in Appendix A. In this way, only a few pages of FORTRAN code are required to compute

all the necessary very high order derivatives of the analytical solution and source terms.

The Method of Manufactured solutions permits any analytical solution to be imposed,

but note that not all solutions and their induced source terms will be explicitly differentiable

as they are here. Therefore, additional care is required in the selection of analytical solutions

for testing very high order methods for otherwise the resulting code will become unwieldy.

NASA]TM--2001-210985 8
6. Stability

Derivinganalyticalstability conditionsfor the nonlinearEuler equationsis generallydifficult 24andnecessarystability conditions(but generallynot suffmient)aretypicallyderived

by linearizingthe conservativegoverningequationsusingtheir constantcoefficientquasilinearform:

OU A OU B OU = o (6.31)

0---(+ Ox + Oy

where A and B are the flux Jacobians which are held constant for a yon Neumann stability analysis. However, in this work, the energy equation is not in conservation form and

therefore this form of analysis is not available.

Instead, engineering rules of thumb were applied here and tested numerically. Let the

Courant number be defined by 25.

cr = + At (6.32)

where Cx = max([lull + c), which is tile magnitude of tile greatest wave velocity in the xdirection over the entire grid (and similarly for cy in the y direction). In addition, for most

single-step explicit nmnerical formulations in two-dimensions the Courant number nmst be

a _< 0.5 to maintain numerical stability 26

Since the Cauchy-Kowalewski procedure is essentially a Taylor series method, a simple

model problem may be used to demonstrate numerical stability limits. For example 27, if

the test equation (for)_ < 0):

!

y = ,_y (6.33)

is time advanced using an n th order Taylor series in time, it will be stable if the following

series (which is generally convergent as n.--+ ec) is satisfied:

I11+Ate+ (1) At2a 2+...+ (I)_ At'_'_[[< 1 (6.34)

At low accuracy the maximum stable time step will vary, but it soon approaches a fixed

value as the accuracy increases. In general, however, high-order explicit schemes usually

have more severe stability restrictions on At 2s and this is evident in Table 3.

The following time step size was found to provide numerically stable solutions for up to

15th order accuracy on a uniform Cartesian grid:

(7

/,t = (6.35)

where the speed of sound is defined by :

/ p(x,y,t) _

c= v

(6.36)

NASA/TM--2001.-210985 9
Numericalinstabilitydid occurwhenthesolutionwasunder-resolvedasshownin Table

6. Here,the 3ra ordermethodbecameunstableby the time t = 799.98, but the 5th, 7 th,

and 9 lh order methods were stable to at least 50 wavelengths of propagation. This suggests

that nonlinear aliasing effects can be delayed significantly by increasing the accuracy of the

simulation. Similar results are seen in Table 3.

Ill practical applications, it will be necessary to filter out. the inevitable unresolved high

frequencies. However, it is encouraging to observe that despite the very low numerical

dissipation, filters were not required in this work.

7. Numerical Results

A comparison of the effects of design accuracy on the efficiency and resolution of the overall

simulation wa._ performed by using a 2 x 2 stencil to produce 28 + 1 order schemes. These

schemes were used to propagate a single wavelength the following distances:

• A single time step as shown in Tables 1-2,

• Five wavelengths as shown in Tables 3-5,

• Fifty wavelengths as shown in Table 6.

The wavenumbers, amplitudes, and displacements defined in Eq. (2.5) were set to ki = 1,

ai = 1 and ci = 2, respectively, except in Table 6 which has ci = 10.

The expected error, Ys, of a method with design accuracy O = 28+ 1 can be extrapolated

using the data in Table 1 since the design accuracy determines the rate of error change as

the grid is refined. Starting with the known numerical error, es, in Table 1 with Xs = 2, the

following relationship:

_ es (7.37)

Ys (20(x__2))

Gridpoints will provide the expected error at Xs = ILog2 IVavelength t" Figure 1 plots this equation.

Note that, due to computer roundoff, errors less than 10 -15 are only possible if quadruple

precision floating point numbers are available.

By solving for Xs in Eq. (7.37), it is possible to estimate the number of grid points

required to meet a particular error level Ys.

log 2 es - log 2 Ys

+ 2 (7.38)

Xs _--- O

The total number of grid points in a two-dimensional unit domain required for propagating

a single wavelength is then given by (2x_ + 1) 2. In addition, each grid point contains the

four prinfitive variables p, p, u, and v as well as their, 0 2, spatial derivatives resulting in a

total memory cost of:

costmemory = 4(28 + 1)2(2 z_ + 1) 2 (7.39)

NASAITM--2001-210985 10
Table 1. Numerical error at t=0.05, cr=0.25, 4 grid points per wavelength, 25 grid points

ci = 2, ai = 1,hi = 1, At = .05, Single Step

s Order Error, es Seconds/Grid Point, cs

1 3 1.57 10 -2 .0012

2 5 2.74 10 .4 .015

3 7 2.5410 .6 .1172

4 9 1.44 10 -8 .678

5 11 5.40 10 -11 3

6 13 1.42 10 -'3 11

7 15 1.78 10 -:_ 36

8 17 3.55 10 -16 95

9 19 roundoff 238

Table 2. Numerical error at t=0.05, a=0.5, 8 grid points per wavelength, 81 grid points

ci = 2,ai = 1,k_ = l,,_kt = .05, Single Step

s Order Error, es Seconds/Grid Point, cs

1 3 9.55 10 .4 .0016

2 5 3.95 10 -6 .0183

3 7 1.18 10 -s .1449

4 9 3.46 t0 -:_" .839

5 11 4.85 10 -t4 3.7

6 13 roundoff 13

Table 3. Numerical error at t=10, a=0.5, 4 grid points per wavelength, 25 grid points

ci = 2, ai = 1, k_ = 1, At = 2.4755 10 -2, 403 Steps, 5 wavelengths

s Order Error, e_ Total Seconds Seconds / Grid Point, cs

1 3 UNSTABLE UNSTABLE UNSTABLE

2 5 1.13 10-: 188.22 7

3 7 1.67 10 -3 1497.53 59

4 9 1.50 10 .5 8610.51 344

5 11 5.21 10 -6 39125.30 1565

6 13 UNSTABLE UNSTABLE UNSTABLE

Table 4. Numerical error at t=10, a=0.5, 8 grid points per wavelength, 81 grid points

ci = 2, ai = 1,k, = 1, At = 1.2377 10 -_, 807 Steps, 5 wavelengths

s Order Error, es Total Seconds Se_'onds / Grid Point, c_

1 3 1.4710 -U 166.90 2

2 5 3.95 10 .3 1392.09 17

3 7 9.73 10 .6 13328.04 164

4 9 1.64 10 -8 64907.56 801

5 ]1 1.45 10 -11 275939.53 3406

6 13 8.29 10 -13 999763.10 12342

NASA/TM--2001-210985 11 ....
Table 5. Numerical error at t=10, 0=0.5, 16 grid points per wavelength, 289 grid points

ci : 2, ui = 1,k, = 1,_t = 6.18810 -3, 1615 Steps, 5 wavelengths

s Order Error, e_ Total Seconds Seconds / Grid Point, c_

1 3 1.48 10-1 935 3.2

2 5 2.54 10-4 10198 35

3 7 2.02 l0 -T 80299 277

4 9 1.29 10-_° 463537 1603

Table 6. Numerical error at t=100, _r=0.5, 8 grid points per wavelength, 81 grid points

ei = 10,ai = 1,k, = 1,At = 5.077 10-a, 19692 Steps, 50 wavelengths

s Order Error, es Total Seconds Seconds / Grid Point, c,

1 3 LLNSTABLE UXSTABLE UNSTABLE

2 5 1.05 10-a 65407 807

3 7 8.65 10-6 406932 5023

4 9 1.08 10-_ 1662013 20518

In addition, tile total time per grid point for each method, cs is shown in Table 1. However,

this time reflects the cost of a single time step. As ms increases (grid density) the maxinnml

stable time step will decrease for explicit methods (to maintain CFL condition). This is

accounted for by including the factor, 2 x*-2, for a total time cost of:

costtime = cs2Zs-2(2 z_ + 1) 2 (7./10)

Fig. 2 compares the relative costs of propagating a wavelength over a range of error tolerances, 10 -15 _< y, _< 10°.

Generally, a direct acoustical sinmlation requires an error of less than 10 -_ to capture

both the acoustical and mean flow phenomena. By this standard, the correct propagation of

a wave cannot be realistically achieved with a 3_a order method because of the exponential

growth in its cost as shown in Fig. 2. In addition, additional error will accumulate at each

time step and as Tables 3-5 indicate, a higher order method is required to propagate out to

a distance of 5 wavelengths. As Table 6 shows, at least 7 th order accuracy is required for

propagating 50 wavelengths with 8 grid points per wavelength.

8. Conclusions

The procedures shown here provide arbitrary accuracy in both space and time for nonlinear

Euler simulations using only a few pages of FORTRAN code. The usefulness of these

procedures hinges on the required accuracy of a given sinmlation.

If an accuracy of 10 -1° is required at time t = 0.05 to maintain a given error tolerance at

some future time, then Fig. 2 indicates the 3rd and 5th order methods could not accomplish

this because of exponential cost growth; Here, the optimal design accuracy is 11 th order.

tn addition, if the error tolerance is decreased to 10-15 because of the desire to run the

sinmlation over a larger time period, it is most optimal to use a 15_'h or higher order

method.

NASAITM--2001-210985 12
0

-10

z

e_

<

E -15

-20

O

-25

2

-1

Relative Wave Resolution by Grid Points Per Wavelength

',:--:?----2, i_.......

"- 2- t

''*, . .

0

l ....

"'''--.., ..

2 3

ILog 21of Grid Points Per Wavelength

-*....

-2/

"4t

z

4 5

Figure 1: Resolution Studies

ORDER

--+- 3rd

--*- 5th

l- 7th

- _ - 9th

-* llth

+ 13th

-*- 15th

These results show that a fixed order algorithm approaches its maximal utility and then

quickly becomes useless. The nonlinear implementation of the MESA schemes described

here can be used to adjust the design accuracy arbitrarily and it offers an opportunity for

solving the more demanding acoustical sinmlations in smooth flows.

NASA/TM--2001-210985 13
p.

3:

o

r..)

[..,

_a

o

3:

E

o

r..)

eD

2000

1500

1000

5O0

6

1.10

800000

600000

400000

200000

Relative Time Cost of Wave Resolution by Error Requirement

/

/

I

/

I

I

/

/

I

/

/

/

/

"1"

0 2 4 6 8 10 12 14

-Log 10 Error

Relative Memory Cost of Wave Resolution by Error Re_ uirement

: I

." /

," /

I

I

/'

I'

/

I

I'

__ .... _-_ -. .-__-._-JC-__----:-_

0 2 4 6 8 10 12 14

-Log 10 Error

ORDER

3rd

*- 5th

i-7th

-A-- 9th

-'- llth

--*- 13th

,* 15th

ORDER

+3rd

-*- 5th

-i 7th

-A--9th

-" llth

-*- 13th

-*- 15th

Figure 2: Time and Memory Efficiency Comparison For Errors Up To 10 -15

NASA/TM--2001-210985 14
Appendix A

SourcestermsQ_, Q2, Q3, and Q4 in Eqs. (2.1, 2.2,

are defined below using the definitions in section 5:

2.3, 2.4) and their explicit derivatives

a a, k _rk •C c°s(kt nt)_/-_cos(k_ n'y) _/-_cos(kx_-x) sin(kz _rx)

2 ",_ "x xyr k "% _'a

2_,2_ t_ --L r-,cos(ktnt)_F_cos(ku_rY)2r_COs(k_rx)sin(k_rx )

'_2 bl r_x n _Vxyt t_ k "J b tJa

4a_ a_ c_ kz T_kx_t C_ °s(kt _-t)_c2OS(k_ _y): CaCOS(G_rx)sin(k_ _'x) _

.... 21_ _I,. f'_c°s(ktrct)3f_c°s(kv_rY) a cos(kx_rx) cos(kzr_x)sJn(kx_rx)

t+ i"_2r_x ,, r_xyt _" k _ b Ca

NASA/TM_2001-210985 15
a,c_c.,k 7rJv .¢_c°s(ktm)_c°s(kz_x)/_sin(kv_Y)

'* "_ 'o y ,xyrt/k tJ u tjb __

I. " .--,cos(ktrrt)-.--,cos(k.rrx) .._,cos{k.,rcy)sin(kvrry )

2a2a'3(q t_yzr_xytt, k t_,_ b b , - _

• . cos kt _rt) 2F, cos(kz zrx) 2 f_cos(kv :,ry) sin(k v _ry) 2alar3(2kyTrkxytC_, ( "-"a _b --

, cos(ktTrt _

2ala2(akvTrk:rytC'_2 ) ccos(t'_X)2CbOS(k_y)sin(l,._u)

3a j. a2a3ky.rrkxytC2°S( k' _t)3Cc°s( k_z)3c2os(_',j _y) cos(J,°v_ry) sin(1,'_ry)

Q4 [c_7rk . ¢_cos('k. rrx) ,.-,cos(k v _ry5_sin (k, rrt)

C_,b,k = --a3cl r- xyr_a _b _"tk --

a . __ ¢-,cos(kzrrx)pcos(kvTry)z-_in(ktTrt )

1(:3 _:t J_r_xyt tJa _ b C,'L_ --

2a _ *. --_ _cos(k_ 7rx) 2 (_TCOS(kun.y)2 ¢_cos(ktwt) sin(kt rrt)

.] U_3tirt ti Exytt_ a "_b t_ k --

(1 C C k 7r]i" cC°s(ktnt)f'c°s(kvlry)_sin(kz 7rx)

_3 '1 2 x 'xyt k tJb t_a --

a2cl C3 kx 751_'xytC; °s(kt nt) (-,cos(k v Try)(_sin(kz _x) "_b w. --

a.p,,c,,b _I,. .g-ycos(ktrrt)pcos(kyrry) sin(kzrrx)

*_z d'_'x', ,_xyt_ k _b _ --

2a2a3cl ]% 7rk'xytC2 °s(ktm)2 c;°S(kv _rY): caCOS(k. _rx) sin(kx ,Tx) _

9o. 0 ..... t. _,t. .f'_c°s(kt_rt)_r-_cos(ku_ry) _ cos(kz_rx)sin kzrrx __ " '* 'd"z".x,, '_xy t "_'k 'J b Ca ( )

2a_ a.,c..k- 7rk . cCOS(k_ _t) _,_cos(k v _rv)_ _cos(k_.x) sin(k_ _z)

• _ o x Xyt k t_" b t-_tt

3a.a_,_..t. ,_t. -(wc°s(kt_rt)3f-'c°s(kvrrY)3pc°s(kzrrx)cos(k.zrx sin(kznx)

z"3'_x .... xyt _l k _'_ b t_ a ) --

- ['-,¢os(kt rr_').--,cos(kz nx) _sin(k ny _

a4 _y 7r_'xyt _,to _,a _ c. b v )

2a c c k 7rk _cos(kt_rt)f_cos(kzwx)f, sln(kvrry)

3 4 3 y xyttJ k t_ a L, b

a_ c_ k.. 7rk .... c.COS(_t,rt) f_cos(k_ _rx) f, sin(kv,_)

o v ._yt k '.-/o _b --

2a_c_ k- 7rk .... _cos(k_ _rt) _,,_cos(k_ _rx) _pcos(k v _ry) sin(k_ry)

"l y *Y" -- If t'_a "-Jb

,, ...... ,, l_ _'b (Tcos(kt nt) _ (-,cos(k_ nx) 2 _cos(k v _ry) sln(k v try)

References

1. Lee, J., Cho, K., Lee, S. "Application of acoustic analogy to automotive engine-cooling fan noise

prediction", AIAA J., Vol. 38 No. 6 pp. 1095-1098 JUN 2000

2. Colonius T, Lele S.K., Moln P. "Sound generation in a mixing layer" J. of Fluid Mechanics, \hi.

330 pp. 375-409 JAN 10 1997

3. Morris P.J., Long L.N., Bangalore A., Wang Q.Z. "A parallel three-dimensional computational

aeroacoustics method using nonlinear disturbance equations" J. Comp. Phys., Voi.133 No.l, pp.

56-74 MAY 1 1997

4. Povitsky, A., and Morris, P.J., "A Higher-Order Compact Method in Space and Time Based On

Parallel Implementation of the Thomas Algorithm", J. comp. Phys., i61, 182-203 (2000)

5. Yee, H.C., Sandham, N.D., and Djomehri, M.J. "Low-dissipative High-Order Shock-Capturing

Methods Using Characteristic-Based Filters", J. Comp. Phys., Vol. 150, 199-238, 1999.

NASA/TM_2001-210985 16
6. Tam, C.K.W. "Computational Aeroacoustics: Issues and Methods", AIAA J., \:ol. 33, No. 10,

Oct. 1995.

7. Harari, I., Avraham, D., "High-order finite element methods for acoustic problems", J. Comp.

Acoust., Vo]. 5: (1) pp. 33-51 .MAIR 1997

8. Goodrich, J.W. "Accurate Finite Difference Algorithms for Computational Aeroacoustics", Computational Fluid Dynamics Review 1998, edited by M.M Hafez and K. Oshima, World SeieT_tific,

Singapore, 1998.

9. Goodrich, J. W., "A Comparison of Numerical Methods for Computational A eroacoustics", A fA A

99-1943, May I999.

10. GoodT"ich, .L W. "A Low Order and a High Order Solution for a Converging-Diverging Nozzle

Problem", NASA/CP-2000-209790, pp. 265-271, Aug. 2000.

11. Dyson, R.W., Goodrich, J.W. "A Very High Order, Adaptable MESA Implementation for Aeroacoustic Computations", NASA/TM-2000-209944, April 2000.

12. Dyson, R.W., Goodrich, J.g_, "Automated Approach to Very High Order Aeroacoustics Computations", A[AA J., March 2001.

13. Roache, P.J., "Verification of codes and calculations", AIAA J., Vol. 36: (5)pp. 696-702, MAY

1998

14. Laney, C.B., Computational Gasdynamics, , Cambridge University Press, 1998, pp. 9-10.

15. Thompson, K.W., Time-Dependent Boundary Conditions for Hyperbolic Systems, H, J. Comp.

Phys., Vol. 89, pp. 439-461, 1990

t6. Kambiz Salari, Patrick Knupp , Code Verification by the Method of Manufactured Solutions,

Sandia National Laboratory Report, SA.N'D2000-1444, June 2000

17. Roache, P.J., Verification and Validation in Computational Science and Engineering, Hermosa

Publishers, I998, pp.65-90

I8. Laney, C.B., Computational Gasdynamics, , Cambridge University Press, 1998, pp. 316-323.

19. Garabedian, P.R., "The Method of Power Series: The Cauchy-Kowalewski Theorem," Partial

Differential Equations, 2"ded., Chelsea, New York, 1986, pp. 16.

20. FuIks, W, Advanced Calculus: An introduction to analysis, John Wiley & Sons, 1978, pp. 109.

21. Powell, M.J.D., "Hermite Interpolation," Approximation Theory and Methods, Cambridge University Press, 1981, pp. 53-57.

22. Botha, J.F., and Pinder, G.F., "Two-Dimensional Polynomial Approximation," Fundamental

Concepts in the Numerical Solution of Differential Equations, Wiley, New York, t983, pp. 17-

22.

93. Akers R.L., Baffes P., Kant E., Randall C., Steinber9 S., Young R.L. "Automatic synthesis of

numerical codes for solvin9 partial differential equations", Math. Comp. in Simulation, Vol. 45

No. 1-2, pp. 3-22, JAN 31 1998

24. Hirseh, (7_..Numerical Computation of Internal and External Flows: Volume 2: Computational

Methods for Inviscid and Viscous Flows, John Wiley & Sons, 1988, pp. 379-384.

25. Thompson, K.W., Lecture Series in Computational Fluid Dynamics, NASA TM I00010, August

1987.

26. Hoffman, K.A., Chiang, S.T., Computational Fluid Dynamics For Engineers - Volume II, Engineering Education System, 1993, pp. 15.

27. Burden, R.L., Faires, J.D. Numerical Analysis, Fourth Edition, PWS-Kent Publishing Company,

Boston, 1989, pp. 304.

28. Fletcher, C.A.J. "Computational Techniques for Fluid Dynamics 1: Fundamental and General

Techniques", Springer- VerIag, Berlin, Heidelberg, 1991,

NASA/TM--2001-210985 17
REPORT DOCUMENTATION PAGE FormApproved

OMB NO. 0704-0188

Publicreportingburdenforthiscollectionofinformationis estimatedtoaverage1 hourperresponse,includingthetimeforreviewinginstructions, searchingexistingdata sources,

gatheringand maintainingthe dataneeded,andcompletingandreviewingthecollectionof information.Sendcommentsregardingthisburdenestimateor any otheraspectofthis

collectionof information, includingsuggestionsforreducingthisburden,to WashingtonHeadquartersServices,DirectorateforInformationOperationsandReports,1215Jefferson

DavisHighway,Suite1204,Arlington,VA 22202-4302,andtotheOfficeofManagementandBudget,PaperworkReductionProject(0704-0188),Washington, DC 20503.

1. AGENCY USE ONLY (Leave blank) I 2. REPORT DATE 3. REPORT TYPE AND DATES COVERED

I July 2001 Technical Memorandum

4. TITLE AND SUBTITLE

Technique for Very High Order Nonlinear Simulation and Vafidation

6. AUTHOR(S)

Rodger W. Dyson

i 7. PERFORMING ORGANIZATION NAME(S) AND ADDRESS(ES)

National Aeronautics and Space Administration

John H. Glenn Research Center at Lewis Field

Cleveland, Ohio 44135- 3191

9. SPONSORING/MONITORING AGENCY NAME(S) AND ADDRESS(ES)

National Aeronautics and Space Administration

Washington, DC 20546-0001

5. FUNDING NUMBERS

WU-781-30-- 11-0

8. PERFORMING ORGANIZATION

REPORT NUMBER

E-12839

10. SPONSORING/MONITORING

AGENCY REPORTNUMBER

NASA TM--2001-210985

11. SUPPLEMENTARY NOTES

Responsible person, Rodger W. Dyson, organization code 5940, 216-433-9083.

12a. DISTRIBUTION/AVAILABILITY STATEMENT

Unclassified - Unlimited

Subject Category: 71 Distribution: Nonstandard

Available eleclronically at http://gltrs._c.nasa.gov/GLTRS

This publication is available from the NASA Center for AeroSpace Information, 301--621--0390.

12b. DISTRIBUTION CODE

13. ABSTRACT (Max/mum 200 words)

Finding the sources of sound in large nonlinear fields via direct simulation currently requires excessive computational

cost. This paper describes a simple technique for efficiently solving the multidimensional nonlinear Euler equations that

significantly reduces this cost and demonstrates a useful approach for validating high order nonlinear methods. Up to

15th order accuracy in space and time methods were compared and it is shown that an algorithm with a fixed design

accuracy approaches its maximal utility and then its usefulness exponentially decays unless higher accuracy is used. It is

concluded that at least a 7th order method is required to efficiently propagate a harmonic wave using the nonlinear Euler

equations to a distance of 5 wavelengths while maintaining an overall error tolerance that is low enough to capture both

the mean flow and the acoustics.

14. SUBJECT TERMS

High order; Nonlinear Euler; Multidimensional; Acoustics

17. SECURITY CLASSIFICATION 18. SECURITY CLASSIFICATION

OF REPORT OF THIS PAGE

Unclassified Unclassified

NSN 7540-01-280-5500

19. SECURITY CLASSIRCATION

OF ABSTRACT

Unclassified

15. NUMBER OF PAGES

23

16. PRICE CODE

20. LIMITATION OF ABSTRACT

Standard Form 298 (Rev. 2-89)

Prescribedby ANSI Std. Z39-18

298-102