John W. Goodrich

Glenn Research Center, Cleveland, Ohio

Error Control With Perfectly Matched Layer or 

Damping Layer Treatments for Computational 

Aeroacoustics With Jet Flows

NASA/TM—2010-216092

January 2010
NASA STI Program . . . in Profi le

Since its founding, NASA has been dedicated to the 

advancement of aeronautics and space science. The 

NASA Scientifi c and Technical Information (STI) 

program plays a key part in helping NASA maintain 

this important role.

The NASA STI Program operates under the auspices 

of the Agency Chief Information Offi cer. It collects, 

organizes, provides for archiving, and disseminates 

NASA’s STI. The NASA STI program provides access 

to the NASA Aeronautics and Space Database and 

its public interface, the NASA Technical Reports 

Server, thus providing one of the largest collections 

of aeronautical and space science STI in the world. 

Results are published in both non-NASA channels 

and by NASA in the NASA STI Report Series, which 

includes the following report types:

• TECHNICAL PUBLICATION. Reports of 

completed research or a major signifi cant phase 

of research that present the results of NASA 

programs and include extensive data or theoretical 

analysis. Includes compilations of signifi cant 

scientifi c and technical data and information 

deemed to be of continuing reference value. 

NASA counterpart of peer-reviewed formal 

professional papers but has less stringent 

limitations on manuscript length and extent of 

graphic presentations.

• TECHNICAL MEMORANDUM. Scientifi c 

and technical fi ndings 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. Scientifi c and 

technical fi ndings by NASA-sponsored 

contractors and grantees.

• CONFERENCE PUBLICATION. Collected 

papers from scientifi c and technical 

conferences, symposia, seminars, or other 

meetings sponsored or cosponsored by NASA.

• SPECIAL PUBLICATION. Scientifi c, 

technical, or historical information from 

NASA programs, projects, and missions, often 

concerned with subjects having substantial 

public interest.

• TECHNICAL TRANSLATION. Englishlanguage translations of foreign scientifi c and 

technical material pertinent to NASA’s mission.

Specialized services also include creating custom 

thesauri, building customized databases, organizing 

and publishing research results.

For more information about the NASA STI 

program, see the following:

• Access the NASA STI program home page at 

http://www.sti.nasa.gov

• E-mail your question via the Internet to help@

sti.nasa.gov

• Fax your question to the NASA STI Help Desk 

at 443–757–5803

• Telephone the NASA STI Help Desk at

 443–757–5802

• Write to:

 NASA Center for AeroSpace Information (CASI)

 7115 Standard Drive

 Hanover, MD 21076–1320
John W. Goodrich

Glenn Research Center, Cleveland, Ohio

Error Control With Perfectly Matched Layer or 

Damping Layer Treatments for Computational 

Aeroacoustics With Jet Flows

NASA/TM—2010-216092

January 2010

National Aeronautics and

Space Administration

Glenn Research Center

Cleveland, Ohio 44135

Prepared for the

15th Aeroacoustics Conference (30th AIAA Aeroacoustics Conference)

cosponsored by AIAA and CEAS

Miami, Florida, May 11–13, 2009
Acknowledgments

This research has been supported by the Supersonic Project of NASA’s Fundamental Aeronautics Program. These results are 

presented with the expectant hope that they will help produce more accurate numerical simulations that can lead to greater 

insight into aeroacoustic phenomena, particularly to jet noise, and subsequently to quieter jet aircraft. A copy of the illustrated 

talk that accompanies this paper can be obtained as a PDF fi le by contacting the author (John.W.Goodrich@nasa.gov).

Available from

NASA Center for Aerospace Information

7115 Standard Drive

Hanover, MD 21076–1320

National Technical Information Service

5285 Port Royal Road

Springfi eld, VA 22161

Available electronically at http://gltrs.grc.nasa.gov

This work was sponsored by the Fundamental Aeronautics Program 

at the NASA Glenn Research Center.

Level of Review: This material has been technically reviewed by technical management. 

This report is a formal draft or working 

paper, intended to solicit comments and 

ideas from a technical peer group.

This report contains preliminary fi ndings, 

subject to revision as analysis proceeds.
1 Abstract

In this paper we show by means of numerical experiments that the error intruduced in a

numerical domain because of a Perfectly Matched Layer or Damping Layer boundary treatment can be controlled. These experimental demonstrations are for acoustic propagation

with the Linearized Euler Equations with both uniform and steady jet flows. The propagating signal is driven by a time harmonic pressure source. Combinations of Perfectly Matched

and Damping Layers are used with different damping profiles. These layer and profile combinations allow the relative error introduced by a layer to be kept as small as desired, in

principle. Tradeoffs between error and cost are explored.

2 Introduction

Aeroacoustics research is conducted with theoretical analysis, experimental observations, and

computation. Computation is being used more often because of improvements in numerical

methods and computer resources, and is increasingly valuable as a source of insight and

as a tool for design [3, 18, 23, 28]. Artificial boundary treatments are a critical element

for Computational AeroAcoustics (CAA), since they are intended to be used to restrict a

numerical domain whithout causing unacceptable error in the numerical solution. A poor

boundary treatment can produce errors in the numerical simulation that grow to the size of

the solution itself, or larger, making the numerical solution of only qualitative or aesthetic

interest. A poor boundary treatment that might produce small errors could be very costly,

making an accurate numerical solution practically unobtainable. A good boundary treatment

is particularly important for simulating time dependant phenomena with detailed physics,

where experimants may be difficult because of cost or impossible because of scale and detail,

experiments such as for investigating jet noise. The recent reviews of the extensive work

done to date on this problem include [2, 5, 6, 13, 15, 16, 24, 27]. In this paper we are

concerned with Perfectly Matched Layer (PML) and Damping Layer (DL) treatments for

acoustic propagation by means of Linearized Euler Equations (LEE) in two space dimensions

with uniform and steady jet flows. A PML boundary treatment is designed for this system by

ensuring that the plane wave solutions for the LEE in the numerical domain perfectly match

the plane wave solutions for the adjacent PML domain along the interface between them.

Error Control With Perfectly Matched Layer or Damping Layer 

Treatments for Computational Aeroacoustics With Jet Flows 

John W. Goodrich 

National Aeronautics and Space Administration 

Glenn Research Center 

Cleveland, Ohio 44135 

NASA/TM—2010-216092 1
Control of the error introduced back into the numerical domain by the PML is attempted by

a combination of numerical domain filtering and interface conditions, and the PML domain

damping terms, damping rates and profiles, layer width, and outer boundary treatment.

In [11] we compared three PML treatments for this problem. Hagstrom [14] has used

transform methods to analytically formulate a PML boundary treatment (PML A) by operator approximation for this system, which has also been further analysed by Motamed [19]

and Fr ̈oling [4]. Since PML A has been developed with linear analysis for constant coefficient

problems, it is not readily apparent how it can be adapted to problems with nonuniform base

flows, or to fully nonlinear problems. Karni [17] has used eigensystem analysis to develop a

PML treatment (PML B) that use directional damping in an artificial boundary layer. Karni

[17] has observed that the damping rate need not be the same for the different eigenfunctions

of the hyperbolic system, but that the damping terms should not alter the eigenvectors of

the system with respect to the direction of damping. In [11] we showed that a Thompson like

characteristic analysis [25, 26] can be used to easily obtain a general damping layer formulation that does not alter the eigenvectors of the Hyperbolic system. The eigensystem analysis

that Karni has used could be applied locally, and is not restricted to constant coefficient

linear problems. The familiar DL or sponge zone treatment does not alter the eigenvectors

of the system, and the damping profile can be chosen to impose a smooth interface between

the numerical domain and the sponge zone or damping layer. A DL boundary treatment

has been recently discussed in [21], where the damping term σU~ is called a relaxation source

term, and the damping coefficient σ = 1/τ , where τ is a relaxation time. Bodony [1] has

also recently given an analysis of sponge zones for computational fluid dynamics. We have

called this familiar, simple, and easily generalized treatment PML C, even though it is not

properly a PML.

In [11] we showed that PML A and PML B were reasonably similar in performance, with

both producing disturbances that are one to three orders of magnitude smaller than those

produced by PML C for similar levels of damping. In this paper we consider both PML B

and C treatments by themselves, and we introduce and test a combination of PML B and

PML C that is intended to combine the best features of both. This combined PML treatment

will be called PML BC. PML BC has PML B terms predominate near the interface, and uses

a smooth transition away from the interface to an outer damping treatment that is purely

PML C. PML BC is essentially the same as a PML B with different damping profiles for

each directional eigenvector, and is intended to incorporate the best features of both separate

PML treatments, the accurate and nondistorting interface of PML B with the nonreflecting

omnidirectional damping of PML C near the outer PML boundary. Comparisons will be

made of the relative errors for combinations of damping profiles and layer widths.

3 PML/DL for the LEE with Uniform Base Flows

The nondimensionalised Linearized Euler Equations (LEE) in two space dimensions with

a uniform base flow (ub, vb)=(Mx, My) and a source can be written for the perturbation

NASA/TM—2010-216092 2
pressure and velocity U~ = (p, u, v)T as

∂U~

∂t + Ac

∂U~

∂x + Bc

∂U~

∂y = G, ~

where Ac and Bc are the constant matrices

Ac =





Mx 1 0

1 Mx 0

0 0 Mx



 , and Bc =





My 0 1

0 My 0

1 0 My



 .

In our numerical experiments we always use the pressure source G~ = (g, 0, 0)T , with

g(x, y, t)=0.01 sin[2πt] exp[−36(x2 + y2

)].

This compact radiating source propagates a time harmonic pressure disturbance with maximum values that are O[10−4]. These equations apply in the numerical domain

ΩN = {(x, y) ∈ [xL, xR] × [yB, yT ]},

where we are concerned with obtaining accurate results. The numerical domain is surrounded

with a combination of PML and DL domains, which we take to be Cartesian with linear

boundaries. The total combined domain is

Ω = {(x, y) ∈ [xL − wL, xR + wR] × [yB − wB, yT + wT ]},

where wL and wR are the PML/DL domain widths on the left and right, and wB and wT are

the widths on the bottom and top. The PML/DL or damping domain is just

ΩD = Ω − ΩN .

We will use the homogeneous LEE in the PML/DL domain ΩD, since the source term g is

virtually zero outside of the numerical domain ΩN . A source in ΩN near the interface with

ΩD could be accomodated merely by including the source in the damping domain equations.

For this constant coefficient case, a PML B treatment [17, 11] perpendicular to the x axis

is

∂U~

∂t + Ac

∂U~

∂x + Bc

∂U~

∂y + δσ(x)AcU~ = 0,

where σ(x) is the damping profile, δ the damping scale, and Ac the damping coefficient

matrix. Similarly, a PML B treatment perpendicular to the y axis is

∂U~

∂t + Ac

∂U~

∂x + Bc

∂U~

∂y + δσ(y)BcU~ = 0.

Note that since PML B is derived from a directional eigenvector analysis [17, 11], δ should

be positive for a layer on the right, and negative for a layer on the left, and similarly for up

and down. A PML C treatment can be written simply as

∂U~

∂t + Ac

∂U~

∂x + Bc

∂U~

∂y + δσU~ = 0,

NASA/TM—2010-216092 3
where σ would depend on the variable across the layer in either coordinate direction, with

δ ≥ 0. In any of the four corners, the PML/DL treatments from the two adjacent sides can

be combined by simply adding them. For PML B, this additive corner treatment is just

∂U~

∂t + Ac

∂U~

∂x + Bc

∂U~

∂y + δ (σ(x)Ac + σ(y)Bc) U~ = 0,

and for PML C it is

∂U~

∂t + Ac

∂U~

∂x + Bc

∂U~

∂y + δ (σ(x) + σ(y))U~ = 0.

A more general approach to corners for PML C is to simply make the damping profile

multidimensional, with

∂U~

∂t + Ac

∂U~

∂x + Bc

∂U~

∂y + δσ(x, y)U~ = 0.

The additive corner treatment is then just a simple example of a multidimensional damping

profile. In [11] we presented an approach for damping in a corner by blending two PML B

treatments from the sides. In this paper we will use additive corner damping for both PML B

and C, and the multidimensional corner damping profile for some cases with PML C. There

appears to be little difference in the errors produced by the various corner treatments, and

each of these three approaches produces no significant extra error than is produced by the

PML/DL treatments on the sides.

A Zero Boundary Data (ZBD) boundary condition is used on all outer boundaries of the

damping layers, with all data set equal to zero along the entire outer boundary. This is a

combined Dirichlet and Neumann condition since we keep and propagate spatial derivative

data at each grid point. Our rationale is that we are primarily concerned with the layer

treatment, and that if there were a good boundary condition, then the entire layer treatment

industry would be unnecessary. In addition, if any signal crosses a damping zone and reaches

an outer boundary with ZBD boundary condition, then the boundary condition will create a

reflected signal with approximately the same strength. This worst of all backscattering from

the outer boundary will show how well the layer treatment is performing. This is particularly

significant for PML B which damps data for outgoing eigenfunctions and amplifies data for

incoming eigenfunctions. Note that we do not use any filtering to eliminate this effect, but

just the damping in the formulation of the PML/DL layer treatment. The backscattering

can be reduced by using a better outer boundary treatment.

There are various ways to implement a PML or DL. One approach is to simply view

the damping terms as ordinary terms in the system, and treat them the same as any other.

For complex algorithms this could lead to significant extra effort. A simple, inexpensive

and general approach is described by Romenski, Titarev, and Toro in [21]. For the vector

variable U~ , consider the evolutionary system

∂U~

∂t + Π(U~ )=0,

NASA/TM—2010-216092 4
where Π is the propagator for U~ , with the associated damped system

∂U~

∂t + Π(U~ ) + δU~ = 0,

where δ is the damping rate. This is approximated in [21] by

U~ n+1 − U~ n

4t = −Π(U~ n) − δU~ n+1, or (1 + δ4t)U~ n+1 = U~ n − 4tΠ(U~ n).

Any explicit algorithm for an undamped evolutionary system can be written as

U~ n+1 = P(U~ n),

where the algorithm could be to any order in time, with any number of time substeps. The

damping in [21] can now be adopted as

(1 + δ4t)U~ n+1 = P(U~ n), or U~ n+1 = P(U~ n)

(1 + δ4t)

.

If the damping terms are not simple, say in the form

∂U~

∂t + Π(U~ ) + DU~ = 0,

where D is a matrix, such as in PML B, then we can write

(1 + 4tD)U~ n+1 = P(U~ n), or U~ n+1 = (1 + 4tD)

−1P(U~ n).

The inversion of the matrix damping term can be done locally at each grid point. These

formulations can be implemented as post processing step after each time step, or after each

time substep, or set of substeps. The damping terms for PML B and PML C have been

implemented with algorithmic treatments that are the same as any other terms with all of

the details required for full accurate time evolution, and also as this simple post processing

step, and there is little difference in the numerical domain solution, but there can be a large

difference in the required effort.

We have generally used four damping profiles. Let ξ = x or y, depending on whether or

not a damping term is for a layer perpendicular to the x or y axis, and let ξI be the location

in x or y of the interface between the numerical domain ΩN and the damping domain ΩD.

We assume that all profiles are set to zero inside the numerical domain. The first damping

profile was used for the constant coefficient cases in [11], with

σ0(ξ, w) = (ξ − ξI )2n

(w2 + (ξ − ξI )2)n ,

for all ξ inside the damping domain, where w is a scaling factor. Note that the order of the

zero of σ0 at ξI and the smoothness of the interface can be controlled by the choice of n.

Note also that σ0 asymptotes to one, instead of growing indefinitely like typical polynomial

NASA/TM—2010-216092 5
damping profiles. The damping profile is multiplied by the maximum damping rate or

amplitude δ. The damping amplitude δ combines with the layer width scale w to control

the total damping in the layer. The second damping profile is

σ1(ξ, w) = 

ξ − ξI

w

!2n 

2w − (ξ − ξI )

w

!2n

, for |ξI |≤|ξ|≤|ξI | + w,

and

σ1(ξ, w)=1, for |ξI | + w ≤ |ξ|.

Note that σ1(ξI , w) = 0, and σ1(ξI + w, w) = 1, and that the first 2n − 1 derivatives of σ1

are 0 at both points. The third damping profile is similar to the second, with

σ2(ξ, w) = 

ξ − ξI

w

!2n 

2w − (ξ − ξI )

w

!2n

, for |ξI |≤|ξ|≤|ξI | + 2w,

and

σ2(ξ, w)=0, otherwise.

Note that σ2(ξI , w)=0= σ2(ξI + 2w, w), and σ2(ξI + w, w) = 1, and that the first 2n − 1

derivatives of σ2 are 0 at all three points. Damping profile σ2 is used for PML B, and σ1 for

PML C, with n = 3 and additive corner treatments for both profiles. The fourth damping

profile is

σ3(ξ, w) = exp[−w2

/ξ2],

where ξ = |x − xI | on the right or left, ξ = |y − yI | on the top or bottom, and in the four

corners ξ =

q

(x − xI )2 + (y − yI )2. The multidimensional corner treatment is used with σ3

by design. Note that σ3 is infinitely smooth at the interface between the numerical domain

and the damping layers.

PML BC is a blend of PML B and PML C using a combination of terms with damping

profiles σ1 and σ2. For a PML BC perpendicular to the x axis on the right, the damping

terms are

D(x)U~ = δBσ2(x, wB)AcU~ + δCσ1(x, wC)IU, ~

where I is the identity matrix. The first term is just PML B for xI ≤ x ≤ xI + 2w, and

the second term is PML C. The use of PML B near the interface is intended to ensure that

lower errors are introduced by the PML, and the final persisting use of PML C is intended to

provide omnidirectional damping with no amplification of errors back towards the numerical

domain. A slightly modified form of PML BC delays the introduction of the third term for

a distance wg from the interface, and can be written as

D(x)U~ = δBσ2(x, wB)AcU~ + δCσ1(x − wg, wC)U. ~

This is the form that we prefer, with appropriate modifications for right and left, and up

and down PMLs. Note that there are five parameters, the two damping scales, δB and δC,

the two width scales wB and wC, and the ”gap” scaling wg before the third pure PML C

term becomes active. This PML BC uses a PML B in the layer or zone near the interface

for accuracy, with greater accuracy from a smaller damping coefficient δB, and PML C away

from the interface to prevent back scattering from the outer boundary, with greater damping

from a larger damping coefficient δC. The interplay of these five parameters will begin to be

explored in the numerical experiments reported below.

NASA/TM—2010-216092 6
4 Numerical Algorithm

All of the experiments that are reported in this paper are conducted with the Hermite/CauchyKovalevsky/Taylor method (see [7], [8], [9], [10], [12], [11]). This method uses two staggered

uniform grids, offset by half a grid spacing in both x and y. The first or basic grid is used

for initial data and for the numerical solution at each full time step. The second or offset

grid is used for the numerical solution at each half time step. The time evolution through

a half time step from one grid to the next starts with multidimensional Hermite spatial

interpolation of the local data on the current grid stencil, with a tensor interpolant for a

rectangular grid. The half time step then proceeds by using the governing equations for

Cauchy-Kovalevsky recursion to produce all of the required time derivatives from the spatial

derivatives produced by the interpolant. The half time step ends by propagating through a

half time step with a Taylor series. This approach is not simply an algorithm, but is actually

a method for developing or specifying numerical algorithms. The recursion routines are dependant upon the governing equations, but they can be written independantly from the order

of the method. Various methods can be realized for a particular system by simply swapping

the interpolation routines, using different stencils and data. This approach has had various

names, but in this paper, we will call it the Hermite/Cauchy-Kovalevsky/Taylor method,

or the HCKT method. In this paper we use the c2o1 method, which is a HCKT method

that uses a four point 2 × 2 grid cell, and interpolates with a bicubic Hermite interpolant at

the cell center, with data for {f,fx, fy, fxy} at each grid point for any evolving variable f.

This is a tensor interpolant which simultaneously and consistently calculates all derivatives

in the local bicubic spatial expansion. The spatial interpolant is used as a local initial data

surface, and all required time derivatives are computed from the local data surface by means

of a Cauchy-Kovalevsky recursion. The new solution data is obtained at the cell center

with a Taylor series in time. The time expansion is to fourth order terms, but because of

the bicubic tensor interpolant, the overall order is third in space and time, since a fourth

order time expansion requires fourth order spatial data which the cubic interpolant does not

provide. Note that the local time evolution is exact up to the order of the method. For two

dimensional hyperbolic systems, the local data is correctly propagated along characteristic

surfaces up to the order of the method. In this sense, the HCKT methods can be said to

correctly extend the method of characteristics to multidimensional hyperbolic systems. By

using the correct Cauchy-Kovalevsky recursion, this has been done for nonlinear problems

as well, in particular for convectively dominated compressible Navier-Stokes equations.

5 Results for the LEE with a Uniform Flow

The numerical experiments in this section are for the LEE with (pb, ub, vb) = (1, 0.4, 0) as

the base flow, and with the pressure source

∂U~

∂t + Ac

∂U~

∂x + Bc

∂U~

∂y = G, ~

where G~ = (g, 0, 0)T = (0.01 sin[2πt] exp[−36(x2 + y2)], 0, 0)T . The propagated variables are

all initially set to zero. All of the results are with the c2o1 HCKT algorithm, which is third

NASA/TM—2010-216092 7
order accurate in space and time. The algorithm is run with a Taylor series time expansion

that has fourth order terms, but the third order Hermite interpolant restricts the overall

algorithm to third order. The data that we report is

kp − pRk∞,ΩN = max{|p(x, y, tF ) − pR(x, y, tF )| : (x, y) ∈ ΩN },

where tF is the final time for the simulation being considered, and where pR is a reference

solution produced with no damping of any kind on a Big Enough Domain (BED) so that

the solution in ΩN is undisturbed at tF by any reflection back from the domain boundary.

Table 1: kp − pRk∞,ΩN at t = 15, from PML B with σ1.

δB wB = 1 wB = 3 wB = 5 wB = 7

10.000 4.2233 × 10−7 1.5443 × 10−7 3.7913 × 10−8 8.7781 × 10−9

1.000 8.8739 × 10−9 2.6081 × 10−9 4.5782 × 10−10 1.5450 × 10−10

0.100 7.2401 × 10−10 1.0766 × 10−10 2.1359 × 10−11 8.5276 × 10−12

0.010 7.3263 × 10−11 9.2124 × 10−12 2.4369 × 10−12 7.8338 × 10−13

0.001 7.3414 × 10−12 9.0569 × 10−13 2.4700 × 10−13 7.7644 × 10−14

Note that kpRk∞,ΩN = O[10−4].

Table 2: kp − pRk∞,ΩN at t = 15, for PML C with σ1.

δC wC = 1 wC = 3 wC = 5 wC = 7 wC = 9

10.000 2.5833 × 10−5 1.6333 × 10−6 8.2610 × 10−7 6.9061 × 10−7 2.9978 × 10−7

1.000 2.9492 × 10−6 5.1488 × 10−7 3.5265 × 10−7 1.2282 × 10−7 6.4461 × 10−8

0.100 2.4761 × 10−7 1.3017 × 10−7 4.9640 × 10−8 2.1429 × 10−8 8.7373 × 10−9

0.010 3.3743 × 10−8 1.3779 × 10−8 5.7895 × 10−9 2.3194 × 10−9 9.0230 × 10−10

0.001 3.4936 × 10−9 1.3810 × 10−9 5.8907 × 10−10 2.3381 × 10−10 9.0522 × 10−11

Note that kpRk∞,ΩN = O[10−4].

The first series of computations are with PML B and PML C by themselves, and are

presented in Tables 1 and 2. These calculations are for the numerical domain

ΩN = [−3, 7] × [−5, 5], with tF = 15,

with algorithm parameters 4x = 1/24 = 4y, and 4t = 1/48. The damping domain layer

widths are wL = wR = wB = wT = 10, the σ1 damping profile and additive corner treatments

are used for both PML B and PML C, and the damping terms are fully propagated. Note

that the propagation speeds of the wave fronts are 1.4 in the +x direction, 0.6 in the −x

direction, and 1 in the ±y directions, and consequently, that the propagating signal never

reaches the outer boundary of ΩD by tF = 15. The data in Tables 1 and 2 clearly shows that

errors introduced by a PML/DL can be controlled, even with a PML/DL on all sides and

corners of the numerical domain. The three controls that effect the results in Tables 1 and 2

are the type of layer that is used, the damping amplitude d, and the damping width scale w.

For every combination of d and w, PML B produces errors that are between two and three

NASA/TM—2010-216092 8
orders of magnitude smaller than those produced by PML C. The effect of both the damping

amplitude and width is linear for both PML B and PML C. The solution scale is O[10−4],

so that the relative error is from O[10−1], or ten percent, to O[10−10]. Any error that is

relatively less than O[10−3] is invisible at the scale of the solution, and this can be readily

achieved. It should be noted that if the damping layer widths wB or wC are greater than 5,

then the effect of the maximum damping profile amplitude is not returned to the numerical

domain ΩN where kp − pRk∞,ΩN is computed. The damping profile can cause an error by

its initial slope, with the layer scale and profile type both playing a role, by the maximum

damping amplitude, and possibly by higher derivatives from the profile shape. The smallest

errors in Tables 1 and 2 are produced by combinations of small damping amplitude and

large damping width scales. The total damping is a product of the amplitude and width,

so that a small damping amplitude needs a corespondingly wide damping layer to provide a

given total damping. The very smallest errors are produced by the very smallest damping

amplitudes, and they provide very little damping except over very wide damping layers.

Table 3: kp − pRk∞,ΩN at t = 25, for PML BC, with σ2 and wG = 3.

wC wB = 1 wB = 3 wB = 5 wB = 7

1 1.2269 × 10−6 1.2225 × 10−6 1.2278 × 10−6 1.2290 × 10−6

3 4.1117 × 10−7 4.1009 × 10−7 4.0892 × 10−7 4.1163 × 10−7

5 1.4603 × 10−7 1.4482 × 10−7 1.4424 × 10−7 1.4563 × 10−7

7 4.4231 × 10−8 4.2385 × 10−8 4.2241 × 10−8 4.3100 × 10−8

9 2.0671 × 10−8 1.9703 × 10−8 1.9697 × 10−8 1.9843 × 10−8

Note that kpRk∞,ΩN = O[10−4].

Table 4: kp − pRk∞,ΩN at t = 25, for PML BC, with σ1 and wG = 5.

wC wB = 1 wB = 3 wB = 5 wB = 7

1 6.8972 × 10−7 6.8669 × 10−7 6.8458 × 10−7 6.8825 × 10−7

3 7.2773 × 10−8 7.1685 × 10−8 7.1769 × 10−8 7.2137 × 10−8

5 1.6277 × 10−8 1.6506 × 10−8 1.6558 × 10−8 1.6707 × 10−8

7 5.4811 × 10−9 4.0522 × 10−9 4.0495 × 10−9 4.2073 × 10−9

9 5.0071 × 10−9 1.5652 × 10−9 1.2252 × 10−9 1.3843 × 10−9

Note that kpRk∞,ΩN = O[10−4].

The second series of computations uses the combined PML BC damping layer, with

results presented in Tables 3 and 4. These calculations are for the numerical domain

ΩN = [−6, 14] × [−10, 10], with tF = 25.

The algorithm parameters are 4x = 1/24 = 4y, and 4t = 1/48, the damping domain

widths are wL = wR = wB = wT = 10, the damping terms are fully propagated, and an

additive corner treatment is used. The propagating signal does not return from the outer

boundary of the damping domain by the final simulation time tF = 25. Recall that σ2 is a

pulse, so that PML B is only active close to the boundary, and that PML C applies only at

NASA/TM—2010-216092 9
a distance greater than wG from the damping domain interface. The damping amplitudes

δB = 1 for the PML B terms and δC = 1 for the PML C terms have been chosen as moderate

but relatively effective damping amplitudes. Note that the data in both tables shows that

the PML BC combination is insensitive to the length scale of the PML B terms. Note also

that the PML C terms decrease the effectivness of the PML B terms, and that the PML BC

errors are from one to three orders of magnitude larger than the comparable errors for PML

B by itself. Comparing the results in Table 3 with Table 2 shows that these PML BC results

are slightly better than the PML C results by a factor of about three, and that the PML

BC results are worse than the PML B results by from one to three orders of magnitude.

This suggests that the PML BC error in these computations is dominated by the effect of

the PML C terms. Comparing the results in Table 4 with Table 1 shows that these PML

BC results are worse than the PML B results by about an order of magnitude, and that the

PML BC results are better than the PML C results by about an order of magnitude. The

best performance of PML BC has PML C fully engaged at a distance of wg + wC = 12 or 14

from the interface. If PML C had been used by itself with wg = 0 and wC = 12 or 14, then

the relative advantage of PML BC would have been significantly reduced. The PML BC

combination can provide some damping layer performance improvement over a pure PML

C, but it appears from these initial results that this combination will work well only with a

relatively large gap before the PML C terms take effect. The best strategy for employing the

PML BC combination seems to be to reduce the solution amplitude by about three orders

of magnitude with the PML B terms, and then to use the PML C terms only after most or

all of this initial damping has been accomplished. The PML BC combination will require

wide damping layers to be effective.

Table 5: kp − pRk∞,ΩN at t = 35, for PML C with PML width 15 and σ3.

δC wC = 1 wC = 5 wC = 9 wC = 13

10.00 8.6856 × 10−6 4.4231 × 10−7 2.6493 × 10−7 1.6355 × 10−7

1.000 5.4827 × 10−7 1.4929 × 10−7 7.8054 × 10−8 3.2457 × 10−8

0.100 6.9008 × 10−8 3.0104 × 10−8 1.3658 × 10−8 6.6209 × 10−9

0.010 1.0729 × 10−8 4.0779 × 10−9 1.7786 × 10−9 7.3515 × 10−10

0.001 1.1433 × 10−9 4.2478 × 10−10 1.8288 × 10−10 7.4310 × 10−11

Note that kpRk∞,ΩN = O[10−4].

The third series of computations is with PML C by itself, and is presented in Table 5.

These computations use the σ3 damping profile, the multidimensional corner treatment, and

the Romenski, Titarev, and Toro [21] implementation of the damping terms. The numerical

domain is

ΩN = [−6, 14] × [−10, 10], with tF = 35,

and the damping domain widths are wL = wR = wB = wT = 15. The algorithm parameters

are 4x = 1/12 = 4y, and 4t = 1/24. The Damping zone is large enough so that the

propagated signal does not return by tF = 35 from the outer boundary of the damping

domain, so that the outer boundary does not effect the solution in the numerical domain.

NASA/TM—2010-216092 10
The data in Table 5 is comparable to the results in Tables 1 and 2, in the sense that

all three tables show the effect of the damping treatment without any reflection from the

outer boundary of the damping layer. Note that the damping profile scales in Table 2 are

wC = 1, 2, 3, 5, and 7, but in Table 5 they are wC = 1, 5, 9, and 13. It is immediately clear that

the results in Table 5 are only slightly different from the comparable results with the same

damping profile scale for PML C in Table 2, with errors that are smaller by no more than

a factor of two. The slightly better results in Table 5 could possibly be due to the greater

decay of the wave front from the initial impulsive start because of the longer simulation

time. On the other hand, for any particular damping amplitude δC, damping profile σ3

increases about three times more gradually than σ1, so that the σ3 profile scale is effectively

three times larger than wC. This more gradual increase in the damping profile will tend to

produce somewhat smaller errors in the numerical domain. Nonetheless, the similarity of

the results in Tables 2 and 5 shows that the damping profile, grid resolution, domain size

and simulation time, corner treatment, and implementation of the damping terms all have

no significant effect upon the accuracy of the damping treatment. Note in this regard, that

all of the damping profiles have been chosen so that, at the interface between the numerical

and damping domains, all of the spatial derivatives of the profiles are zero up to at least the

order of the numerical method. In particular, σ3 is C∞ with every spatial derivative equal to

zero at the interface, so that the interface between the numerical domain and the damping

layer provides a smooth transition between the governing systems in each domain. One of

the implications of the data that has been presented so far is that relatively small errors

are introduced into the numerical domain solution if the transition from the numerical to

the damping zone is gradual as well as smooth. The data in Table 5 also serves to validate

the algorithm mix that is used here, and to show that the algorithm mix and simulation

parameter selection is representative of the other results that have been discussed.

The fourth series of computations is with PML C by itself, and is presented in Tables 6,

7, 8 and 9. These computations use the same algorithm mix that produced Table 5, with

the same damping profile, multidimensional corner treatment, and implementation of the

damping terms. The numerical domain is the same,

ΩN = [−6, 14] × [−10, 10],

with the same algorithm parameters, 4x = 1/12 = 4y and 4t = 1/24. However, the final

simulation times and the damping layer widths are all different. Table 6 is for final simulation

time tF = 55 with damping domain widths wL = wR = wB = wT = 15, Table 7 is for tF = 65

with wL = wR = wB = wT = 20, Table 8 is for tF = 75 with wL = wR = wB = wT = 25, and

Table 9 is for tF = 85 with wL = wR = wB = wT = 30. In each case, the final simulation

time is large enough so that the transient effects from the impulsive start can pass out of the

numerical domain, and so that the propagated signal can reflect from the outer boundary

of the damping domain back through the entire numerical domain. The results in Tables

6-9 are designed to show the effect of reflection from the outer boundary as mediated by the

damping profile and the damping zone width. For almost all combination of δC and wC in

Tables 6-9, kp−pRk∞,ΩN decreases as the damping layer width increases, with exceptions that

appear to be minor variations around a limiting value. In these four tables, insufficient total

damping for a fixed δC is indicated if kp − pRk∞,ΩN increases with wC, which stretches the

dampng profile and therefore decreases the total damping in the layer. Similarly, insufficient

NASA/TM—2010-216092 11
Table 6: kp − pRk∞,ΩN at t = 55, for PML C with PML width 15 and σ3.

δC wC = 1 wC = 5 wC = 9 wC = 13

10.0 8.4234 × 10−6 9.9517 × 10−8 5.9050 × 10−8 9.2317 × 10−8

1.0 5.3772 × 10−7 3.8008 × 10−8 5.0535 × 10−8 6.8419 × 10−8

0.1 1.9056 × 10−7 6.2012 × 10−7 1.3412 × 10−6 2.1143 × 10−6

Note that kpRk∞,ΩN = O[10−4].

Table 7: kp − pRk∞,ΩN at t = 65, for PML C with PML width 20 and σ3.

δC wC = 1 wC = 5 wC = 9 wC = 13

10.0 8.3917 × 10−6 1.6884 × 10−8 2.2645 × 10−8 3.7814 × 10−8

1.0 5.4330 × 10−7 1.5659 × 10−8 2.6035 × 10−8 2.3734 × 10−8

0.1 7.2112 × 10−8 1.6864 × 10−7 4.3244 × 10−7 8.3129 × 10−7

Note that kpRk∞,ΩN = O[10−4].

Table 8: kp − pRk∞,ΩN at t = 75, for PML C with PML width 25 and σ3.

δC wC = 1 wC = 5 wC = 9 wC = 13

10.0 8.3890 × 10−6 2.0426 × 10−8 1.0227 × 10−8 8.0675 × 10−9

1.0 5.4134 × 10−7 9.5706 × 10−9 6.1355 × 10−9 1.1910 × 10−8

0.1 4.7353 × 10−8 4.7189 × 10−8 1.3228 × 10−7 2.9477 × 10−7

Note that kpRk∞,ΩN = O[10−4].

Table 9: kp − pRk∞,ΩN at t = 85, for PML C with PML width 30 and σ3.

δC wC = 1 wC = 5 wC = 9 wC = 13

10.0 8.3996 × 10−6 2.9749 × 10−8 1.2372 × 10−8 8.4349 × 10−9

1.0 5.4403 × 10−7 1.3398 × 10−8 6.3281 × 10−9 3.4250 × 10−9

0.1 4.6392 × 10−8 1.5651 × 10−8 4.0343 × 10−8 9.9061 × 10−8

Note that kpRk∞,ΩN = O[10−4].

total damping for a fixed wC is indicated if kp − pRk∞,ΩN increases with a decrease in δC,

which diminishes the damping profile, decreasing the total damping. Note from the data for

wC = 1 in Tables 6-9 that the errors for δC = 1 and 10 are the same for any PML width,

from 15 to 30, and that these errors are close to the comparable errors in Table 5 where

the outer boundary has no effect. We may conclude that a PML width of 15 is sufficient

to produce the best result possible with wC = 1 and δC = 1 or 10, or that the profile itself

creates an error level that is not diminished by further damping farther from the interface.

By the time the solution has propagated through a layer of width at least fifteen, it has

been damped sufficiently so that whatever error is reflected from the outer boundary is less

than the error produced by the damping layer with these profiles. Similarly, for wC = 1 and

δC = 0.1, a damping layer width of about 20 or slightly more is sufficient to produce the best

that is possible with this profile. Table 10 presents crude estimates of damping layer widths

that suffice to produce the best possible results with the damping profile and parameter

NASA/TM—2010-216092 12
combinations that were used to produce the error data in Tables 8-9. A total damping can

Table 10: Layer width for smallest errors with PML C and σ3.

δC wC = 1 wC = 5 wC = 9 wC = 13

10.00 15 20 25 25

1.000 15 25 25 > 30

0.100 20 30 > 30 > 30

be defined as

T D(wD, δC, wC) = δC

Z wD

0

σ3(ξ, wC)dξ.

If the results in Tables 6-9 are sorted into smaller and greater errors than the comparable

cases in Table 5, then

T D(wC, wD) ≥ 1.7 for smaller errors,

and

T D(wC, wD) ≤ 1.4 for greater errors.

This in principle permits at least the crude estimation of damping layer parameters that

would maintain a specified error bound. It must be remembered that relatively small errors

are produced if the transition from the numerical domain to the damping zone is gradual as

well as smooth, or with large damping profile scaling wC and small damping amplitudes δC.

6 DL for the LEE with a Nonuniform Base Flow

The LEE in two space dimensions with a nonuniform base flow requires two thermodynamic

variables, and we choose to use the disturbance pressure p and specific volume vs = 1/ρ,

where ρ is the density. In terms of the perturbation quantities

U~ = (p, u, v, vs)

T ,

the LEE can be written as: the pressure equation

∂p

∂t +  ̄u

∂p

∂x +  ̄v

∂p

∂y + γp ̄(

∂u

∂x +

∂v

∂y) + u

∂p ̄

∂x + v

∂p ̄

∂y + γp(

∂u ̄

∂x +

∂v ̄

∂y) = 0;

the velocity equations

∂u

∂t +  ̄u

∂u

∂x +  ̄v

∂u

∂y +

1

γM2

R

v ̄s

∂p

∂x + u

∂u ̄

∂x + v

∂u ̄

∂y +

1

γM2

R

vs

∂p ̄

∂x = 0,

and

∂v

∂t +  ̄u

∂v

∂x +  ̄v

∂v

∂y +

1

γM2

R

v ̄s

∂p

∂y + u

∂v ̄

∂x + v

∂v ̄

∂y +

1

γM2

R

vs

∂p ̄

∂y = 0;

NASA/TM—2010-216092 13
and the equation for the specific volume

∂vs

∂t +  ̄u

∂vs

∂x +  ̄v

∂vs

∂y − v ̄s(

∂u

∂x +

∂v

∂y) + u

∂v ̄s

∂x + v

∂v ̄s

∂y − vs(

∂u ̄

∂x +

∂v ̄

∂y) = 0;

where the base flow is ( ̄p, u,  ̄ v,  ̄ v ̄s)T , and where MR is a reference Mach number. The dimensionless variables for the base flow should satisfy

 ̄a2

M2

R =  ̄p/ρ ̄ =  ̄pv ̄s,

so that if we want the nondimensional speed of sound to be  ̄a = 1, then we must have

MR =

q

( ̄p/ρ ̄) = q( ̄pv ̄s).

For this paper we assume a constant base temperature, and choose

p ̄ = 1, and  ̄vs = 1, so that MR = 1.

In this case, we shall write the LEE with nonuniform base flows as

∂U~

∂t + Av

∂U~

∂x + Bv

∂U~

∂y + CvU~ = 0,

where Av, Bv, and Cv are

Av =





u γ  ̄ p ̄ 0 0

v ̄s/γ u ̄ 0 0

0 0 ̄u 0

0 −v ̄s 0  ̄u





, Bv =





v ̄ 0 γp ̄ 0

0  ̄v 0 0

v ̄s/γ 0  ̄v 0

0 0 −v ̄s v ̄



 ,

and

Cv =





γ

∂u ̄

∂x + ∂v ̄∂y ∂p ̄∂x

∂p ̄

∂y 0

0 ∂u ̄

∂x

∂u ̄

∂y

1

γ

∂p ̄

∂x

0 ∂v ̄

∂x

∂v ̄

∂y

1

γ

∂p ̄

∂y

0 ∂v ̄s

∂x

∂v ̄s

∂y −

∂u ̄

∂x + ∂v ̄∂y





,

and where U~ = (p, u, v, vs)T are the unsteady perturbation quantities. Note that these

coefficient matrices could be nonconstant in time as well as in space, and that a minor

modification would permit a nonuniform base temperature. We considered the LEE with

the pressure source

g(x, y, t)=0.01 sin[2πt] exp[−36(x2 + y2

)],

in the pressure equation. With this pressure source, the LEE with nonuniform flow can be

written as

∂U~

∂t + Av

∂U~

∂x + Bv

∂U~

∂y + CvU~ = G, ~

NASA/TM—2010-216092 14
where G~ = (g, 0, 0, 0)T .

In the case of a nonuniform base flow, a PML B treatment perpendicular to the x axis

will be adapted in the form

∂U~

∂t + Av

∂U~

∂x + Bv

∂U~

∂y + CvU~ + δσ(x)AvU~ = 0,

where σ(x) is a variable damping profile, δ is a damping scale, and the damping coefficient

matrix is just the x propagation matrix Av. A similar formulation can be written down for

a PML B treatment perpendicular to the y axis,

∂U~

∂t + Av

∂U~

∂x + Bv

∂U~

∂y + CvU~ + δσ(y)BvU~ = 0.

Note that since PML B is derived from a directional eigenvector analysis, δ should be positive

for a layer on the right, and negative for a layer on the left, and similarly for up and down.

A PML C treatment can be written simply as

∂U~

∂t + Av

∂U~

∂x + Bv

∂U~

∂y + CvU~ + δσU~ = 0,

where σ would depend on the variable across the layer in either coordinate direction, and

δ ≥ 0.

7 Results for the LEE with a Parallel Jet

The next set of experiments are for the LEE with a nonuniform base flow and a pressure

source. The pressure source is the same time harmonic oscillation of a narrow Gaussian

spatial profile that was used above,

g(x, y, t)=0.01 sin[2πt] exp[−36(x2 + y2

)].

The nonuniform base flow that we use here is a simple parallel jet flow, with

U~ = ( ̄p, u,  ̄ v,  ̄ v ̄s)

T = (1, 0.4+0.4 exp[−36y2

], 0, 1).

The numerical domain is

ΩN = {(x, y) ∈ [−3, 7] × [−5, 5]},

and the third order c2o1 HCKT algorithm is used, with 4x = 1/24 = 4y, and 4t = 1/60.

Solutions are computed for t ≤ 10. A PML C is used with damping profile σ1 and zone width

wR = wL = wB = wT = 5. The damping zone uses the PML C formulation for nonuniform

flows in the inflow and outflow damping zones

ΩDI = {(x, y) ∈ [−8, −3] × [−5, 5]} and ΩDO = {(x, y) ∈ [7, 12] × [−5, 5]},

NASA/TM—2010-216092 15
Table 10: kp − pRk∞,ΩN at t = 10, for PML C with σ1 and PML width 5.

δC wC = 1 wC = 3 wC = 5 wC = 7

10.000 1.7287 × 10−5 2.4439 × 10−6 4.1319 × 10−7 1.2125 × 10−7

1.000 2.4723 × 10−6 4.1704 × 10−7 6.4439 × 10−8 1.5573 × 10−8

0.100 3.2806 × 10−7 4.5465 × 10−8 7.2817 × 10−9 1.5973 × 10−9

0.010 3.5859 × 10−8 4.5992 × 10−9 7.3718 × 10−10 7.3718 × 10−10

0.001 3.6219 × 10−9 4.6049 × 10−10 7.3809 × 10−11 1.6018 × 10−11

Note that kpRk∞,ΩN = O[10−4].

and the PML C formulation for uniform flows in the top and bottom damping zones

ΩDT = {(x, y) ∈ [−8, 12] × [5, 10]} and ΩDB = {(x, y) ∈ [−8, 12] × [−10, −5]}.

The corners are treated additively and as if both side zones used the uniform flow PML C

formulation. The data that will be reported is the maximum absolute error in the numerical

domain kp − pRk∞,ΩN , where pR is a reference numerical solution at t = 10. The reference

pressure solution pR is O[10−4] in ΩN at tF = 10.

The data in Table 10 for PML C with a nonuniform flow is comparable to the data in

Table 2 for PML C with a uniform flow. Both sets of data are for the same numerical domain,

and the same spatial grid resolution. Because of the higher maximum velocity of the parallel

jet the time resolution has been increased from 4t = 1/48 to 4t = 1/60. Note that this jet

is superimposed upon the same uniform flow as was used for the experiments that produced

Table 2. Recall that the data in Table 2 is with a damping zone width wD = 10 and final

time tF = 15, while the data in Table 10 is with a damping zone width wD = 5 and final

time tF = 10. The same damping profile is used in both sets of calculations, with the same

space scales for the profiles, and the same set of damping amplitudes (except that wC = 9

has not been used for the parallel jet case). The damping terms have been implemented

with full time evolution, and not with the simpler postprocessing algorithm, just as for the

simulations that produced the data in Table 2. We note that almost all of the data for

kp − pRk∞,ΩN in Table 10 is smaller than the corresponding data in Table 2 by up to almost

an order of magnitude. As in Table 2, the errors in Table 10 range in absolute terms from

O[10−5] to O[10−11], or in relative terms from O[10−1] to O[10−7]. Note here as well that the

PML error decreases linearly with the damping amplitude, and linearly with an increase in

the damping profile width scale. As in the uniform flow case, the error induced by a PML C

damping zone with a nonuniform flow is controllable. Overall, as demonstrated by the data

in Table 10, the performance of PML C with a nonuniform flow is remarkably close to the

performance of PML C with a uniform flow.

8 Conclusions

Computations for the Linearized Euler Equations (LEE) have been reported for acoustic

propagation with a uniform flow and a parallel jet. Two damping layer treatments for the

NASA/TM—2010-216092 16
Linearized Euler Equations have been considered, Perfectly Matched Layer (PML) PML B

based upon directional damping [17, 11], and Damping Layer (DL) PML C [21, 11]. The

multilayer PML BC has been introduced, with a damping matrix that combines PML B

for accuracy near a PML interface with PML C for damping in an outer layer. A series of

numerical experiments have been conducted with all three boundary layers, primarily for

the LEE with a uniform base flow, but also with a jet flow. The intent of these experiments

has been to try and understand what are the critical parameters that effectively control

the errors produced by a PML/DL, and what is the practically achievable level of error

that can be obtained. Various combinations of numerical domain size, simulation time,

grid resolution, damping layer size, treatments for layer corners, implementation of damping

terms, damping profiles types, and damping profile parameter mix have been used in the

numerical experiments. The data from the numerical experiments shows that:

1. The errors introduced into a numerical solution by either a PML B or a PML C can

be controlled. Results from the numerical experiments show relative error reductions

by as little as O[10−1] to as much as O[10−10]. It seems clear that smaller errors can

be obtained with sufficient effort.

2. Control of the error from a PML/DL appears to be effected most strongly by the layer

width, and the amplitude and width scale of the damping profile. If the transition from

the numerical domain to the damping layer is sufficiently smooth, then effect of both

the damping amplitude and width is linear. Differences in the damping profile, grid

resolution, domain size and simulation time, corner treatment, and implementation

of the damping terms all have much less effect upon the accuracy of the damping

treatment. If the damping treatment offers insufficient total damping, then reflection

from the outer boundary of the damping layer can become the dominant source of the

error produced by the layer.

3. PML B produces approximately two orders of magnitude less disturbance in the numerical domain than PML C, but PML B amplifies distortions back towards the numerical

domain, while PML C dampens omnidirectionally.

4. The best strategy for employing the PML BC combination seems to be to reduce the

solution amplitude by about three orders of magnitude with the PML B terms, and

then to use the PML C terms only after most or all of this initial damping has been

accomplished. Because of this, the PML BC combination appears to require wide

damping layers to be effective.

5. Estimation of the total damping as a function of layer width, damping amplitude,

and damping spatial scale permit estimation of layer parameters needed to maintain

a specified error bound. Total damping is a product of the amplitude and width, so

that a small damping amplitude requires a corespondingly wide damping layer.

Relatively small errors are introduced into the numerical domain solution if the transition

from the numerical to the damping zone is gradual as well as smooth. Accurate damping

layers with O[105] relative errors are not difficult to achieve. Very accurate damping layers are

difficult to achieve and appear to require wide layers and substantial computational effort.

NASA/TM—2010-216092 17
The treatment of the outer boundary of the damping layer and its effect on controlling

the error produced by the layer as a whole has not been considered here. Every order of

magnitude decrease in the error produced by the outer boundary treatment will have a

significant effect on the error produced by the damping layer, or on the cost of producing

any specified error level.

References

[1] Daniel J. Bodony, “Analysis of sponge zones for computational fluid dynamics,” J.

Comp. Phys. 212 (2006) pp681-702.

[2] Tim Colonius, “Modeling Artificial Boundary Conditions for Compressible Flow,”

Annu. Rev. Fluid Mech. 36 (2004) pp315-345.

[3] Tim Colonius, Sanjiva K. Lele, “Computational aerocoustics: progress on nonlinear

problems of sound generation,” Prog. Aerospace Sciences 40 (2004) pp345-416.

[4] Maria Fr ̈oling, A Numerical Investigation of a Perfectly Matched Layer for Computational Aero Acoustics, TRITA-NA-E04041, Master’s Thesis, Royal Institute of Technology, Stockholm Sweden (2004).

[5] Dan Givoli, “Non-reflecting Boundary Conditions,” J. Comp. Phys. 94 (1991) pp1-29.

[6] Dan Givoli, “High-order local non-reflecting boundary conditions: a review,” Wave

Motion 39 (2004) pp319-326.

[7] J. W. Goodrich, “An Approach to the Development of Numerical Algorithms for first

Order Linear Hyperbolic Systems in Multiple Space Dimensions: The Constant Coefficient Case,” NASA TM 106928 (1995), and

[8] J. W. Goodrich, “Accuret Finite Difference Algorithms,” NASA TM 107377 (1995).

[9] J. W. Goodrich, “High Accuracy Finite Difference Algorithms for Computational Aeroacoustics,” AIAA 97-1584, The 3rd AIAA/CEAS Aeroacoustics Conference, (Atlanta,

GA, May, 1997).

[10] J. W. Goodrich, “A comparison of Numerical Methods for Computational Aeroacoustics,” AIAA 99-1943, The 5th AIAA/CEAS Aeroacoustics Conference, (Bellevue, WA,

May, 1999).

[11] J. W. Goodrich, “A comparison of Three PML Treatments for CAA (and CFD),”

AIAA 2008-2922, The 14th AIAA/CEAS Aeroacoustics Conference, (Vancouver, British

Columbia Canada, May, 2008).

[12] J. W. Goodrich, T. Hagstrom, J. Lorenz, “Hermite Methods for Hyperbolic InitialBoundary Value Problems,” Math. Comp. 75 (January, 2006) pp595-630.

[13] Thomas Hagstrom, “Radiation Boundary Conditions for the Numerical Simulation of

Waves,” Acta Numerica 8 (1999) pp47-106.

NASA/TM—2010-216092 18
[14] Thomas Hagstrom, A New Construction of Perfectly Matched Layers for Hyperbolic

Systems with Applications to the Linearized Euler Equations, in “Proceedings of the

Sixth International Conference on Mathematical and Numerical Aspects of Wave Propagation,” Finland, E. Heikkola, P. Neittaanmki, G. C. Cohen, P. Joly, eds., pp. 125-129,

Springer Verlag (2003).

[15] Fang Q. Hu, “Absorbing Boundary Conditions,” Int. J. Comp. Fluid Dyn. 18 (2004)

pp513-522.

[16] Moshe Israel, Steven A. Orszag, “Approximation of Radiation Boundary Conditions,”

J. Comp. Phys. 41 (1981) pp115-135.

[17] Smadar Karni, “Far-Field Filtering Operators for Suppression of Reflections from Artificial Boundaries,” SIAM J. Numer. Anal. 33 (1996) No.3, pp1014-1047.

[18] Sanjiva K. Lele, “Computational Aeroacoustics - A review,” AIAA-97-0018 (1997).

[19] M. Motamed, PML Methods for Aero Acoustics Computations, TRITA-NA-E03108,

Master’s Thesis, Royal Institute of Technology, Stockholm Sweden (2003).

[20] Peter G. Petropoilos, Li Zhao, and Andreas C. Cangellaris, “A Reflectionless Sponge

Layer Absorbing Boundary Condition for the Solution of Maxwell’s Equations with

High-Order Staggered Finite Difference Schemes,” J. Comp. Phys. 139 (1998) pp184-

208.

[21] E. Romenski, V. A. Titarev, and E. F. Toro, “Perfectly Matched Layers with High Rate

Damping for Hyperbolic Systems,” CFD J. 15 (2006) No. 2, pp240-246.

[22] James A. Sochacki, Absorbing Boundary Conditions Using Regions with Exponential

Decay, in “Mathematical and Numerical Aspects of Wave Propagation Phenomena,”

Proceedings of the First International Conference on Mathematical and Numerical Aspects of Wave Propagation, (Strasbourg, France, April 23-26, 1991). Finland, Gary

Cohen, Laurence Haloern, Patrick Joly, eds., pp. 293-302, SIAM (1991).

[23] Christopher K. W. Tam, “Computational Aeroacoustics: Issues and Methods,” AIAA97-1774 (1997).

[24] Christopher K. W. Tam, “Advances in Numerical Boundary Conditions for Computational Aeroacoustics,” AIAA J. 33 (1995) No. 10, pp1788-1796.

[25] Kevin W. Thompson, “Time Dependent Boundary Conditions for Hyperbolic Systems,”

J. Comp. Phys. 68 (1987) pp1-24.

[26] Kevin W. Thompson, “Time Dependent Boundary Conditions for Hyperbolic Systems,

II” J. Comp. Phys. 89 (1990) pp439-461.

[27] Semyon V. Tsynkov, “Numerical solution of problems on unbounded domains. A review.,” Applied Num. Math. 27 (1998) pp465-532.

[28] Meng Wang, Jonathan B. Freund, and Sanjiva K. Lele, “Computational Prediction of

Flow-Generated Sound,” Annu. Rev. Fluid Mech. 38 (2006) pp483-512.

NASA/TM—2010-216092 19
REPORT DOCUMENTATION PAGE Form Approved 

OMB No. 0704-0188 

The public reporting burden for this collection of information is estimated to average 1 hour per response, including the time for reviewing instructions, searching existing data sources, gathering and maintaining the 

data needed, and completing and reviewing the collection of information. Send comments regarding this burden estimate or any other aspect of this collection of information, including suggestions for reducing this 

burden, to Department of Defense, Washington Headquarters Services, Directorate for Information Operations and Reports (0704-0188), 1215 Jefferson Davis Highway, Suite 1204, Arlington, VA 22202-4302. 

Respondents should be aware that notwithstanding any other provision of law, no person shall be subject to any penalty for failing to comply with a collection of information if it does not display a currently valid OMB 

control number. 

PLEASE DO NOT RETURN YOUR FORM TO THE ABOVE ADDRESS. 

1. REPORT DATE (DD-MM-YYYY)

01-01-2010 

2. REPORT TYPE

Technical Memorandum 

3. DATES COVERED (From - To)

4. TITLE AND SUBTITLE

Error Control With Perfectly Matched Layer or Damping Layer Treatments for 

Computational Aeroacoustics With Jet Flows 

5a. CONTRACT NUMBER

5b. GRANT NUMBER

5c. PROGRAM ELEMENT NUMBER

6. AUTHOR(S)

Goodrich, John, W. 

5d. PROJECT NUMBER

5e. TASK NUMBER

5f. WORK UNIT NUMBER

WBS 984754.02.07.03.17.02 

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 

8. PERFORMING ORGANIZATION

 REPORT NUMBER

E-17137 

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

National Aeronautics and Space Administration 

Washington, DC 20546-0001 

10. SPONSORING/MONITOR'S

 ACRONYM(S)

NASA 

11. SPONSORING/MONITORING

 REPORT NUMBER

NASA/TM-2010-216092 

12. DISTRIBUTION/AVAILABILITY STATEMENT

Unclassified-Unlimited 

Subject Category: 71 

Available electronically at http://gltrs.grc.nasa.gov 

This publication is available from the NASA Center for AeroSpace Information, 443-757-5802

13. SUPPLEMENTARY NOTES

14. ABSTRACT

In this paper we show by means of numerical experiments that the error intruduced in a numerical domain because of a Perfectly Matched 

Layer or Damping Layer boundary treatment can be controlled. These experimental demonstrations are for acoustic propagation with the 

Linearized Euler Equations with both uniform and steady jet flows. The propagating signal is driven by a time harmonic pressure source. 

Combinations of Perfectly Matched and Damping Layers are used with different damping profiles. These layer and profile combinations 

allow the relative error introduced by a layer to be kept as small as desired, in principle. Tradeoffs between error and cost are explored.

15. SUBJECT TERMS

PML; CAA; Jet flows 

16. SECURITY CLASSIFICATION OF: 17. LIMITATION OF

 ABSTRACT

UU 

18. NUMBER

 OF 

 PAGES

25 

19a. NAME OF RESPONSIBLE PERSON

STI Help Desk (email:help@sti.nasa.gov) 

a. REPORT

U 

b. ABSTRACT

U 

c. THIS 

PAGE

U 

19b. TELEPHONE NUMBER (include area code)

443-757-5802 

Standard Form 298 (Rev. 8-98)

Prescribed by ANSI Std. Z39-18