1 ^ tiov, 985

RECEIVE0

ti

i

'^0

RAY W HERRICK LABORATOR IES

A Graduate Research Facility

of The School of Mechanical Engineering

f

d

Purdue University

West Lafayette, Indiana 47907

?VjheSA-CR-176295) PREDICTION Of SOUND FIELDS N86-13057

IN ICOUSTICIL CIVITIz;S USING THE BOUNDARY

ELEMENT !METHOD H.S. Zbesis (Purdue Univ. )

163 p HC A08/MF A01 CSCL 20A Unclas

G3/71 15955

M f
1 

TRANSMISSION OF SOUND THROUGH LIGHT AIRCRAFT

Purdue Research Foundation 520-1288-0353

i

PREDICTION OF SOUND FIELDS IN 

ACOUSTICAL CAVITIES USING THE

BOUNDARY ELEMENT METHOD

Sponsored by

NASA

0353-1 HL 85-26

Submitted by:

Carl R. Kipp, Graduate Research Assistant

R.J. Bernhard, Principal Investigator

Approved by.

Raymond Cohen, Director

Ray W. Herrick Laboratories

August, 1985

-_ii
E

t

PREDICTION OF SOUND FIELDS IN ACOUSTICAL

CAVITIES USING THE

BOUNDARY ELEMENT METHOD

i

a•

F Thesis

Submitted to the Faculty

of

Purdue University

by

Carl R. Kipp

In Partial Fulfillment of the

Requirements for the Degree

of

Master of Science in Mechanical Engineering

August 1985
ACKNOWLEDGEMENTS

could like t.c thank Dr. Robert J. Bernhard for

as my major professor and for providing advice

g ut this entire project. I would also like to thank

Arles M. Krousgrill, Dr. James F. Hamilton and Dr.

t Bolton for serving on the advisory committee.

s research was funded by the Langley Research

of the National Aeronautics and Space AdministrarSA) and I thank that organization for doing so.

Finally, but most importantly, 1 extend my most sincere heartfelt thanks to my fiancee, Ms. Paula Lewis.

Throughout this endeavor, she has provided me with her

.Love and support. During my times of bitter frustration,

she was always ready with words of encouragement. She made

it all worthwhile. I think I'll keep her.
1rF TABLE OF CONTENTS

Page

LIST OF FIGURES ........................................ v

LIST OF NOMENCLATURE .................................. ix

ABSTRACT ...............................................xii

CHAPTER 1 - INTRODUCTION ............................... 1

CHAPTER 2 - LITERATURE REVIEW .......................... 4

2.1 Finite Element Methods in Cavity Acoustics..... 5

2.1.1 Duct Acoustics............................. 6

2.1.2 Cavity Acoustics for Automotive

Applications ............................... 10

2.2 Boundary Element Acoustics ...................... 13

2.2.1 Investigations of Exterior Acoustical

Problems ................................... 16

2.2.2 Investigations of Interior Acoustical

Problems ................................... 28

2.3 Conclusions ..................................... 31

CHAPTER 3 - THEORY ..................................... 33

CHAPTER 4 - NUMERICAL IMPS .EMENTATION ................... 42

4.1 Element Definition .............................. 42

4.2 Equation Discretization ......................... 46

4.3 Primary Numerical Integration Procedure......... 50

4.4 Singularity Numerical Integration Procedure..... 55

4.5 Solution Prodecure .............................. 70

CHAPTER 5 - CASE STUDIES ............................... 72

5.1 Spherical Cavity Model .......................... 72

5.2 Spherical Mudei For Radiation .................. 98

5.3 Rectangular Cavity Model ........................102

5.4 Summary ......................................... 141
iv

CHAPTER 6 CONCLUSIONS ................................142

LIST OF REFERENCES .....................................145
v

LIST OF FIGURES

F igure Page

3.1 - Problem Description .............................. 34

4.1 - Local Isoparametric Element ...................... 43

4.2 - Global Isoparametric Element ..................... 44

4.3 - Integration Point Mesh ........................... 53 -3

4.4 - Vectors For A Flat Element ....................... 57

4.5 - Singularity Cases ................................ 59

4.6 - Circular Domain Mesh ............................. 65

4.7 - Transformed Circular Domain Mesh ................. 67

5.1 - Analytical Spherical Pressure Distribution

Fork-8 .......................................... 75

5.2 - Spherical Pressure Distribution For Pressure

Boundary Conditions At k-8 ....................... 77

5.3 - Spherical Pressure Distribution For Velocity

,• Boundary Conditions At k=8 ....................... 78

5.4 - Analytical Response For The Pressure At The

Wall Of A Pulsating Spherical Cavity ............. 82

5.5 - Numerical And Analytical Response For The

Pressure At The Wall Of A Pulsating Spherical

Cavity...•

. • • .

84

5.6 - Spherical Cavity Pressure Distribution For

Velocity Boundary Conditions At kr=2 ............. 86

5.7 - Spherical Cavity Pressure Distribution For

Velocity Boundary Conditions At kr=5.4........... 87

5.8 - Spherical Cavity Pressure Distribution For

Velocity Boundary Conditions At kr =7.1........... 88
vi

kEk

5.9 - Pressure Distribution For A Sphere With Radiation Impedance Boundary Conditions At

k-0.583 .......................................... 90

5.10 - Pressure Distribution For A Sphere With

Radiation Pressure Boundary Conditions At

k-0.583 ..................................... ... 92

5.11 - Pressure Distribution For A Sphere With

Radiaton Velocity Boundary Conditions At

r k-0.583 ......................................... 93

5.12 - Numerical Response For The Pressure At The

'M Wall Of A Rigid-Walled Spherical Cavity......... 95

5.13 - Pressure Distribution In A Rigid-Walled

Spherical Cavity With An Applied Source At

kr-4.5 .......................................... 96

5.14 - Analytical Pulsating Sphere Boundary Pressure Prediction For Radiation ...................100

5.15 - Numerical And Analytical Pulsating Sphere

Boundary Pressure Predictions For Radiation..... 101

5.16 - Analytical Pressure Distribution For A

Driven/Closed Cavity At k-1 .....................105

5.17 - Analytical Pressure Distribution For A

Driven/Closed Cavity At k-2 .....................106

5.18 - Analytical Pressure Distribution For A

Driven/Closed Cavity At k-3 .....................107

5.19 - Pressure Distribution For First Box Model At

k-1 .............................................109

5.20 - Pressure Distribution For First Box Model At

k-2 ................ ............................110

5.21 - Pressure Distribution For First Box Model At

k-3 .............................................111

5.22 - Pressure Distribution For Second Box Model

Atk-1 ..........................................112

5.23 - Pressure Distribution For Second Box Model

Atk-2 ..........................................113

5.24 - Pressure Distribution For Second Box Model

At k-3 ..........................................114

r
vii i

5.25 - Pressure Distribution For Third Box Model At

' k-1 .............................................115

5.26 - Pressure Distribution For Third Box Model At

k-2 .............................................116

5.27 - Pressure Distribution For Third Box Model At

k- 3 .............................................117

5.28 - Rectangular Cavity Pressure Response - X

Direction.......................................121

5.29 - Rectangular Cavity Pressure Response - Y

Direction .......................................122

5.30 - Rectangular Cavity Pressure Response - Z

Direction .......................................123

5.31 - Rectangular Cavity Pressure Distribution At

k-1.34 (XY Plane) ...............................126

5.32 - Rectangular Cavity Pressure Distribution At

k-1.34 (XZ Plane) ...............................127

5.33 - Rectangular Cavity Pressure Distribution At

k-1.34 (ZY Plane) ...............................128

5.34 - Rectangular Cavity Pressure Distribution At

k-1.65 (XY Plane) ...............................130

5.35 - Rectangular Cavity Pressure Distribution At

k-1.65 (XZ Plane) ...............................131

5.36 - Rectangular Cavity Pressure Distribution At

k-1.65 (ZY Plane) ...............................132

' 5.37 - Rectangular Cavity Pressure Distribution At

k-2.06 (XY Plane) ... ............................134

5.38 - Rectangular Cavity Pressure Distribution At

k-2.06 (XZ Plane) ...............................135

5.39 - Rectangular Cavity Pressure Distribution At

k-2.06 (ZY Plane) ...............................136

5.40 - Rectangular Cavity Pressure Distribution At

k-1 (XY Plane) ..................................138

5.41 - Recta► :-1121ar Cavity Pressure Distribution At

k-1 (XZ PI-et ) ..................................139

V.
I

r

viii

ity Pressure Distribution At

............ .....................140
Ix

LIST OF NOMENCLATURE

a,

U

Symbol Description

D Problem domain

X Vector locating a point in the domain

B Problem boundary

C Vector locating a point on the boundary

A

Normal vector at a point on the boundary

* Applied distributed source within the domain

t i Applied monopole source strength

X Vector locating an applied source

m Velocity potential

k Wave number

8 Dirac delta function

p Pressure boundary conition

t —

v Velocity boundary condition

j The complex number \1 -1

F Circular frequency

' P Fluid density

r id Distance between two points

t

o Source density distribution

c Singularity integration 'free term'
x

B Boundary excluding the region surrounding a

singularity point

p Fundamental pressure solution

v Fundamental velocity solution

T Solid domain angle

z Impedance boundary condition

`i

• A general variable

N i Lagrangian shape function

f Local element coordinate

71 Local element coordinate

N Number of boundary elements

M Number of nodes per element

K Number of applied sources

a il Numerical integration matrix term

d Applied sou*_ .: vector term

al Boundary condition vector term

g General function

NG Total number of integration points

( i Rectangular integration point coordinate

71 Rectangular integration point coordinate

w i Integration weighting factor

J Jacobian

r Radial integration point coordinate

8 Angular integration point coordinate

i j 0

Bessel's function of the first kind (0 order)

YO

Bessel's function of the second kind (0 order)

!3
xi

r0

A

B

a

b

Spherical surface radius

Boundary condition constant

Boundary condition constant

Boundary condition constant

Boundary condition constant

3
xii

ABSTRACT

Kipp, Carl R. MSME., Purdue University, August 1985.

Prediction of Sound Fields in Acoustical Cavities Using

the Boundary Element Method. Major Professor: R. J.

Bernhazd, School of Mechanical Engineering

A method is developed to predict sound fields in

acoustical cavities. The method is based on the indirect

boundary element method. An isoparametric quadratic boundary element is incorporated. Using this method, either

pressure, velocity and/or impedance boundary conditions

may be applied to a cavity. The capability to include

acoustic point sources within the cavity is implemented.

The method is applied to the prediction of sound

fields in spherical and rectangular cavities. All three

boundary condition types are verified. Cases having a

point source within the cavity domain are also studied.

_ Numerically-determined cavity pressure distributions and

responses are presented. The numerical results correlate

well with available analytical results.
I

CHAPTER 1

INTRODUCTION

Sound fields in cavities affect many aspects of

everyday life. It is well known that distribution ofsound in an interior work environment can markedly influence the performance of those working within that environment. Most modern modes of transportation place the traveler in some form of cavity. The perceived comfort level

of the traveler may again be influenced by the acoustical

characteristics of the transporting cavity. Consequently,

noise control engineers in the architectural and transportation industries become involved in optimizing acoustical

cavities so as to minimize the possible negative effects

on persons residing within those structures.

y
i

2

Computer aided design (CAD) procedures are very beneficial to individuals participating in noise control

activities. Using CAD procedures, the acoustic characteristics of a cavity can be determined analytically, usually with less time and expense than required for

experimentally-based techniques. Noise source identification can be accomplished with computerized modeling techniques. In addition, CAD procedures can readily be applied

in optimizing the acoustical characteristics of a cavity

after the noise sources have been identified.

The objective for the research presented here is the

t.evelopment of a numerical tool based on boundary element

theory to predict the sound field within an acoustical

cavity. The theory incorporates an indirect boundary element method utilizing isoparametric quadratic boundary

elements. In addition to the usual boundary conditions of

acoustic pressure and velocity, impedance boundary conditions are also formulated and implemented. Furthermore,

capability to model acoustical point sources within the

cavity domain is included. With all of the aforementioned

features, the computer program developed during this

research will provide a very useful tool for noise control

activities.

The content of this thesis is organized in the following manner. Chapter 2 contains an extensive literature
method for acoustics is presented in chapter 3. Numerical

implementation of the indirect boundary element theory is

presented ii

completed

ing from th

chapter 6.

t

1 AII,. .

3

review on numerical predictions of sound fields. The

{ theoretical background of the indirect boundary element
.t`

4

CHAPTER 2

LITERATURE REVIEW

Many researchers have been investigating techniques

for numerical prediction of acoustic fields. Consequently,

a significant body of technical and academic literature

' has been published. In contrast to numerical techniques

_ developed for specific problems, such as Green's function

pF.

methods and series solutions, the emphasis in this work

will be generalized techniques. Such techniques are versatile and powerful and will someday play an important role

in computer-aided noise control design. This chapter

presents a review of the literature applicable to the

numerical prediction of sound fields in acouPtical cavities using generalized techniques. Two methods are used

primarily, the finite element method and the boundary element method. The existing applications of each will be

discussed.

F

Y„ir

-
5

2.1 Finite Element Methods in Cavity A,g_Qustigs

The finite element method is the most widely-known

numerical tool available for engineering analysis of complex structures. A recent survey of finite element technology [1] contains a list of 155 textbooks and monographs

relating to finite element technology. This list does not

inci»3e the numerous journal articles relating to the

topic. It is obvious from the survey . that the finite element method (FEN) is now firmly rooted in various fields

of engineering analysis.

In acoustics, the finite element method is an established analytical tool. Most commercially-available finite element codes (e.g., NASTRAN, ANSYS, SAP IV) are

designed to be flexible, thereby providing the ability to

solve various classes of problems. With an understanding

F'

L of the specific code being utilized and the necessary

acoustical fundamentals, a finite element code may be used

to solve acoustical problems as shown by Bernhard [2]. 

The majority of finite element acoustical studies can

be classified into two primary categories based on the

application. These two categories are (1) duct acoustics

and (2) cavity acoustics for automotive applications.

3
:i

6

2.1.1 Dilct Acoustics

A common analytical method used to analyze acoustical duct systems is four-pole network theory. In fourpole network theory, a transmission matrix can be used to

define the relationships between the input and output

variables of acoustic pressure and particle velocity for

an acoustic device (e.g., a pipe, an expansion chamber, a

Helmholtz resonator). The sound attenuation or transmission loss of a device may be determined using the

transmission matrix. The transmission characteristics of

an acoustical system comprised of a series of well--defined

devices can be found by combining the respective transmis•..

sion matrices.

Young and Crocker [3] employed the finite element

method to determine the four-pole parameters of transmission matrices for expansion chambers. Since the initial

studies were for simple expansion chambers, twodimensional, four-node rectangular elements were used. 

Each element had twelve degrees-of-freedom : a pressure

and two pressure gradients, one in each coordinate direction, at each of the nodes. Utilizing fourth-order Hermitian shape functions, the transmission loss of a rigidwall expansion chamber was numerically evaluated and found

to correlate well with theoretical results. The chamber

was driven at the inlet with a harmonic velocity and
7

terminated at the outlet first with a rigid-wall boundary

condition and then with a pressure release boundary condition.

As an extension of this work, Young and Crocker studied complicated muffler systems having multiple chambers

(4]. Furthermore, the effects of non-rigid walls were

included. Numerically, a multiple chamber system is

assimilated by the cascade connection of the transmission

matrices for the individual components. For a more complex system, the finite element method was well-suited for

determining the four-pole parameters of the individual

components to be used in calculating the transmission loss

for the total system.

Astley and Eversman examined acoustical propagation

abcve the cut-off frequencies in ducts in terms of the

eiclenvalue problem [5]. The research included the effects

of steady mean flow or sheared flow within the duct. Capability was incorporated to study lined ducts by considering wall admittance. Two types of elements were considered. The primary element was a three-node element

with quadratic shape functions. The variables in the formulation were the perturbed velocities in the x and y

directions and the perturbed pressure. With the primary

element, extraneous eigenvalues were found to result at

relatively high frequencies. It was felt that the
8

I; extraneous eigenvalues were due to slope discontinuities

in the quadratic shape functions. To overcome this difficulty, two-node elements incorporating shape functions

having slope continuity were considered. The extraneous

eigenvalues were eliminated although the solution accuracy

was reduced at the higher eigenvalues. The finite element

method proved to be applicable to the solution of an

eigenvalue problem for ducts.

In addition to using acoustical finite elements,

Craggs studied the use of absorptive finite elements for

acoustical applications (6]. The implicit advantage of 

absorptive elements is that extended reaction of linings

can be included. Most techniques supply the lining admittance as a boundary condition thereby only allowing for

local reaction of the lining. Using the absorptive elements, the absorption coefficients for various resistivities of a constant thickness material were determined and

found to compare quite well with analytical results. The

acoustic and absorption finite elements were also coupled

to study a lined expansion chamber. The transmission loss

of the chamber as a function of lining resistivity was

determined. Results for the lined expansion chamber would

be very difficult to determine analytically. Thus, the

advantage of the absorption finite element is evident.

F-.

mom., arr sAMML'.E ^ 
9

ti The finite element method may also be used to examine

the sound field within a duct as was shown by Cederfeldt

[7). Towards that objective, Cederfeldt studied a rightangle bend, a lined rectangular duct and an expansion

chamber. The models were two-dimensional and comprised of

four-node quadrilateral acoustic elements. Contour maps

of the sound fields within the structures were presented.

A variation of the right-angle bend created by adding a

lining on the exterior corner provided quite a different

sound field as compared to the bend without the lining.

Even with the somewhat simplistic elements, dramatic

representations of the sound fields were obtained.

A more sophisticated element was used by Ling, Hamilton and Allen to determine the sound field in axisymmetric

ducts [8). The element was a two-dimensional, isoparametric, axisymmetric element with cubic polynomials as

the shape functions. With this element, the sound field in

a bottle-like duct was determined numerically and compared

with one-dimensional theoretical results and experimental 

measurements. Two types of lined ducts were also studied:

(1) a duct with segmented lining (i.e., hard wall and

impedance boundary conditions on different segments) and

(2) a duct with an exponential distribution of P.bsorbent

material along its length. Both of the above lining conditions produce a pressure distribution within the duct

which is difficult to obtain analytically but, as shown by
10

those researchers, can be obtained with the finite element

method.

Whereas Ling, Hamilton and Allen used an isoparametric element, Doyle and Faulkner applied a subparametric element for acoustic pressure analysis [9]. The

geometry variation within the element was linear while the

pressure variation was made to be cubic under the philosophy that the subparametric: element would be more computationally efficient than an isoparametric element. At each

grid point, the four independent variables of pressure and

the pressure gradient in each of the three orthogonal

? coordinate directions were used. An additional feature

highlighted in this research is the ability to assign different fluid densities and speeds of sound at any grid

point. With the method as described above, the eigenvalue

problem for a uniform hard-walled duct was solved to find

_ both the natural frequencies and the mode shapes of the

' duct. A second type of problem was analyzed where a harmonic velocity was assigned to one wall of the duct. The

pressures in the duct were determined for a hand-walled

case and cases with acoustic impedances on one or more

walls_

2.1.2 Cavity Acoustics for Automotive ApplicAUQp&

Cavity acoustics problems for automotive applications are generally more complex than those for duct

3 

T
i

s

11

acoustics problems. The geometry of the structures are not

simple rectangular or cylindrical shapes. In addition, the

interior of the cavity may contain obstacles affecting the

sound field. Forcing functions can not normally be defined

with ease when forced response problems are considered.

For these reasons, the research presented below has been

separated from the duct acoustics research presented in

the previous section.

Craggs studied the sound field in a passenger car

cavity when the cavity is coupled to the engine compartanent through a flexible boundary [10). The boundary

motion was formulated in terms of the normal modes of the

enclosure. The governing equations for a single enclosure

with a flexible boundary and the structural coupling

between two enclosures were developed. An example of a

passenger car compartment coupled through a flexible boundary with the engine enclosure was then studied using the

numerical formulation. The model could be used to study

the sound field induced within the cavity due to an excitation within the' engine compartment. From the results,

some generalizations were made regarding the situations in

which the passenger cavity either was or was not greatly

influenced by the excitation in the engine compartment.

Sung studied finite element applications to the

acoustics of a passenger car compartment and an engine
12

combustion chamber (11]. The NASTRAN finite element code

was used for both cases. Three-dimensional models of the

passenger compartment were constructed for two cases: (1)

a compartment with a bench seat and (2) a compartment

with bucket seats. The resonant frequencies of the cavities were determined and the corresponding mode shapes

e

presented with an acoustic isobar representation. Inclusion of the seat in the model provided a more realistic

situation than a model without the seat. Pox the combustion chamber model, a piston head with an asymmetric bowl

was included in the chamber model. Again, the resonant

cavity frequencies and mode shapes were determined from

the finite element model. These resonant frequencies were

compared with experimentally determined resonant frequencies and found to have good agreement. The final observation was that the acoustic finite element models were successful in predicting cavity boom frequencies in passenger

compartments and knock-induced frequencies in a combustion

chamber.

In conjunction with their work in duct acoustics,

Doyle and Faulkner extended the research by studying the

sound field in an automobile interior (9]. The first analyses performed were similar to those done by Sung (11] in

that the resonant frequencies and mode shapes were quantified. However, more interesting results were obtained for

a forced response problem. Experimentally determined

-'kf-
13

pressures at a sun roof opening were supplied to the finite element model as the forcing function. The acoustic

field produced by the forced model was compared with a

discr*te number of points at which experimental data had

been collected. Although the numerical pressures were

alightiy less than the experimental pressures, there was

good correspondence between the two sets of data.

The above-cited references exemplify the types of

acoustic cavity problems which may be studied with finits

element techniques. Moreover, they indicate the popularity of finite element acoustics. The following section

presents acoustical applications of the Boundary Element

Method. Although not yet as popular as finites element

techniques, boundary element techniques share many of the

attractive features of the finite element method and are

probably better adapted to certain problems.

V 2.2 Boundary Element Acoustics

There are two basic classifications of Boundary Element Methods: (1) Direct 3oundary Element Methods (DBEM)

and (2) Indirect Boundary Element Methods (IBEM). The distinction between the two methods must be made at this

point for clarification of the following discussion. For

they current project, an IBEM formulation is utilized as

will be shown in Chapter 3.
14

The difference between the two classif icati( , ns is

evident in the integral equations. For the DBEM, "the unknown functions appearing in the integral equations are the

actual physical variables of the problem" [12]. The physical variables for an acoustics formulation are the acovatic pressure and the particle velocity. A general form of

the DBEM integral equation is

PM - I[p*(X„)v(X) - p(X)v *(X,F)]dB

B

+ ;p*(x,f)V(x)dD (2.1)

D

_i where p(f) is the pressure at any point, p(x) is the pressure at a boundary point, v(x) is the normal velocity a, a

- 1

boundary point, p (x,f) and v (x,E) are the respective

fundamental pressure and velocity functione from a boundary point to any point and V(x) represents applied source

strengths over the domain of the problem.

The formulation for an IBEM is significantly different from that for a DBEM. For the IBEM,

"The integral equations are expressed entirely in

terms of a unit singular solution of the original

differential equations distributed over the boundaries of the region of interest. The density funcLi ins themselves have no specific physical significance but once they have been obtained from a numerical solution of the integral equations the values of

the solution parameters anywhere within the body can

be calculated from them by simple integration

processes."[12]

Or, in other words, the hosndary is replaced by a distribution of sources which reproduce the specified boundary
•

t

15

solution and allow computation of the physical variables

anywhere in the problem domain or boundary. Separate

equations are required for portions of the boundary

depending upon the specified boundary conditions on the

particular portion of the boundary. The general form of

f the IBEM integral equation for a pressure boundary condition is

PM - fa (x)P * ( x , f )dB + fPx (x, t)f(x)dD (2.2a)

B D

whereas the general equation for a velocity boundary condition is

V(f) - *CO(x) + fc(x)v * (x,E)dB + fv *(x,E)f(x)dD(2.2b)

B D

where o(x) represents the fictitious source density functions at the boundary points and c represents an Integra-

` tion constant due to integral singularities which will be

further explained in chapter 3.

Although the details of the two formulations are not

critical at this point, an understanding of the existence

of differing formulations is important. Some of the following reviewed works use a DBEM while others use an IBEM.

This point will be emphasized within each review.

7

t
16

2.2.1 Investigations of Exterior Acoustical Prnhlramn

Most of the boundary element research heretofore

published has dealt with the acoustic field of a domain

exterior to a vibrating or scattering structure. Primarily, this is a consequence of the capability of boundart' element methods to model domains extending to infinity, a task not easily accomplished with finite element

models. However, the integral formulations of interior

regions are nearly equivalent to those for exterior

regions. Hence, much of the development for exterior

regions may be extrapolated to developments for interior

regions.

Chen and Schweikert (13) were the earliest contemporary researchers to apply boundary element techniques to

acoustic radiation problems. A numerical implementation of

Huygen's principle, which is essentially a primal IBEM

formulation, was utilized whereby a distNr...bution of simple

sources was imagined over the bounding surface. The surface itself was considered to be an array of triangular

elements over which a uniform velocity with a constant

phase angle was prescribed on each element. Each element

was analogous to a rigid piston vibrating with a constant

harmonic behavior. After defining a model by its geometry

and boundary conditions, the boundary integral equations

were numerically evaluated to calculate the simple source

J s
17

i

distribution over the model boundary. The field pressures

were quantified from the source distribution on the boundary. Two examples were given. The first example consisted

3

of a piston in a rigid sphere for which the exterior field

' pressure radiation pattern was calculated. Likewise, similar exterior pressure patterns were determined for a stiffened cylinder immersed in water. The fluid-structure

interface effects were included for the second example.

A DBEM was incorporated into a study of acoustic

radiation from surfaces of revolution by Chertock (14].

The variables in the integral equations for a DEEM

represent the acoustic pressure and particle velocity. So

by specifying one type of boundary condition, the remain- 

ing unknown boundary condition may be determined with the

Helmholtz integral equations. Furthermore, the field pressures can be calculated based on the boundary parameters.

Chertock used this procedure to numerically evaluate the

surface pressures and field pressures from the velocity

boundary conditions for rigid-body vibration of a sphere, 

quadrapole vibration of a prolate spheroid (i.e., a

cigar-shaped surface) and quadrapole vibration of an

axisymme, hric surface which was unsymmetric about all axes

except tAe major axis. Surfaces of revolution were studied because theoretical results may be obtained for comparison with numerical results. Good correlation was found

`s
18

between the numerical and theoretical solutions for the

cases in question.

i

Copley [15] also applied a DBEM to the study of

acoustic radiation from axisymmetric radiators. However,

Copley (15) used an approach which differed from that of

Chertock (14]. Copley utilized the Interior Helmholtz

integral equation whereas Chertock utilized the Surface

Helmholtz integral equation. In the interior Helmholtz

integral approach, the field points of integration lie

` within the bourdary of the vibrating body and, in faces,

are located along the axis of symmetry for axismmetric

radiators. The corresponding field points lie on the

boundary in the surface Helmholtz integral approach.

Hence, the appropriate titles of each approach refer to

the location of the field integration points. Using the

Interior Helmholtz integral equation, the acoustic radiation from a spheroid and a finite cylinder with a capped

end were studied. The f arfield directivity pattern of the

finite cylinder was compared to results obtained by Williams, et.al . using a method based on expansion in spherical harmonics [16]. Generally, the results correlated

well with the exception of some discrepancy which was

thought to be due to differences in which the velocity

distributions were prescribed.
19

U Taking a more general point of view, Copley (17) also

studied failures which occur with integral formulations

for acoustic radiation problems. At certain frequencies,

the integral formulations are unable to represent the

steady-state harmonic radiation from a finite, smooth,

closed surface on which normal velocity is prescribed.

These failures occur at the interior Dirichlet eigenfrequencies for both the IBEM and the Surface Helmholtz

integral equation which is a DBEM. Copley [17] presented

proofs that the failures were due to nonexistence of

sources for the IBEM and as nonuniqueness of sources and

doublets for the Surface Helmholtz integral equation. The

failures are inherent in the formulations and do not

represent physical attributes of the problem. However,.

the Interior Helmholtz integral formulation used by Copley

[15] does not suffer from these failures.

Proofs similar to those of Copley [17] were presented

by Burton and Miller [18]. Moreover, suggestions to overcome the failures were made. One alternative consisted of

overdetermining the system of equations by using both the

Surface Helmholtz integral equations and its differentiated form. The second but more economical alternative

was to combine the two types of equations into a single

set of equations which would not overdetermine the system.

No computational examples were provided for these

approaches.

-s

.r

y._. 
20

A combined Helmholtz integral equation formulation

was developed by Schenck [19) to overcome the deficiencies

of the above mentioned formulations. The technique initially applyed the Surface Helmholtz integral equations

and then overdetermined the system of equations with the

Interior Helmholtz integral equations. The overdetermined

set of equations was solved by a least squares orthonormalizing procedure to evaluate the surface pressures. From

the surface pressures, the field pressures were calculated

with the Helmholtz integral equation. Essentially, the

described procedure is a DBEM as all the equations are of

the direct formulation type. A number of numerical examples were presented for the combined Helmholtz integral

equation formulation of which only three will be discussed. The first example is a uniformly vibrating sphere.

The problem illustrates the advantage of the combined

equation formulation. At a characteristic wave number, the

Surface Helmholtz integral equations fail to yield the

correct surface pressure. However, the combined Helmholtz

integral formulation was capable of correctly solving for

the surface pressure on the sphere by specifying only one

point within the sphere. A second example of a right circular cylinder with rigid ends on which velocity boundary

conditions were prescribed was studied. The resulting farfield pressure pattern was found to be consistent with

that obtained by Chertock (14) but inconsistent with that
obtained by Williams, et.al . [16]. As a variation on this

model, the velocity distribution was revised to reflect

the velocities that would occur at the boundary if a

source were at the center of the cylinder and the cylinder

was not physically present. The acoustic pressure radiation pattern was equivalent to a pattern for a simple

source as was expected. A rectangular parallelepiped was

the third model considered. Again, the velocity distribution due to a simple source was assigned on the boundary.

The resulting radiation pattern was the expected simple

source radiation pattern. For both of the latter cases,

the Surface Helmholtz integral equation without the interior points failed at the characteristic wave numbers.

The concept of overdetermination of the system equation4 was also addressed by Paiszczyk and Klosner [20].

The overdetermining equations were developed with the

Exterior Helmholtz integral equations unlike Schenck's

research which used the Interior Helmholtz integral equations for the overdetermination [19]. An iterative pro- 

' cedure is required in the method proposed by Piaszczyk and

Klosner [20]. At the surface, an approximate acoustic

fluid impedance was assumed. The Exterior Helmholtz

integral equations were used to solve for approximate

pressures at selected field points. The approximate field

pressures were then applied in conjunction with the Exterior and Surface Helmholtz integrals to solve for the
r

a

22

surface pressures. A least squares procedure was required

to calculate the surface pressures. The surface pressures

were then substituted into the overdetermined set of equations to solve for the field pressures which were considered to be the next iterative value of the field pressures. The entire iterative process continued until a convergence was established. For this exterior overdetermination scheme, the selection of the overdetermining points

is not critical due to its iterative nature. Thus, the

scheme may be applied without loss of generality to structures of arbitrary shape. A number of computational examples were presented displaying the versatility and accuracy of the exterior overdetermination method.

The research of Schenck [19] was used in combination

with finite element analysis to study acoustic radiation

from sonar transducers by Smith, Hunt and Barach [21].

Although the combined Helmholtz integral equation formulation was applied, only the Surface Helmholtz integral

equations were utilized in the sonar transducer study. A

specific acoustic impedance matrix was derived with the

integral equations. Definition of consistent mass and

stiffness matrices for the complex structures came from

the finite element analysis. The total structuralacoustical system of equations combined the mass and

stiffness matrices with the acoustic loading included as

surface loading forces. Finally, the pressure radiation

d
23

patterns were determined using the results of the

structural-acoustical analysis and the Helmholtz integral

equations to evaluate the pressures at points in the field

surrounding the transducer. Numerical results for a

piezoelectric sphere and a piezoelectric free-flooded

cylinder were compared with experimental results and found

to agree within 5 per cent. Some of the error was attributed to the piezoelectric nature of the material in the

transducers.

Engblom and Nelson [22) expanded upon the work of

Smith, Hunt and Barach [21] by the addition of two

features. The first additional feature involved allowing

quadratic variation of the acoustic variables and linear

variation of the geometrical variables over the surface of

an element. Prior to this work- a constant variation of

the acoustic variables had been assumed. In doing sc,

discontinuities of acoustic variables may arise at the

boundary of neighboring elements. A quadratic variation

ensures that the acoustic variables become continuous at

the element boundaries. The second feature added was a

coordinate transformation used to integrate around the

singularity which exists in the Helmholtz integral formulations. Since the elements were triangular, the coordj.-°

nate transformation involved a change from area coordinates to polar coordinates. In doing so, the singularity

was removed thus allowing the integration to be completed

6
i

24

without the singularity complications. The common examples

W! acoustic radiation from a sphere and a right circular

cylinder were numerically evaluated with good results.

Unlike the approach originally developed by Schenck

[19] to overcome the difficulty of nonuniquaness, Meyer,

Bell, Zinn and Stallybrass [23] implemented an approach

suggested by Burton and Miller [18]. That is, the system

of equations consisted of the Surface Helmholtz integral

equation and its differentiated form, both of which were

direct formulations. The system of equations was not an

overdetermined set but a combination ot the two types of

integral equations. A difficulty arose in that the differentiated form of the Surface Helmholtz integral equation contained a strong singularity which could not be

directly numerically integrated. Nevertheless, the authors

proved through rijorous mathematics and computational considerations that the integration was possible. For an

example of a piston set in a rigid sphere, the error for

the numerical farfield pressures remained less than 10 per

cent as compared with exact analytical solutions.

One further refinement of the DBEM for acoustic radiation problems was achieved by Seybert, Soenarko, Rizzo

and Shippy [24]. The authors provided further sophistication in the discretization process through the use of an

isoparametric element. With the isoparametric element,
25

both the geometry and the acoustic variables were int:erpo-

.t

lated with quadratic shape !unctions. Although not immediately apparent, the isoparametric element is somewhat more

refined then the element used by Engblom and Nelson [22]

which used a linear geometrical variation. Consequently,

efficiency and accuracy were increased as lose elements

were required to model curved or other irregular surfaces.

To illustrate the isoparametric element capability, the

examples of a pulsating sphere and an oscillating sphere

were considered and found to compare well with theoretical

results except at the characteristic wavenumbers. The

nonuniqueness of the formulation was not accounted for by

these authors.

The preceding reviews in this section represent the

advances which have occurred in boundary element technology for exterior domain acoustics, particularly for radiation problems. Many of the same authors have also studied

scattering problems. Thus far, the theoretical implications have been discussed. The remaining portion of this

section will review applications of the aforementioned

theory to somewhat more complicated problems.

A DBEM has been used by Seznec (25) to study the diffraction of sound around barriers. The barriers were any

general type or shape used for noise abatement (e.g., a

barrier around a roadway to hinder traffic noise). The
' 2 6

F

capability to include reflectivity and absorptivity

effects of both the ground and the barrier was achieved

through the application of impedance boundary conditions

on the surfaces. Prior to the work of 3eznec, researchers

had been unable to adequately account for the reflectivity

and absorptivity effects. The diffraction of sound around

barriers of various shapes were determined. For one barrier, that of a infinitely thin reflecting barrier, the

numerical solution was favorably compared with solutions

already existent in the literature. overall, the DEEM was

found to be effective as a design tool in studying the

diffraction of sound around barriers.

3adek and his coworkers developed a procedure to

C

predict the acoustic emission from a machine using the

^. design drawings for the machine (26]. A finite element

_

model of the machine structure was used to predict the

structural modes of vibration and this modal response to

typical operating loads for the machine. The modal

response yielded the surface velocities which were used as

boundary conditions for a DBEM formulation to predict the

acoustic emission. Two applications of the procedure were

presented. The first example was for a forging machine

structure under impulse loading characteristic of the

machining operation. Reflections of sound from the ground

_ could be and were included in determining the acoustic

emission. Comparison of the numerical results with

wr
27

experimental results showed that each exhibited the same

general trends and magnitudes although there were variac

tions at specific points. a , second application involved 

f

studying the acoustic emission from one component of a

R

hydraulic hammer [27]. The component considered to be the

most offending noise source was studied. The acoustic

emission analysis procedure was applied to the initial

design and possible redesign configurations. As a reault

of these analyses, one redesign configuraticei quantitatively produced lass acoustic emission than the initial

:i design and the other redesigns. The process followed in

this work emphasizes the manner in which a boundary element method may be used in the design process.

A method to predict acoustic intensity on the surface

of a vibrating body utilizing a DBEM formulation has been

developed by Benner and Koopmann [28,29]. Much like the

previously presented research, a finite element code was

used to define the dynamics of the structure. Since determination of acoustic intensity requires both the pressure

and velocity, both of these quantities were calculated

with the boundary element computer erode. In addition,

acoustic pressure and velocity, and hence acoustic; intensity, were calculated at any field point. The boundary

element acoustic intensity method was applied to the study

of crack propagation along the boundary of a flat plate

[30]. To simulate crack propagation, successive collinear

f

:i

t

3

d
degrees-of-freedom were released within the finite element

model. The acoustic surface intensity was calculated for

various modes of vibration as each additional degree-offreedom was released. A definite change in the intensity

distribution occurred as the crack propagated. This application highlighted the possible use of a boundary element

program as a tool in analyzing cracked or flawed structures for non-destructive testing. A second application

was to study the radiation characteristics of a slender

box-type structure representative of a box girder in a

rapid transit transportaticn system [31]. The intent was

to isolate high radiation efficiency modes so that damping

treatments could concentrate on :educing those modes_

Relatively high radiation efficiency modes were identified

using the boundary element formulation but analysis of the

structure with damping was not presented. This application

highlights one manner in which boundary element formulations-may be used for design modifications.

2.2.2 Investiaations of Interior Acoustical Problems

The preceding section and the vast majority of

literature relating to boundary element technology for

acoustical applications consider exterior domain problems.

The following discussion reviews literature applied to

acoustical problems in the interior domain.
if

i

k

s

29

Tai and Shaw [ 32] have used a DBEM formulation to

determine eigenvalues and eigenmodes of the Helmholtz

equation for arbitrary domains. The technique is applicable to both two- and three -dimensional domains. However,

the numerical example presented was for a two-dimensional

domain. Within an element, the acoustic variable was considered to be constant. The numerical example presented

consisted of a right triangle along whose hypotenuse the

eigenmodes were determined after the eigenvalues had been

found. The first two modes were found for three cases having varying magnitudes of the acute angles. This problem

wds Jo chosen to emphasize that the integral formulation

is capable of providing results when a separation of variables approach is not. The separation of variables technique can be used when the boundaries are along constant

coordinate lines thereby allowing the equations to

separate into individual coordinate directions. For the

case of a right isoceles triangle, which has a known

analytical solution. the numerical solution with the DBEM

formulation and analytical solution showed excellent

agreement.

Tanaka, Fujikawa, Abe and Utsuno developed a method

incorporating a DBEM formulation to calculate the transfer

matrices used in analyzing muffler systems [33]. This

method was found to yield good accuracy and better efficiency than the finite element approach applied by Young

w 
30

f and Crocker [3,4]. Boundary element theory was applied in

regions where it is difficult to apply plane-wave theory.

f Plane-wave theory was applied in regions where that theory

was still applicable. The transfer matrix for the entire

system was obtained by multiplying the matrices from each

component and used to calculate the insertion loss.

Impedance boundary condi^ions were included for some

cases. The method combining boundary element and planewave theory gave good results and also decreased both

modeling and computational time as compared to other

numerical techniques.

Sestieri, Del Vescovo and Lucibe : lo have onsidered

the problem of structural-acoustical coupling in cavities

by application of a DBEM [34]. The effects of the acoustic

loading were included as forces on the structure in an

inhomogeneous Helmholtz formulation (i.e., a forced

acoustical problem). An overdetermination procedure similar to that used by Schenck [ 19] was used whereby the 

resulting pressures were a least squares solution of the

problem. Constant acoustic variable variation was implied

on an element surface. Impedance boundary conditions were

also included. A series of studies involving a cylindrical

cavity excited by a harmonic driver at one end produced

very good agreement with experimentally obtained results.

Other studies investigating the quantity and placement of

s
31

1

the interior overdetermining points were also presented

but will not be discussed here.

4 

2.3 Conclusions

This chapter has presented a number of references

= related to prediction of sound fields in cavities and 

boundary element applications in acoustics. It reviews the

status of techniq'.,es being applied to the problems within

the co cern of the current research. In addition to the

technical literature, there are good academic references

available relevant to boundary element tect.niques for both

general and specific topics [12,35,36].

A few developments from the literature review should

be emphasized due to their significance to the current

9

i

research. IBEM's require less information than DBEM's.

Equations 2.1 and 2.2 illustrate that less information is

required in the integral equations. More importantly, only

one quantity, the surface source distributiot ► , is initially computed at the boundary whereas, for DBEM's, both 

pressure and its gradient, depending on the boundary condition, are resultants. Hence, the IBEM may be prefered

for noise source identification since interpretation of

the results would be more straightforward. Impedance

boundary conditions have been applied in past research and

f

' will be applied in the current research. The information

relating to the failures of the boundary element
32

formulations may provide guidance in case difficulties of

this nature are encountered. The literature review has

shown that isoparametric representation, which is to be

utilized in the current research, has been previously

utilized successfully. One very Important concept

discovered during this review is that of domain transformation to integrate around the singularity inherent in the

Helmholtz integral formulations. A similar approach will

be used and is detailed in chapter 4. Finally, the literature review shows a definite lack of boundary element

applications to acoustical cavity problems. Of the 21

boundary element references cited, only 3 were concerned

with cavities.
THEORY

The problems to be considered are those dealing with

acoustical behavior within cavities. A description of the

problem is shown in Figure 3.1. A point within the domain

D is located by the vector x. Likewise, a point on the

boundary B is defined by the vt:tor C. The vector normal

to the boundary is signified by n and is considered to be

directed outwsrd from the cavity. An applied distributed

source of strength IP within the domain is located by the

vector x .

s

s
IL;

34

I LN.

Figure 3.1 - Problem Description
35

For a problem of this type, the governing equation is

the familiar non-homogeneous Helmholtz equation.

o20 + k20 - V(x5

) (3.1)

1 where x is the vector locating the applied source. Equation 3.1 is the linearized, lossless Helmholtz equation

i

formulated in terms of velocity potential, It. The variables of acoustic pressure and particle velocity can be

related to the velocity potential through the relationships (37)

p (3.2)

v - off. (3.3)

A Huygen's principle formulation which may be used to

solve the homogeneous form of equation 3.1 has been

developed by Chen and Schweikert (13]. In this f.)rmulation, a distribution of simple sources is considered to

exist at the boundary of the cavity in question. The boundary source magnitudes are determined in such a way that

the boundary conditions for the problem are satisified.

Chen and Schweikert's method is essentially a simple

indirect boundary element method. '"o develop a more complete indirect boundary element method, a variable distribution of sources may be assumed as will be shown in

Chapter 4. The velocity potential at any point due to the

assumed source distribution is

,_r,;
Wr

1

36

e -j kr

'*(x) - fa (C) r dB (3.4)

B

where r is the distance between the points x and C and

O(C) is the source strength. Although equation 3.4 is

written for a point in the domain, it applies equally as

well for a point on the boundary. Similarly, the gradient

of the velocity potential is

-jkr

vm(x) - 47TcQ(x) + f o(C)v 

e 

r dB (3.5)

B

where c results from a singularity integration which will

be defined further in this chapter and B is the boundary

excluding a smail region at the singularity.

At this stage, a definition of two quantities will be

made to ease notational difficulties. Combination of equations 3.2 and 3.4 allows an integral representation of the

pressure at a point to be expressed as

e-jkr

P(X) - - jc' pfc (C) = dB. (3.6)

B

A quantity known as the fundamental pressure solution, p

is defined by

* _ _ e-jkr

P (C,x) - -jwp r (3.7)

Hence, equation 3.6 may be more concisely written as

P(X) - fo(C)p *(C,x) d8. (3.8)

B

t
37

Likewise, the fundamental velocity solution, v * , may

t

be defined by considering . the partial derivative within

the integral portion of equation 3.5, that is

E . r

v r (3.9)

e-jkr

Since, for a well-posed boundary value problem, the component of velocity of interest is the normal velocity, the

following relationship from Kaplan [ 38] is used.

ao v() n (3.10)

an

Using equations 3.9 and 3.10, the fundamental velocity

solution may be defined as

v

* (C,x) - n• or —2 - z ^e-3kr

(3.11)

r

Equation 3.5 may now be written more concisely as

v(x) - 4nca(x) + f a(C)v*(C,x) dB (3.12)

B

where v is the normal velocity component of the velocity

• potential gradient in equation 3.3. 

The quantity c in equation 3.12 results from integration around a singularity point. From examination of the

fundamental solutions, it is evident that as the distance

between the two points of concern approaches zero the fundamental solutions become infinite. However, the fundamental pressure solution is a i singularity while the fundamental velocity solution is a 

r 2

n si gularity. The r

k

f

f
r

38

singularity is weak and may be integrated but the 2

r

singularity is strong and must to treated as a Cauchy

principal value [40]. The basic premise is to exclude a

small region surrounding the singularity from the integration and account for the exclusion with the added 'free

term' (i.e., the quantity c) [41,42]. Evaluation of the

'free term' for an interior domain can be donee with the

following rules.

y c 0 for an exterior point

• c - -1 for an interior point

c - -Z for a point on a smooth boundary

c - -4n for a point: on a corner

A smooth boundary is any boundary with a unique tangent

plane. For a corner point, 0 is the value of the eDlid

angle included in the domain at the corner expressed in

steradians.

The previous discussion is for the homogeneous problem but the inhomo"±eneou- problem can be rolved by augmenting equations 3.8 and 3.12 with a domain integral

[39]. The domain integral represents the contribution of

an applied distributed source of strength *. After inclur

ino the domain integral, the equations become

s t

t

ZIV_ :
39

`

p ( x ) fO(C)p (C,x) dB + f*(xs) p (xs ,x) dD (3.13)

B D

and

v(x) - 4nca(x) + f v(c)v*(C,x) dB +

Bf*(xs )v (x s ,x) dD (3.14)

"t D

where x s indicates the location of the applied source.

If the applied sourcea are considered to be monopole

sources, evaluation of the domain integral is significantly simplified. For a point source, the domain integral

is easy to evaluate as the point Pource may be represented

by a Dirac delta function multiplied by the source

strength. Thus, the integral exists only at the point of

application of the source. The domain integral becomes

_ * nsource fv( x s )P (x s ,x) dD - E fti O(x si - x)p (x s ,x) dD =

D i-1 D

nsource

E 4rlip (xsi'x) (3.15)

i-1

for a pressure boundary condition and

:source

fv( xs ) v (x s ,x) dD - E ft i a(x si -

x)v (x s ,x) dD -

D i-1 D

nsource

E 4nt i v (x si ,x) (3.16)

isl

for a velocity boundary condition where nsource is the

F

total number of applied sources. Distributed internal

sources are not difficult to evaluate but are rare.

i
f

40

Vibrating sources which do not behave as point sources can

be treated as additional boundaries.

Impedance boundary conditions may also be applied

with this IBEM. Specific acoustic impedance is defined as

the ratio between pressure and particle velocity (43] or

z -P. (3.17) V'

Manipulation of equation 3.17 yields an equation which is

applicable to the integral equations.

p - zv - 0 (3.18)

In terms of the integral equations, impedance boundary

conditions may be formulated as

-49zca(x) + fo(C)Ip * ( C,x) - zv * (C,x), dB +

d

necurce

E 4n'Pi(xs)(p*(x ai l

'

) - zv (x ail X)] - 0 . (3.19)

i-1

Equation 3.19 can be solveC for the unknown fictitious

source distribution Q.

All the integral equations have now been developed 

for this indirect boundary element formulation. Depending

on Lhe type of boundary condition applied, equations 3.13,

3.14 and/or 3.19 are used to solve for the assumed source

density distribution, o, on the boundary of the body in

question. Once the source dersIty distribution has been

evaluated, the same three equations may be applied to

solve for both the unknown boundary conditions and the
41

acojatic variables at desired points within the domain.

The numerical implementation of the precceding theory is

presented in the next chapter. For the purpose of thR

current research, only the acoustic pressures at the boundary and field points will be evaluated as those quantities adequately describe the sound field within the cavity.
42

CHAPTER 4

NUMERICAL IMPLEMENTATION

ng chapter describes the theory of the indirect

ement formulation of the cavity acoustics probplement the theory into a usable computer pronumber of numerical analysis techniques were

is chapter presents the details of the numerical analysis techniques which were used for the implementation.

4.1 Element Definition

The boundary element incorporated into this research

is an isoparametric element with quadratic shape functions. The term isoparametric means the same interpolation

function is employed for both the geometric approximation

and the acoustic variable approximation. Pox an isoparametric element, the number of nodes is equal in both
43

`i

,f

--.

the local and the global coordinate systems. Consequently,

a node in the local coordinate system correlates directly

to a node in the global system. The local element is shown

in figure 4.1. kn example of a global element is shown in

figure 4.2.

4 7 3

f

8 6

9

1 5 2

Figure 4.7 - Local 1soparametris Element
44

Y

X

Figure 4.2 - Global Isoparametric Element

The important feature of the element is the shape functions. The shape functions are used to interpolate the

variation of the variables within an element. For the isoparametric element, one set of functions describe the

variation of both the element geometry and the solution

variaoles. The element in figure 4.1 has a 3x3 nodal distribution in a two-dimensional space (i.e., the E and the

n directions). Consequently, a quadratic variation can

occur. The interpolation procedure for a general variable,
4 5 

41 , is

9

- E N1 0 i

X4.1)

1-1

where i is an index which varies with the node number, Ni

is the shape function for the ith node and 0 i is the value

of the variable at the ith node. The Lagrangian shape

functions for the quadratic, isoparametric element on the

parent shape are (44]

N l - 4( 4 2 - E)( 71 2 - 71) (4.2a)

N2 = 4( E 2 + E)(,12

- 71) (4.2b)

N3 - 4(E2 +

E)(712

+ 71) (4.2c)

N4 -

4(E2 - 

4)(11 2 + 71) (4.2d)

N 5 - 2( 1

- E2) (,12 - 11) (4-2e)

N6 - 2(E2 + E)(1 - 11 2 ) (4.2f)

N7

 - 2( 1

- 4 2 ) ( 11 2 + n) (4.2g)

14 8

= 2( E 2 - E)(1 - n 2 ) (4.2h)

N9 - (1 - E 2 )(1 - 11 2 ). (4.2i)

The Lagrangian shape function, N i , is a continuous function having a value of unity at node i and a value of zero

at Node j*i. In the numerical iwplementation, unknown

variable distributions will be replaced by an approximation with the form of equation 4.1. By this approximation,

a finite number of parameters may then be solved.
46

4.2 Equation Discretization

Having defined the element, the next stage of numerical implementation is to discretize the integral equations. The discretization concept in boundary element

methods is similar to that in finite element methods. The

concept is that by subdividing the domain of the problem

into smaller subdomains, or elements, the integrals for

the complete domain may be approximated on a piecewise

basis. Furthermore, by enforcing certain boundary conditions between the subdomains, an approximate solution over

the entire problem domain may be determined.

The matrix equations are assembled by satisfying

boundary conditions at an appropriate number of discrete

points. The specific points used in this development are

the nodes of the model. If the given boundary condition

is a pressure boundary condition, the pressure at that

point in terms of the boundary source distribution, a, is,

by equation 3. 13,

* K

P 1 - Jap (r) dB + E 4v*kpkl. (4.3)

B k=1

where p l is the pressure at the 1 t boundary condition

point and pkl is the fundamental pressure solution from

the kth applied source to the 1 t boundary condition

point. The unknown boundary source distribution, a, on

each element is «pproximated by

P

t
47

9

a = lElaiNi (4.4)

_

utilizing equation 4.1. Furthermore, the boundary integral

term will be evaluated on a piecewise basis (i.e., on an

elemental basis). Hence, equation 4.3 can be rewritten as

N M * K

p l : E f E a i N i p (r l ) dBn + E 4n*kpkl (4.5)

n=1 B i=1 k=1

where N is the number of boundary elements, M is the

number of nodes per element, K is the number of applied

sources, r l is the variable distance from the 1 t boundary

condition point to the n th element and B represents the

surface of the nth element. Since the boundary source

values, a i , are constants, the integral in equation 4.5

may be rewritten as

N M

E f ( E a i N i )p (r l ) dBn =

n=1 Bn i-1

M N #

iEl n lai E B 

i N p (r l ) dBn E aiail (4.6)

n

The numerical evaluation of the a il term will be discussed

in sections 4.3 and 4.4. Finally, equation 4.3 can be more

simplistically rewritten in approximate form as

K

p l : E aiail + E 4n*kpk1 (4.7)

k=1

A similar approach can be used for the velocity and

impedance boundary conditions. With these two types of
' 48

t .

boundary conditions, an additional Lerm due to the singularity integration is present. To show the consequence of

the 'free terms', the equation discretiza%ion for the

velocity boundary condition will be presented. In equation

3.1,4, the velocity in terms of the boundary source distria buti(n, is given as

vl - 4nco + f av * (r 1 ) dB + 

K 

4n*kvkl . (4.8)

B k=1

Evaluating the integral on a piecewise basis yields

N M

vl _ 4nc i1° i + E f E a iN iv ( r l ) dBn

n=1 B i=?

n

K

+ E 4"1 vkl

(4.9)

k=1

where the variables maintain the same meanings as in equation 4.5. The 'free term', c il , is zero when i 0 l. When

i=1, c il is -2 for a point on a smooth boundary and -42

for a point on a corner where n is the solid interior

angle at the corner. The constant source distributions can

be factored out of the integral.

N M

477c ila i + E f ( E a i N i )v (r l ) dBn -

n=1 B i-1

n

M N

E 

o i (4nc. + E f N i v (r l ) dBn =

i=l n=1 B 

E °i (41c il + a il ) (4.10)

The final approximation equation for a velocity boundary

condition is
49

M K x

v

i

E o i ( 4 9c i1 + a il )

+ E 49*kvkl (4.11)

k=1

Likewise, for impedance boundary conditions, the integral

equation is

.t

f oi(px(r 1) - z lv x ( r l )) dB - 41zlcil°i

t Bn

K x x

+ E 4n'rk (Pkl - z lvkl ) ' 0 (4.12)

k-1

on a piecewise basis or

K x x

E Q i (-412 1c il + a il ) + kEl4 ntk (Pkl - z lvkl ) = 0(4.13)

for the approximation equation.

For all boundary conditions, a general matrix equation can be wriiien as

1 A+ A' 1 19 1 + I D I j a l (4.14)

where either equation 4.7, 4.10 or 4.13 are used as each

row of the matrix. The evaluation of the matrix terms for

each boundary condition type are:

1. For Pressure Boundary Conditions

a il - f Ni p x ( r l ) dBn

(4.15a)

B 

al 0 (4.15b)

dk - E 4n`Pkpkl

(4.15c)

k=1

a 1 P1

(4.15d)

^
R

Y,

n.:

W

50

2. For Velocity Boundary Conditions

ail = f N i v (r l ) dBn

(4.16a)

B

n

aril - 4nc il

(4.16b)

K

dk - E 49tkvkl

(4.16c)

k-1

al - V (4.16d)

3. For Impedance Boundary Conditions

a il = f ( P * (r l ) - z lv * (r l ) ) dBn

(4.17a)

B

n

a@ il - -4nz i c; l

(4.17b)

K

dk M

k=E 1

4ntk ( Pkl - z lvkl ) (4.17c)

a l

0 (4.17d)

In the D matrix terms, the fundamental solutions are from

the kth applied source to the 1 t node. The theory can be

implemented into a computer format using equations 4.14

through 4.17.

4.3

To numerically evaluate the integrals in equations

4.15a, 4.16a and 4.17a, a numerical integration procedure

must be utilized. A numerical integration procedure using

a Gauss-•Legendre quadrature was selected. The numerical

evaluation of the integral equations is performed on a

e

F

Q
51

piecewise basis by mapping the integral in the global element coordinate system into the local element coordinate

system. The domain of the local element is defined between

the f and V coordinates of -1 and 1. The limits are consistent with those of a Gauss-Legendre quadrature as the

quadrature formula is (45]

1 NG f g(f) di = E g(f i )w i (4.18)

-i i=1

for a one-dimensional integration of a general function,

g(f). The values E i and w i are the integration point and

weighting factor, respectively, and NG is the total number

of integration points used in the integration. A boundary

integration is essentially an area integration (when considering three-dimensional problems) which may be

represented by

1 1 f g(E,71) dB -

n f f g((,71) d071 (4.19)

A -1-1

n

For a two-dimensional integration, the quadrature formula

is (46]

NGE 

f 

NGn

B

g(f 71) dB - 

mEl 

I nEl 9(F n ► '1 m ) wn I

J 

wm

(4.20)

n l 

Equation 4.20 was used to numerically evaluate the boundary integrals everywhere except where a singularity

exists in the fundamental solutions as will be dicussed in

section 4.4.

4
52

The location of the integration points and the values

of the weighting factors must be determined to finalize

the local element integration procedure. A 4x4 Gaussian

quadrature having coordinates

a

E 1 - 71 1 - -0.861136311594053 (4.21a)

f2 - 71 2

-

-0.339981043584856 (4.21b)

E 3 - 71 3 - 0.339981043584856 (4.21c)

E 4 - 7) 4

- 0.861136311594053 (4.21d)

was selected. The mesh is shown in figure 4.3. The weighting factors corresponding to the integration points are

w - 0.347854845137454 (4.22a)

w2 - 0.652145154862546 (4.22b)

w3 - 0.652145154862546 (4.22c)

w4 - 0.347854845137454. (4.22d)

The abovementioned integration quadrature is exact

for a seventh-order or lower polynomial, g(E). The functions being integrated in the boundary integrals are the

* s

fundamental solutions, p and v times the shape functions

,N i , which are quadratic polynomials. The fundamental

pressure solution is

* e-jkr

p (C,x) - - jWP r (4.23)

The series expansion of the fundamental solution is
« ^i

j

53

Wr

t I

Figure 4.3 - Integration Point Mesh

* 2 4 3 6 5 8 7 1 k r k r k r k r p (C,x) - -jcip ( ( r - 2! + 4^ - 6, + y—) + ...

_ k32 k r54 kr7 6

-j ( k 31

r

+ 51 - 7! ) + ... ) (4.24)

Hence, from equation 4.24, some of the limitations of this

integration can be illustrated. The function will be accurately integrated only if the first few terms of the

series approximation dominate. The r term poses a problem
54

near a singularity point and will be addressed in the next

t

section.

The preceding discussion applies equally as well for

the fundamental velocity solution but will not be

presented here.

The final point to be made for the integration procedures involves transformation from the global domain.

Since the integration is for the global element, a

transformation between the domains is necessary as follows

(47)-

fG(x,y•z) dA ' fg(4,71),'J(f,71); dfdn (4.25)

In equation 4.25, the function G(x,y,z) is a general function written in terms of the global element coordinates,

x, y and z. The function g(f,n) represents the transformation of G(x,y,z) into the local coordinates, f and n. The

quantity ;J(f,n)1 is the determinant of the Jacobian

matrix J(f,n). For the application at hand, a threedimensional global element is mapped into a twodimensional local element. The determinant of the Jacobian

for this transformation is [48]

1

fJ(f,n)f - ( a l + a2+ a3 )2

where

a y az

` a l ( - z a ) of an an of

(4.26)

(4.27a)

^
^^.^

_ 

55

a 2 - ( aE a - d a^ ) (4.27b)

a 3 - ( aE a - at a^ ) (4.27c)

The partial derivatives in equation 4.27 are determined

with interpolations such as

9 am.

alEl xi at

as an example.

(4.281

The procedure described above can be applied to any

element which does not contain a singularity point. when

the value of r in the fundamental solution approaches

zero, special considerations are necessary for the

integration. These considerations are explained in than

following section.

4.4 Singglarity Numerical Integration Procedure

A singularity occurs in the fundamental solutions

whenever the particular boundary condition under consideration is located on a node of an element over which

the boundary integral is being evaluated. The variable r

in the fundamental solutions is the dia'ance between the

boundary condition point and the integration point. The

primary Gaussian quadrature is unable to adequately evaluate the integrals whenever this singularity occurs. Hence,

a new quadrature had to be developed to handle the singularity situation.

P
56

In some cases, the singularity can be removed from

the integration domain. As detailed in the previous section, the integration domain is a quadrilateral region in

a rectangular coordinate system. However, if the region

were considered in polar coordinates, the differential

area is

dA - rdrd8 . (4.29)

The radius appearing in equation 4.29 may be used to cancel the r factor in the fundamental solution. Consider,

for example, the fundamental pressure solution.

- j kr

-jwP f 

e 

r UJ5 = -jwP ff e -3kr drde (4.30)

n n

A similar result occurs for the fundamental velocity solution.

f A.

rl - jk )e 

]kr drde (4.31)

B

r.

In equation 4.31, it is apparent that a r singularity

still exists even after transforming to a polar coordinate

system. However, the dot product of the vector normal to

the surface and the gradient of the radius is zero for a

flat place element and tends to zero as the radius of curvatu-e approaches infinity for curvilinear elements. For

example, in figure 4.4, the fundamental solution is

evaluated near the u pper right corner. The n *vr will tend

to zero. Thus, the integral in equation 4.31 is integrable
57

as a Cauchy principle value. (Hence, the 'free term' mentioned in chapter 3.)

Figure 4.4 - Vectors For A Flat Element

There are three singularity integration cases which

arise. Those cases are:

1. Boundary condition at a corner node

2. Boundary condition at a midside node
58

3. Boundary condition at the center node

One polar domain transformation can be developed to accommodate all three cases. The transformation for a quadrant

will apply for a corner node singularity. In addition, the

remaining two cases can be considered as multiple quadrants and utilize the transformation for the corner node

singularity. The transformation can be applied to both

quadrants neighboring a midside node singularity and all

four quadrants neighboring the center node singularity.

Figure 4.5 shows the three cases. The shaded areas

represent the quadrants transformed to a polar domain.

t

5 F

Y
59

a

Corner node singularity

Midside node singularity

Center node singularity

Figure 4.5 - Singularity Cases
In general, a Gauss-Legendre quadrature similar to

that explained in the previous section can be used for the

polar domain. The difference lies in the integration variables, in this case r and 6.

The integration points in the polar coordinate system

can be transformed from the rectangular coordinate system.

The values of the f and n coordinates range from -1 to 1.

In the circular sector, the radial coordinate ranges from

0 to 1 and the angular coordinate from 0 to 2. Mapping

functions can be established which allow the polar coordinates to be obtained from the rectangular coordinates.

These mapping functions are

1

r - 2 (f + 1) (4.32a)

for the radial coordinate and

e 4( 71 + 1) (4.32b)

for the angular coordinate. With equations 4.32, the polar

integration points become a function of the rectangular

integration points.

Determination of the polar weighting factors requires

more work than for the determination of the polar integration points. Weighting factors are evaluated by the equation

b n (f - fj)

w i = J II df (4.33)
{

Ir,''

for a Gauss-Legendre quadrature. To illustrate the determination of the polar weighting factors, conaidex a two

point quadrature with points E 1 and E 2 in the rectangular

system and points r 1 and r 2 in the polar system. The equation for the first weighting factor in the rectangular

system is

1 E - E2

wE - f - E 

U46. (4.34)

1 -1 1 2

For the first weighting factor in the polar system, the

equation is

1 r - r2 f _ dr. (4.35)

Wr 1 0 r 1 r2

Using equation 4.32 with equation 4.35, the weighting factor may be expressed as

1 2(E + 1) - 2( E 2+ 1)

w

1`c

r f 1 1

(2a.) (4.36)

1 -1 2(E1 + 1) - 2 ( E 2 + 1)

or

1 1 (E + 1) - (E 2 + 1)

wr i - 2-1 ( .C l

+ 1) - (E 2 + T)- d (4.37)

or

1 1 E E2

1

wr - 2f 1 C l - 2 dE

(4.38)

Equation 4.38 may now be written as

wr - 2-wE (4.39)

1 1

by substituting in equation 4.34. The resulting
relationship is independent of the number of integration

points but dependent. solely on the multiplicative factor

in the - mapping functions (i.e., equations 4.32). Consequently, the weighting factors for the polar domain are

wr - 2 (4.40a)

e

for the radial weighting factor and

w9 4 71 (4.40b)

for the angular weighting factor. The weighting factors

and the integration points in the rectangular domain may

be found in any reference containing Gaussian quadratures

(49].

In the circular domain, a 4x16 mesh of integration

points was selected; 4 radial coordinates and 16 angular

coordinates. The integration points and weighting factors

for the four point scheme are shown in equations 4.19 and

4.20. For the lc point scheme, the integration points and

weighting functions are

7? 1 - -0.98940093 0,991649 (4.41a)

'1 2 - -0.944575023073232 (4.4lb)

7? 3 - -0.865631202387831 (4.41c)

71 4 - -0.755404408355003 (4.41d)

7? 5 - -0.617876244402643 (4.41e)

71 6 - -0.458016777657227 (4.41f)

71 7 - -0.281603550779258 (4.41g)

E 
71 8

- -0.095012509837637 (4.41h)

71 9 - 0.095012509837637 (4.411)

71 10 ' 0.281603550779258 (4.41j)

7111 0.458016777657227 (4.41k)

71 12 - 0.617876244402643 (4.411)

7113 0.755404408355003 (4.41m)

7114 0.8656312n2387831 (4.41n)

7115 0. 8 44575023073232 (4.41o)

71 16 ' 0.989400934991649 (4.41p)

w l - 0.027152459411754 (4.42a)

w2 0.062253523938647 (4.42b)

w 3 = 0.095158511682492 (4.42c)

w4 = 0.124628971255533 (4.42d)

w5

-

0.149595988816576 (4.42e)

w6 - 0.169156519395002 (4.42f)

w 7 ' 0.182603415044923 (4.42g)

w8 = 0.189450610455068 (4.42h)

w 9 = 0.189450610455068 (4.42i)

wl0 = 0.182603415044923 (4.42j)

wll ' 0.169156519395002 (4.42k)

w12 ' 0.149595988816576 (4.421)

w 13 ' 0.124628971255533 (4.42m)

w 14 = 0.095158511682492 (4.42n)

w 15 ' 0.062253523938647 (4.420)

w 16 ' 0.027152459411754. (4.42p)

• 
and

The mesh for the circular quadrant is illustrated in figure 4.6.

Examination of figure 4.6 reveals a gap between the

circular domain and the original rectangular domain. As

shown, the boundary integration for the circular domain

would exclude that portion of the boundary outside the

circular domain yet still within the original rectangular

domain. This omission was corrected with a second domain

transformation from the circular domain to a rectangular

domain. The values of the f' and 71 1 coordinates range

from -1 to 0 since a quadrant of the element is being considered. Note that primes are used to differentiate the

rectangular coordinates in the second transformed domain

from the original rectangular coordinates. The mapping

functions for the integration points are

E' r - 1 for 9 < n4

(4.43a)

f' - rcote - 1 for 9 > 4 (4.43b)

rtane - 1 for 9 < 4 (4.43c)

r - 1 for 9 > n4 (4.43d)

For the weighting factors, the mapping functions are

f' r

(4.44)

and

we

-

►1' Cosa for 6 4 n4 (4.45a)
V.

Figure 4'6 - Circular Domain Mesh

- '-
66

we

w7l ain9 for 8 > 4 (4.45b)

These mapping functions were determined in the same manner

as the mapping functions for the first domain transformation (i.e., rectangular to polar). Figure 4.7 shows the

mesh for the transformed circular domain.

WC 
I

0

Figure 4.7 - Transformed Circular Domain Mesh

i,j
fit- r

The final term required for the numerical integration

is the determinant of the Jacobian. For a normal integration procedure, the determinant of the Jacobian is computed using equations 4.26 and 4.27. The determinant of

the Jacobian can also be viewed as a ratio of the global

area to that of the local area.

dxdy

(4.46) didn

In polar coordinates, equation 4.46 is rewritten as

r dr de

r _g (4.47)

where the subscripts indicate the global and local domain

variables. For the singularity integration., the integral

is of the form

ff(rg ,eg )dr 9 do9 -

where ;J; r is the determinai

the singularity integration

the required determinant of

fF(r l . e l )iJ= r d r l d e l (4.48)

it of the Jacobian required for

procedure. From equation 4.48,

the Jacobian will be

dr de

r = dr99 (4.49) 1de 1

and can be calculated from the normal determinant of the

Jacobian by

r

;J; r - rl ;J; (4.50)

9

If the mapping between the global and local domains is

simply a scaling of the same angular dimension sector, the
69

6, mapping equations can be written as

r g - art (4.51a) 

for the radial component and

eg - 9 1 (4.51b)

for the angular coordinate since it is assumed that the

element is only varying in size and not in shape. Thus,

the differentials are expressed as

dr g - adr 1

(4.52a)

and

deg - de l

(4.52b)

For this case, the normal determinant of the Jacobian is

IJ, - 

r9dr9de9 - (ar1)(adr1)del - al

r ldr lde l rldrld8l

(4.53)

and the required determinant of the Jacobian is

dr de (adrI)del

z dr-de drld9l - a (4.54)

Therefore, the required oetermina-..t of the Jacobian may be

calculated from the normal determinant by

j ;J;r (4.55)

The relationship in equation 4.55 is an approximation if

tE. the shape of the global element is distorted from that of

the local element. Nevertheless, this relationship is

utilized for the singularity integration. Hence, the quadratuxe formula becomes
70

ff(x,y)dxdy - fF ( r l re 1

) -,J ;J;dr l d9 l . (4.56)

The numerical integration procedures detailed above

cd^ be utilized to evaluate the discretized integral equations to establish the system matrices for a boundar y element model. Although this numerical integration procedure

is successful, it is not the only possible procedure to

apply. An alternative integration procAdure has been suggested by Zimmerle and Bernhard (50) but was not att,::i-+ted

in this research.

4.5 Solution Procedure

Equation 4.14 is used to solve for the assumed source

densities once the numerical procedures described above

have been used to evaluate and assemble the system

matrices. To solve for the source densities, a full

matrix, complex equation solver is required. The matrix

equation is neither banded nor symmetric so a robust equation solver is required. The subroutine used, the IMSL

subroutine LEQTIC, is a linear equation solver for equations having complex quantities. It uses a full matrix

storage mode. This storage mode is required since thmatrices for boundary element methods, unlike finite eiameat methods, are fully populated matrices. The IMSL subroutines are readily accessible as a commercial code which

was the primary reason for the application in this
71

research. It is possible that a more efficient equation

solving routine is available or may be generated for use

in solving the boundary integral equations.

The goal for this research was to predict the sound

fields in acoustical cavities. To accomplish this, the

source densities, after quantization, were used to solve

for Lhe acoustic pressures at field points interior to the

cavity domain. Equations 4.14 and 4.15 were again utilized. The unknown variables, the pressures at the field

points, can now be solved utilizing the source density

distribution found earlier. The only major difference in

the formulation is V.- .he distance r in the fundamental

pressure solution is the distance from the source density

to the field point at which the pressure is desired. Using

this procedure, the pressures at the interior field points

may be evaluated.

The procedures required for numerical implementation

of the boundary integral equations have been discussed. To

illustrate the techniques, a number of examples were considered. The numerical results of these examples are

presented in the following chapter.
4

l

CHAPTER 5

CASE STUDIES

A number of case studies were undertaken to illustrate the theory and numerical implementation discussed in

chapters 3 and 4. Two models were considered in these studies: a spherical cavity model and a rectangular cavity

model. This chapter presents the numerical and, where

available, the analytical results of the case studies.

5.1 Spherical Cavity Model

A spherical cavity was chosen as the first case study

for a number of reasons. A sphere with a velocity distribution which is uniform with respect to angle can be

treated as a one-dimensional problem as will be demonstrated in the derivation of the analytical solution.

Consequently, the analytical solution for the spherical

cavity can be easily determined. Another advantage for the

sphere is the absence of corners and edges. Corners and

r

^ ter
73

create difficulties with boundary element methods as

)e shown for the rectangular cavity.

'he analytical solution for the spherical cavity is

!d from the wave equation for spherically symmetric

ire fields [511.

art EaP 1a P+ _

2

ar r - c2 

at2

(5 1)

ing the time dependence from equation 5.1 yields the

Lon

2

2--2 + ? -^R + k 2 - 0 (5.2)

a r t r ar

represents the time-independent wave equation in a

Lcally symmetric geometry without applied sources.

equation is a special case of the spherical Bessel's

ion

2

r2 ^ + 2r^ + k2 (r 2 - n2 )p - 0 (5.3)

ed by r 2 and with n- 0. Thus, the solution to equation

ill involve spherical Bessel functions and will be of

)rm

p(r) - Aj 0 ( kr) + By0 (kr) (5.4)

the constants A and B are determined by the boundary

:ions. Alternatively, a second general solution to

, )herical Bessel's equation of equal validity is

,- a -jkr _ ejkr

p(r) - a r + b r (5.5)

E:
} s

74

where a and b are determined by the boundary conditions.

For the first verification of the numerical code, the

specific solution of

P(r) - 100j 0 (kr) (5.6)

was selected.The spherical model has a radius of 2 meters.

Also, an arbitrary wave number of 8 was selected. This

corresponds to a frequency of 437 Hz in air with a speed

of sound of 343 meters/second. Figure 5.1 shows the

analytical solution of the pressure distribution within

the the cavity. The sphere center corresponds to a radial

position of zero while the sphere boundary corresponds to

a radial position of 12'

do
75

L, S

160 .0

i tA

m 0.000

L

CL 

--180.0

0.00 .050 .100 .150 .200 .250 .300 .350 .400 .450 .500

120.00

100.00

0) 80.000

7

•+ 60.000

C

a,

40.000

20.000

0.0000

0.00 .050 .100 .150 .200 .250 .300 .350 .400 .45J .500

Radial Position

Figure 5.1 - Analytical Spherical Pressure Distribution

For k-8
5

76

Appropriate boundary conditions of various types were

derived to verify the solution in the cavity interior. The

solution of equation 5.6 at the boundary is p( 2) - -18.92

with the selected wave number. For the first example,

this pressure was applied at the outer boundary of the

sphere. The numerical results are shown with the analytical results in figure 5.2. The numerical solution is

represented by the solid line marked with Y's at the

discrete solution points. The dashed line represents the

analytical solution.

A sphere with uniform velocity boundary conditions

was the second numerical example. The velocity boundary

conditions were determined by

V(5.7)

which relates the velocity of a spherical wave to the

pressure (51]. For the numerical example, the velocity

boundary condition was v - -0.0456 + j0.0114. Figure 5.3

shows the numerical and analytical results.

D
180.0

a

IA

0.000

t

Q_ -180.0

0.00 .050 .100 .150 .200 .250 .300 .350 .400 .450 .500

120.00

100.00

U 80.000

7

.+ 60.000

C

a,

E 40.000

20.000

0.0000x----r— x--

0.00 .050 .100 .150 .200 .250 .300 .350 .400 .450 .500

Radial Position

77

Figure 5.2 - Spherical Pressure Distribution For Pressure

Boundary Conditions At k-8
78

E

160.0

QJ I

^00.000 - -

C

a. 

-180.0

0.00 .050 .100 .150 .200 .250 .300 .350 .400 .450 .500

120.00

100.00

Qr 80.000

7 ^

«a

.. 60.006

C

O+

E 40.000

20.000

+7.0000

0.00 .050 .100 .150 .200 .250 .300 .350 .400 .450 .500

Radial Position

Figure 5.3 - Spherical Pressure Distribution For Velocity

Boundary Conditions At k-8
79

The analytical and the num-tr ical results in both f igures 5.2 and 5.3 are in excellent agreement. In fact, the

correlation is so good for figure 5.2 that the two results

are nearly coincident. One note must be made regarding the

phase plot. When examining the phase, it must be understood that a phase of 180° is equivalent to a phase of

-180°. Hence, the apparent jumps in phase in the region

from r - 0.4 to r - 0.5 do not actually represent phase

differences as it may appear at first glance. For the

velocity boundary condition model represented in figure

5.3, a slight amount of difference can be detected between

the analytical and numerical solutions. This slight error

is due to the interpolation of the radial distances

between the boundary condition points and the integration

points. As previously noted, the shape functions are quadratic functions. However, the radial distance is determined by

I - ( x2 - x 1 ) 2 + (y2 - y l ) 2 + (z 2 - z 1 ) 2

(5.8)

which is not a purely quadratic function due to the radical involved in its determination. Consequently, some

interpolation error can be expected and, in fact, is evident. The interpolation error is more prevalent when the

velocity boundary conditions are used because of the 1

r

2

term in the fundamental velocity solution as opposed to
80

j

the term in the fundamental pressure solution. Regardless, for cases where derived boundary conditions are

specified, the prediction of the solution in the domain is

very good.

The next case study was that of a pulsating spherical

cavity. The analytical solution for this problem can be

represented by equation 5.5 where the boundary conditions

for this problem are

ar . 0 (5.9a)

ap t the center of the sphere and

a r -J(JPV (r o ) (5.9b)

at the surface of the sphere. Boundary condition 5.9a

respresents the condition that the pressure is finite at

the sphere center. Equation 5.9b is derived from Eular's

equation. Using equation 5.5 with the boundary conditions

of equations 5.9, the analytical solution for the pressure

is

pckr2V(ro

sin kr

)

P(r) ^(sin(kr o ) - kr 0 cos(kr o )) r

(5.10)

For the model in question, the radius at the sphere surface is r - 1

Resonant frequencies can be determined for the spherical cavity using equation 5.10. Resonance occurs whenever

the impedance vP becomes infinite. From equation 5.10,

...^-. a -

_ _ .

_ _ -- •^-^ -fir'.- .
^i

81

resonance for the spherical cavity occurs whenever

sin(kr 0 ) - kr 0cos(kr C; = 0 (5.11)

The transcendental equation for the natural frequencies is

tan(kr 0 ) - kr 0 (5.12)

for the spherical cavity. Thus, the first three natural

frequencies occur at kr o - 0, kr o - 4.493 and kr o - 7.725.

The analytical solution of equation 5.10 at r - r 0 is

shown in figure 5.4 as a function of the frequency. The

natural frequencies are evident in figure 5.4 as peaks in

the pressure magnitude. In addition, figure 5.4 shows zero

pressures and 180 0

phase shifts occuring at

kr o - n and xr o - 2n. The zero pressures are expected

since p - 0 whenever sin(kr 0 ) - 0 in equation 5.10. The

velocity at the boundary of the sphere was chosen to be 1

to generate the analytical solution in figure 5.4.

To numerically evaluate the response of the spherical

cavity, a sweep of frequencies from 0 to 873 Hz was used 

tL generate results comparable to those in figure 5.4. The

frequency range corresponds to kr o values from 0 to 8. A

uniform velocity boundary condition of 1 was applied to

the spherical cavity model. Both the analytical and numerical results are shown in figure 5.5. The solid line

marked with Y's represents the numerical results and the

dashed line is the analytical solution.

t
_. -_ Orr

82

180.0

a

m 0.000

L

CL

90.0-r

0.00 1.00 2.00 3.00 4.00 5.00 6.00 7.00 8.00

5000.0

4500.0

4000.0

0) 3500.0

'a :) 3000.0

4j

-4 2500.0

c

m 

CP 2000.0

1500.0

1000.0

500.00

0.0000

0.00 1.00 2.00 3.00 4.00 3.00 6.00 7.00 8.00

Frequency (kr)

Figure 5.4 - Analytical Response For The Pressure At The

Wall Of A Pulsating Spherical Cavity
all 

83

Figure

analytical a

range studied

numerical r

approach zerc

approach inf

reason for tt

fied. Howeve

drawn upon tc

such as CoE

presented evi

with interior

mulations to

indirect bot

from a nonexi

those eigenf

mentioned in

utilized by

direct bounda

eigenfrequenc

not prop*en,

resolve the

5. 5 , an appr c

research. Thf

` in a very lov

F most problemF

.5 shows excellent correlation between

d numerical results throughout the frequency

except in the very low frequency regime. The

suits from kr - 0.25 to kr - C tend to

whereas the analytical results tent: to

nity within that same region. A definitive

s phenomenon has not been clearly identi-

, work published by other researchers may be

formulate possible explanations. Researchers

3y [17] and Burton and Miller [18] have

once that the coincidence of frequencies

,-:-igenfrequencies cause boundary element foriil for acoustic radiation problems. For the

iary element methods, the difficulty stems

dente of surface source distrihutions at

3quencies as shown by Copley [17]. As was

napter 2, overdetermination procedures •+ere

henck [19] and Paiszczyk and Kloaner [20] in

y element methods to overcome the interior

difficulties. Although it is possible, but

that an overdetermination proeradure may

low frequency diecrepanciets shown in figure

ch of this type was not used in the current

primary reason is that the phenomenon occurs

frequency regime which is not of interest in
180 .0 --

4J

M 0.000 -

Alj^

L

CL -180.0

lolswm

0.00 1.00 2.00 3.00 4.00 5.00 6.00 7.00 8.00

5000.0

4500.0 j

4000.0

Oi- 3500.0

7

3000.0

.^ 1

-+ 2500.0

C

0' 2000.0 t

1500.0

1000.0

500.00

0.0000

0.00 1.00 2.00 3.00 4.00 5.00 6.00 7.00 8.00

Frequency ( k r )

are 5.5 - Numerical And Analytical Response For The

x. Pressure At The Wall Of A Pulsating Spherical

Cavity

a
To further investigate the spherical cavity response,

pressure distributions within the cavity were numerically

determined at three frequencies. The three frequencies are

kr o - 2, 5.4 and 7.1 corresponding to 1372, 3704 and 4871

Hz, respectively if r o - Z meters and c - 343 m/s. The

spherical cavity pressure distributions for each of these

frequencies are shown in figures 5.6, 5.7 and 5.8. All

three distributions exhibit the expected spherical

Bessel's function-type behavior. The particular frequencies were selected to show the behaviors occuring before

and after the zero pressure frequencies, i.e.,

kr - A and 29 . The analytical solutions shown in these

figures are from equation 5.10. As in figures 5.2 and 5.3,

excellent correlation is shown in figures 5.6 through 5.8.
86

1@0.0 N

tn

m 0.000 L

0- -180.0

0.00 .030 .100 .150 .200 .250 .300 .330 .400 .450 .500

1000.0

800.00

j 600.00

C

400.00

E

200.00

0.0000

0.00 .050 .100 .150 .200 .230 .300 .350 .400 .450 .500

Radial Position

Figure 5.6 - Spherical Cavity Pressure Disttibution For

'f Velocity Boundary Conditions At kr-2

I
87

150.0

d

m o. 000

-180.0

0.00 .050 .100 .150 -i:.i' .2-50 .300 .350 .400 .450 .500

3000.0

2500.0

0U 2000.0

D

4j

•.+ 1500.0

C

01

m

E 1000.0

500.00

0.0000

0.00 .050 .100 .150 .200 .250 .300 .350 .4U0 .450 .500

Radial Position

Figure 5.7 - Spherical Cavit} Pressure Distribution For

Velocity Boundary Conditions At kr-5.4
88

180.0

VI

N 0.000

L

OL

-180.0

0.00 .050 .100 .150 .200 .250 .300 .350 ."0 .4: .500

IV*

6000.0

5000.0

W 4000.0

U

7

•• 3000.0

C

0%

N

E: 2000.0

1000.0

0.0000

0.00 .050 .100 .150 .200 .250 .300 .350 .400 .45C .500

Radial Position

Figure 5.8 - Spherical Cavity Pressure Distribution !,r

Velocity Boundary Conaitions At kr=7.1

{

k

4

i
89

1 Although the results presented thus far have only

been for pressure and velocity boundary conditions,

impedance boundary conditions have also been implemented.

Models having only impedance boundary conditions must contain an applied source within the cavity. Otherwise, the at

vector in equation 4.14 is a zero vector (see equation

4.17d) and the resulting source distribution is trivial.

In other words, there is no excitation of any sort in the

problem. To verify the impedance boundary condition implementation, a spherical cavity with purely impedance boundary conditions and an applied source of strength 1

located at the sphere center was utilized. The impedance

boundary conditions were established such that they would

simulate the impedance experienced by a point source radiating into an infinite free space. The impedance relationship can be derived from equation 5.6 as

Z = V - Pc(5.13)

(1 - k )

With a wave number of k - 0.583 (200 Hz) and r - 2 at- the

boundary, the resulting impedance is z - 32 .5 + j111.5.

Figure 5. 9 is the pressure distribution obtained numerically for the described model.
90

r

t

:e

180.0

ar

M 0.000 L

(_

-180.0

0.00 .050 .100 .150 .200 .250 .300 .350 .40f) .450 .500

;)uuu.0

1

1

4000.0

ar ^

3000.0

C

2000.0

E

1000.0

0.0000

0.00 .06 .100 .150 .200 .250 .300 .350 .400 .450 .500

Radial Position

Figure 5.9 - Pressure Distribution For A Sphere With Radiation Impedance Boundary Conditions At

k-0.583

4
91

For comparative purposes, similar distributions for

models with pressure and velocity boundary conditions were

obtained. The pressure solution and, consequently, the

boundary conditions may be determined in a number of ways.

In this case, the pressure boundary conditions were determined from equation 3.13 and 3.15.

P( x ) - ja(C)p (t,x) + E 4nt i p ( xsi ,x) ( 5.14)

B i-1

The pressure at any point is due entirely to he applied

source if the-boundary simulates an infinite domain problem. Hence, the source density distribution is zero at

the m-del boundary (i.e., o(C) - 0). Equation 5.14 becomes

e- j kr

p(r) - -JwP4n'. 1 r (5.15)

after the fundamental pressure solution is insertf . a. With

an applied source strength of 1, the equation to determine

the pressure boundary condition is

e- j kr

p(r) - - jwp4m r (5.16)

For the model with k - 0.583 and r - 2 .at the boundary,

the pressure boundary condition is p - -139 . 07 - j463.51.

The numerical results obtained are shown in figure 5.10.

Equivalent velocity boundary conditions found using equation 5.6 axe v - -4.167 + j0.033. Figure 5.11 shows the

pressure distribution numerically obtained using velocity

boundary conditions.
92

{

4

t

180.0 . 

m 0.000 L

CL #

-180.0

0.00 .050 .100 .150 .200 .250 .300 .330 .400 .430 .500

Z)uuu .0

4000.0

o,

^ 3000.0

c

2000.0

1000.0

0.0000

0.00 .050 .100 .150 .200 .250 .300 .330 .400 .450 .500

Radial Position
93

190.0

t

m 0.000

-190.0

0.00 .050 .100 .150 .200 .250 .300 .350 .400 .450 .800

8000.0

1

1

4000.0

1

10a

D 3000.0

C

2000.0

E

1000.0

0.0000

0.00 .080 .100 .150 .200 .250 .300 .350 .400 .450 .500

Radial Position

a

Figure 5.11 - Pressure Distribution For A Sphere With

Radiation Velocity Boundary Conditions At

k-0.583
94

To evaluate the validity of the impedance boundary

condition capability of the program, figure 5.9 may be

compared with figures 5.10 and 5.11. Both pressure and

velocity boundary conditions have previously been shown to

produce excellent resulte in the spherical cavity model.

Figures 5.9, 5.10 and 5.11 are identical pressure distributions. Regardless of the type of boundary , condition

being applied, the numerical results for this model are

unchanged.

The final study using the spherical cavity was to

numerically generate a forced response prediction for the

case of a point source in a rigid-walled cavity. velocity

boundary conditions of zero were applied to the sphere. A

source having a strength of 1 was placed at the center of

the cavity to excite the model. The swept frequency band

ranged from 0 to 873 Hz (or kr from 0 to 8) as was the

case for the previous response in figure 5.5. The response

for the rigid-walled spherical cavity is ehown in figure

5.12. In addition, a pressure distribution within the cavity for a frequency of 3087 Hz (kr - 4.5), near the first

resonant frequency, is shown in figure 5.13.
1

i

95

t

180.0

G

1n

.0 0.000

1!

-180.0

CL

0.00 1.00 2.00 3.00 4.00 5.00 6.00 7.00 8.00

1.5000

1.3500

1.2000

0 

1.0500

5.4

X

..^ .90000

.75000

60000

C

.45000

E .30000

15000

.00000 T---^

0.00 1.00 2.00 3.00 4.00 5.00 6.00 7.00 8.00

Frequency (kr)

Figure 5.12 - Numerical Response For The Pressure At The

Wall Of A Rigid-Walled Spherical Cavity

E

V,

- - --''
180.0

0.000

L

D` -180.0

0.00 .050 .100 .150 c^00 .230 .300 .330 .400 .450 .590
97

The rigid-walled spherical cavity response compares

.t

well with the response in figure 5.5 in terms of the

resonant frequency irlentification. Both responses identify

peaks occuring at approximately 3080 Hz and 5300 Hz (kr's

of 4.493 and 7.725). For the rigid-walled cavity, the

pressure magnitudes are much higher due to the effect of

t the applied source. The phase relationships differ between

the two responses for the same reason. The cavity pressure

distribution on figure 5.13 exhibits the expected Bessel's

function-type behavior. However, the effect of the applied

source can be seen near the radial position of zero. Since

the pressure at the applied source is infinite, the nearfield effect close to a radius of zero is expected. Otherwise, the pressure distribution is dominated by the

response of the cavity near the resonant frequency.. This

behavior could be expected to be much different if the

cavity was being excited at a frequency away from a

resonant frequency. Note that the low frequency regime in

figure 5.12 is charact'sristic of the low frequency regime

i

in figure 5.5. Once again, the low frequency results are

affected by the formulation and do not approach infinity

as they should.

In general, the indirect boundary element method has

done an excellent job for the variety of spherical cavity

I

problems tested. The models for boundary excitation,

internal point source excitations and with impedance

r
1.._1\ 98 3

IF,

boundary conditions all perform very well with the one

exception at low frequency.

5.2 Spherical Model For Radiation

Although the research being presented deals primarily

with cavity acoustics, it is possible to deal with acoustical radiation problems applying the same formulation. To

do so, the models need only have negated 'free terms'

(i.e., the quantity c in equations 4.16b and 4.17b). In

other words, instead of having the 'free terms' be negative (as is done for cavity problems), the 'free terms'

must be positive.

for the sphere was

the formulation would

as for the exterior

the previous section,

been demonstrated by

The acoustical radiation model

studied only to determine whether

fail at the interior eigenfrequenci,

radiation problems. As mentioned in

this type of formulation failure has

other researchers.

Radiation from a pulsating sphere can be modeled by

applying uniform velocity boundary conditions to the

sphere and applying the 'free term' negation. The analytical solution for the pressure on the surface of a pulsating sphere is (52)

kr 0

p(ro) - PcVo

2 (kro + (5.16)

((kr 0 ) + 1)

ow
99

Equation 5.16 was utilized to generate the analytical

boundary pressure response in figure 5.14 for a velocity

boundary condition of 1.

The sphere model was used to numerically determine a

boundary pressure response plot for a pulsating sphere. A

uniform velocity boundary condition of 1 was applied to

the sphere. The 'free terms' were negative of those used

for the spherical cavity studies. Figure 5.15 contains the

analytical solution represented by the dashed line and the

numerical solution represented by the solid line marked

.` with Y's at discrete frequency points.
100

x

180.0

aJ

m 0.000

.0

Q_ -180.0

0.00 1.00 2.00 3.00 4.00 5.00 6.00 7.00 8.00

1000.0

900.00

800.00

700.00

600.00

-• 500.00

C

400.00

E 300.00

200.00

100.00

0.0000

0.00 1.00 2.00 3.00 4.00 5.00 6.00 7.00 8.00

Frequency (kr)

i 

Figure 5.14 - Analytical Pulsating Sphere Boundary Pressure Prediction For Radiation
1J1

AAl1

180.0

LAVI

(0 0.000

.c

a_ 

-180.0

0.00 1.00 2.00 3.00 4.00 5.00 6.00 7.00 8.00

1000.0

900.00

800.00

700.00

600.00

500.00

C

400 .00

ra

300.00

200.00

100.00

o.00bo

0.00 1.00 2.00 3.00 4.00 5.00 6.00 7.00 8.00

Frequency Ckr>

Figure 5.15 - Numerical And Analytical Pulsating Sphere

Boundary Pressure Predictions For Radiation

s

1
3

i

f

102

4

Figure 5.15 very clearly shows that the formulation

.. failures for the radiation problem are present in the formulation being applied in this research. At the interior

eigenfrequencies of kr equal to n and 2n, the boundary

_ pressures diverge from the analytical solution. Koopmann

and Benner [29] present results for a pulsating sphere

which portray the same behavior as shown in figure 5.15.

Thus, the conclusion to be drawn is that an overdetermination procedure should be applied if this formulation were

to be used for radiation problems. Also, these results

lend support to the hypothesis that the difficulties

obtained in the low frequency regime for the cavity problems may be due to a formulation failure. Near a frequency

of zero, rigid-body modes exist which this formulation is

evidently incapable of determining for cavity problems.

5.3 Rectangular Cavity Model

All the results presented heretofore in this chapter

have been for a spherical model. For boundary element

methods, corners and edges are known to create difficulties for boundary element techniques due to the boundary

discontinuities. These difficulties are described in the

literature in some detail. Banerjee and Butterfield [12]

have devoted an entire chapter in their text to the subject. Although the problem is well-documented, there does

not exist one ultimate solution for all cases.
103

F To investigate the consequences of corners and edges,

a rectangular cavity with dimensions of 4x5x3 was studied.

Three different methods of modeling the corners and edges

were attempted. The first method was to model the corners

and edges as intersections of adjacent sides. Essentially, the first method is the intuitive aproach one

would assume having no knowledge of the inherent difficulties. The normal vectors for the nodes on the corners and

edges were given components perpendicular to -each intersecting wall. The second method involved leaving gaps

between the intersecting walls. In essence, the nodes at

the corners and edges were not physically attached to

adjacent walls of the box. Hence, compatibility of adjacent node points is not enforced. (Element formulations of

this type are referred to as 'noncompatible' elements.)

The final method involved rounding off the corners and

edges to avoid the boundary discontinuities.

To evaluate the three methods of modeling the corners

and edges, one wall of the cavity was driven while the

remaining five walls were rigid. A simplistic analytical

approach for this situation is to consider the problem on

a one-dimensional basis as a driven/closed cavity. The

pressure solution for this problem is

p Ae jk(L - x) + Be-jk(L - x) (5.17)

where L is the length of the cavity in the direction
104

normal to the driven end, x is the length variable ranging

from 0 to L and the constants A and B are determined by

the boundary conditions. The boundary conditions are that

the velocity at the closed end is zero (V(L) - 0) and the

velocity at the driven end is equal to the driving velocity (V(0) - Vo ). With these boundary conditions and equa- 

tion 5.17, the approximate analytical pressure solution is

cosk(L - x)

P - -jpuv sin(kL) (5.18)

To quantify equation 5.18, a driving velocity of 1 and a

dimension of L 5 (corresponding to the y dimension of

the cavity) were used. The pressure distribution along

the y coordinate for wave numbers of k equal to 1, 2 and 3

(or frequencies of 343, 686 and 1092 Hz with a speed of

sound of 343 m/s) are shown in figures 5.16, 5.17 and

i

5.18, respectively. 
r

i

W

L

a.

105

a

180.0

0.000

-180.0

0.00 .500 1.00 1.50 2.00 2.50 3.00 3.50 4.00 4.30 5.00

300.00

450.00

400.00

Oa 350.00

j 300.00

250.00 C

200.00

E 150.00

100.00

50.000

0.0000

0.00 .500 1.00 1.30 2.00 2.50 3.00 3.50 4.00 4.50 5.00

Y-Coordinate

Figure 5.16 - Analytical Pressure Distribution For A

Driven/Closed Cavity At k=1
106

malytical Pressure Distribution For A

► riven/Closed Cavity At k-2

180.0 d

m 0.000 L

CL -180.0

0.00 .500 1.00 1.50 2.00 2.50 3.00 3.30 4.00 4.50 5.00

1000.0

900.00

800.00

700.00

j 600.00

500.00 C

400.00

E 300.00

200.00

100.00

0.0000

0.00 .500 1.00 1.50 2.00 2.50 3.00 3.50 4.00 4.50 5.00

Y-Coor d i nat e
107

m 0.000

c

0.00 .500 1.00 1.50 2.00 2.50 3.00 3.50 4.00 4.30 5.60

900.00 -

700.00 -

600.00 -

'0 500.00

400.00

c

cp

M 300.00

200.00

100 - Oc

0.0000

0.00 .500 1.00 1.50 2.00 2.50 3.00 3.50 4.00 4.50 5.00

Y-Coordinate

Figure 5.18 - Analyt'cal Pressure Distribution For A

Driven/Closed Cx^rity At k-3
Boundary element models incorporating the

W 

corner and edge methods were used to numerically evaluate

a driven/closed cavity. However, the driving velocity was

not uniform as assumed in the analytical solution. An

attempt to produce such a situation is impossible in the

first and third models because the normal velocity at the

corner nodes must be compatible on both walls. Imposing

such a boundary condition would cause the adjacent walls

to be non-rigid in the first and third models. Consequently, a velocity profile having a magnitude of 1 at the

center but varying linearly to zero at the edges was

applied to the driving wall (i.e., the wall at y equal to

zero). Pressure distributions were calculated for wave

numbers of k equal to 1, 2 and 3 just as for the analytical solution. The numerical results are shown in figures

5.19, 5.20 and 5.21 for the first model (the true box),

figures 5.22, 5.23 and 5.24 for the second model (the

gapped box) and figures 5.25, 5.26 and 5.27 for the third

model (the rounded box).

i
109

180.0

t0 0.000

L

CL -180.0

0.00 500 1.00 1.50 2.00 2.50 3.00 3.50 4.00 4.50 5.00

500.00

450.00

400.00 

350.00 

300.00

250.00

c

c 200.00

E: 150.00

100.00

50.000

0.0000

0.00 .500 1.00 1.50 2.00 2.50 3.00 3.50 4.00 4.50 5.00

Y-Coordinate

Figure 5.19 - Pressure Distribution For First Box Model At

k-1

i 
110

A

180.0

W

(0 0.000

L

CL _ion. o

0.00 .500 1.00 1.50 2.00 2.50 3.00 3.50 4.00 4.50 5.00

iuuv.v

900.00

800.00

0) T00.00

600.00

500.00

C

400.00

E 300.00

200.00

100.00

0.000f)

0.00 .500 1.00 1.50 2.00 2.50 3.00 3.50 4.00 4.50 5.00

Y-Coordinate

Figure 5.20 - Pressure Distribution For First Box Model At

k-2

t

I
180.0

(b

LA

m 0.000

.c

a.

-180.0

lift -0

ill I 

0.00 . 500 1.00 1.50 2.00 2-50 3.00 3.50 4.00 4.50 5.00

800.00

700.00-

600.00

500.00

400.00

C

300.00

200.00

100.00,

0.0000

0.00 510 1.60 1-50 2-00 2.50 3.00 3.50 4.00 4.50 5.60

Y-Coordinate

Figure 5.21 - Pressure Distribution For First Box Model At

k-3
117.

4

1

180.0

w

U

r0 0.000

L

OL -180.0

0.00 .500 1.00 1.50 2.00 2.50 3.00 3.50 4.00 4.50 5.00

zIuu.uu

450.00

400.00

N 350.00

300.00 t

-^ 250.00

C

U 200.00

150.00 s

i

1G0.00

50.000

0.0000

0.00 .500 1.00 1.50 2.00 2.50 3.00 3.50 4.00 4.50 5.00

Y-Coordinate

Figure 5.22 - Pressure Distribution For Second Box Model

At k-1

a.
iuuu . u

900.00

800.00

N 700.00

7

600.00

-. 500.00

C

^0 400.00

E: 300.00

200.00

100.00

0.0000

0.00 .500 1.00 1.50 2.00 2.50 3.00 3.50 4.00 4.50 5.0(

Y-Coordinate

bution For Second Box Model
114

(a 0.000

-CCL 100.0

0.00 .500 1.00 1.50 2.00 2.50 3.00 3.50 4.00 4.50 5.0 0

,tUUU . U -

3500.0 -

3000.0

as

,a 2500.0

2000.0 -

C

rp 1500.0 -

1000.0 -

500.00

0.0000

0.04 .500 1.00 1-ZQ 2.00 2.50 3.00 3.50 4.60 4.50 5.00

Y-Coordinate

Figure 5.24 - Pressure Distribution For Second Box Model

At k-3
115

180.0

Ln

m 0.000

L

CL -180.0

0.00 .500 1.00 1.50 2.00 2.50 3.00 3.50 4.00 4.50 5.00

5

500.00

450.00

400.00

350.00

7 300.00

N

•+ 250.00

C

200.00

1: 150.00

100.00

50.000

0.0000

0.00 .500 1.00 1.50 2.00 2.50 3.00 3.50 4.00 4.50 5.00

Y-Coordinate

Figure 5.25 - Pressure Distribution For Thi=d Box Model At
116

180.0

m 0.000

L

a 

-180.0

0.00 .500 1.00 1.50 2.00 2.50 3.00 3.50 4.00 4.50 5.00

1000.0

900.00

800.00-

al 700.00

600.00

•-+ 500.00

C

C 400.00

Z: 300 .00

200.00

100.00

0.0000

0.00 .500 1.00 1.50 2.00 2.50 3.00 3.50 4.00 4.50 5.00

Y-Coordinate

Figure 5.26 - Pressure Distribution For Third Box Model At

k-2

I
A. 180.0

t r W

i

^O 0.000

1

a -180.0

0.00 500 1.00 1.50 2.00 2.50 3.00 3.50 4.00 4.50 5.00

Si

200.00

100.010

0.0000

0.00 .500 1.00 1.50 2.00 2.50 3.00 3.50 4.00 4.50 5.00

800.00

700.00

600.00

OJ

500.00

-. 400.00

C

r^0 300.00

I

)x Model At

f
118

Two points must be made regarding the numerical

models. Firstly, the excitation at the driving wall is not

uniform as in the analytical solution. Secondly, the

numerical models are three-dimensional. Thus, one can not

expect purely one-dimensional results, especially when the

excitation at the driving wall is not purely ones

dimensional. These points need be highlighted to provide

explanations for the differences between the analytical

and numerical solutions.

z 

.. With the previous points in mind, some qualitative

judgements can be made on the three numerical models. The

analytical solutions in figures 5.16 through 5.18 can be

used as guidelines to evaluate the three models. For all

r

three frequencies, the analytical solutions are

cosinusoidal in nature. The numerical results for the

first model appear to exhibit some form of sinusoidal

behavior for the second and third wave numbers but not as

clearly defined as would be expected either in magnitude

r^ 

or phase. The second model shows more clearly cosinusoidal

behavior in the pressure distribution. Phase shifts of

9e

o

180 are evident at positions of zero pressure magnitude.

There is some variability in the peak magnitudes within an

f

individual pressure distribution. However, it is quite

possible that the variability is due to cross-coupling

with modes in the other two orthogonal directions since

the absolute maximums occur near the center of the cavity
119.

(e.g., figure 5.24). Finally, the third model, the

rounded corner model, shows cosinusoidal behavior at the

first wave number but not at the second or third wave

numbera. The pressure magnitudes are not regular peaks nor

are any clearly defined phase shifts evident. As a result

of this numerical study, the second method of modeling the

Fs

corners and edges, that of providing gaps, appears to be

the best of the three methods and is also the one often

E

recommended in the literature (e.g., [12]). It should

also be noted that another solution which should be investigated for acoustical problems is the use of 'noncompatible' elements. In 'noncompatible' elements, the node

r;. points lie in the interior of the element rather than on

the edges. Many of the problems encountered here will not

be a consideration with 'noncompatible' elements.

After determining that the gapped edge model provided

the beat results of the three models studied, the gapped

edge model was utilized to study the internal source excitation problem of the rectangular cavity. The boundary

conditions for the response study were rigid-wall boundary

conditions. Excitation of the model was supplied by a

source placed in the center of the rectangular cavity. The

strength of the source was 1. For a rigid-walled rectangular cavity, the pressure distribution in terms of

modal response will be of the form (53)

f

I,
120

A x ) p lmn cos( Lx )cos(L^)cos(LZZ

(5.19)

1 m n Y

where Lx , L and L are the dimensions of the rectangular

cavity in the x, y an6 z directions, respectively and Alm

is the modal amplitude. The modal amplitude will depend on

the closeness of the driving frequency to a natural frequency and also on the coupling of the exciter to the

mode. Since an applied source is placed at the center of

the cavity for excitation purposes, the odd modes will not

be excited. (Odd modes means that the variables 1, m or n 

in equation 5.19 are odd numbers.) The applied source

creates pressure antinodes in the center of the cavity.

Hence, only even modes should be excited.

The numerically determined response plots for the

rigid-walled :ectangular cavity are shown in figures 5.28,

5.29 and 5.30. The pressures shown in the response plots

are associated with a point located in the center of a

wall of Cie cavity. For Example, the x direction response

is for a point located in the center of a wall having a a

constant x coordinate (i.e., a wall parallel to the yz

plane). The response plots indicate resonant frequencies

by the significant peaks in the magnitudes and shifts in

the phase. Incidentally, the difficulties in the low frequency regime are apparent again in these response plots

below a wave number of about 0.25.
AA,,

180.0

'W

'w

IA

(0 0.000

.0

0` -180.0

0.00 .250 .500 .750 1.00 1.25 1.50 1.75 2.00 2.25 2.50 

8000.0

7000.0

6000.0

ap

5000.0

7

.+ 4000.0

C

CP

m 3000.0

2000.0

1000.0

0.0000

0.00 .250 .500 .750 1.00 1.25 1.50 1.75 2.00 2.25 2.50

Frequency (k)

Figure 5.28 - Rectangular Cavity Pressure Response - X

Direction
W

180.0 

m 0.000

.0

a -180.0

0.00 .250 .500 .730 1.00 1.25 1.50 1.75 2.00 2.2w 2.50

8000.0

7000.0

6000.0

0)

5000.0

7

.+ 4000.0

c

r^0 3000.0

2000.0

1000.0

0.0000

0.00 .250 .500 .750 1.00 1.25 1.50 1.75 2.00 2.25 2.50

Frequency C k >

Figure 5.29 - Rectangular Cavity Pressure Response - Y

Direction
a,

180.0

0.000

.0

OL -180.0

0.00 .250 .500 .730 1.00 1.23 1.30 1.75 2.00 2.23 2.50

8000.0

7000.0

6000.0

"0 5000.0

40?0.0

C

r^0 3000.0 E

2000.0

1000.4

0.0000

0.00 .250 .500 .750 1.00 1.25 1.50 1.75 2.00 2.25 2.50

Frequency Ck)

Figure 5.30 - Rectangular Cavity Pressure Response - Z

Direction

.123
124

Pressure distributions within the rectangular cavity

were generated from the numerical results to identify the

pressure response shapes occuring at rssonant frequencies.

Near the resonant frequencies the response should be dominated by the respective mode shape. Between resonance, it

would be expected that the response would be lower and

would be a superposition of all the near-frequency

resonant mode shapes. Since the cavity and pressure distributions are three-dimensional, the pressure distributions were determined for planes perpendicular to one of

the orthogonal directions and passing through the center

of the rectangular cavity. Thus, for each response at

which a distribution is determined, three plots are generated; one in each of the three orthogonal planes. Furthermore, the magnitude at the center of all the magnitude

plots is infinite and the phase at the center of the phase

plots is -90 0 because of the applied source. As a result,

the response of the center point is truncated to the level

of the maximum pressure magnitude.

From the cavity response, the first resonant frequency occurs at a wave number of 1.34 (460 Hz). The pressure distribution within the rectangular cavity was generated at that wave number. The pressure distribution of

the variations in the z direction is shown in figure 5.31.

Likewise, figure 5.32 is the y variations and figure 5.33

is the x variations. Figure 5.31 indicates a possible

s

j

f

i
125

cosinusoidal variation in magnitude along the y coordinate

which is confirmed by the clearly defined 180 0 phase shift

along the y coordinate. Relatively constant pressures are

displayed in figure 5.32. The cosinusoidal variation in

the y coordinate is evident again in figure 5.33. It is

readily apparent that the (0,2,0) mode is the resonant

mode shape at 460 Hz. Analytically, the (0,2,0) mode shape

occurs at k = 5n or 432 Hz. The conclusion is that the

boundary element model successfully located the first

resonant frequency with about: a 6% upward shift in the

frequency.

.z:

M 
126

180.0

90.00

d

L

0.000

n.

-90.00

-180.0

5

k..Coos

1narP

5

1.600

CD 1.400

X 1.200

1.000

.8000

.6000

r+ .4000

C

.2000

E .0000

LQ

O^ G°

U

Figure 5.31 - Rectangular Cavity Pressure Distribution At

k-1.34 (XY Plane)

F
180.0

90.00

a)

N

m 0.000

.0

d

-90 .00

-180.0

u

3

t

• 

3

^. 1.600

CD 1.400

X 1.200

v 1.000

.8000

.6000

.4000

C

0 .20000

E .0000

4 p a.e

ICporcl1

aa< dl

rya tP ^,G

Figure-5.32 - Rectangular Cavity Pressure Distribution At

k=1.34 (XZ Plane)
5

180.0

90.00

QJ

L

0 .000

a-

-90.00

-180.0

u

^. 1.600 0 1.400

X 1.200

`e 1.000

U

.8000

.6000

.4000 C

CT .2000

^ .0000

5

1

V U

Coord

Inar P

Ot Go

Figure 5.33 - Rectangular Cavity Pressure Distribution At

k-1.34 (ZY Plane)
129 

The second resonance condition occurs at a wave

•' s

number of about 2.65 ( 566 Hz). The pressure distributions

within the rectangular cavity are given' in figures 5.34

through 5.36. A cosinusoidal variation along the x direction is indicated by the magnitudes and by the 180 0

 phase

shifts in figures 5.34 and 5 . 36. Hence, the (2,0,0) mode

shape appears to be the dominant mode at 566 Hz. The

(2,0,0) mode occurs analytically at k - 4n or 539 Hz. For

the second resonant condition, the boundary element model

located the resonance with a 5% upward shift in the frequency.

s

i

P.

1

• 3

:t

 e 

f

4

k

{[
180.0

90 .00

m 0.000

.0

(L

-90.00

-180.0

U

4.800

4.200

X 3.600

3.000

2.400

D 1.800

1.200

C

.6000

m

Z: .0000

3^^ 10 ,,e

CO 0(

0^-

Figure 5.34 - Rectangular Cavity Pressure Distribution At

k-1.65 (XY Plane)

L
180.0

z

90.00

aJ

r0 0.000 L

OL

r

-90.00

• 3

3

4.800

p 4.200

X 3.600

.r

3.000

2.400

1.800

•+ 1.200

C

.6000

Z- .0000

u

r 

4 -w 0 a^, e

t

Opf„d1 L

nd t

/GpO d

P 2

Figure 5.35 - Rectangular Cavity Pressure Distribution At

k-1.65 (XZ Plane)

Ed

M'
180.0

90.00

0) W

r

0.000

0L

-90.00

-160.0

5

Coot-

-

5

4.800 0 4.200

9-4

X 3.600

3.000

40 2.400

1.800

1.200 C

CP

N

.6000

E .0000

t,e

< a1% o

G°

132

U

Figure 5.36 - Rectangular Cavity Pressure Distribution At

k-1.65 (ZY Plane)
133

The third resonant pressure distribution is at a wave

number of 2.06 (707 Hz). In the response of figures 5.28

through 5.30, two closely-spaced peaks appear between the

wave numbers of 2 and 2.25. The first peak is the one at

707 Hz. Figures 5.37 through 5.39 are the pressure distributions of the third resonant condition. In figure 5.37,

it seems as though a cosinusoidal variation exists along

both the x and y directions. However, it is not entirely

clear from the phase relationship that this is indeed the 

case. By examining figure 5.38, it becomes more obvious

that just such a variation is present in the x direction.

Similarly, the y direction variation is better viewed in

figure 5.39. The resulting conclusion is that the (2,2,0)

mode is the dominate mode at this frequency. Theoretically, the (2,2,0) mode occurs at a wave number of

(4 2 + (5n) 2 - 2.01 (689 Hz). For the third resonance

in the rectangular cavity response, the boundary element

model located the resonant frequency with only a 3% upward

shift. 
I

134

+00

ly

,

u

,COOrd

1^arP

5

C) 0.100

X 1.800

`* 1.500

Ib 1.200

.9000

C .6000

O' .3000

.0000

t,e

GO 0t

• u

5

0.0

.00

000

.00

D.0

Figure 5.37 - Rectangular Cavity Pressure Distribution At

k-2.06 ( XY Plane)
3

U

p 2.100

X 1.800

1.500

1.200

.9000

-^ .6000

C

.3000

.0000

3

180.0

90.00

L

0.000

d

-90.00

-180.0

Li

0 tie ^a

^Cooraj n ooh d1

ar 0 2.^

Figure 5.38 - Rectangular Cavity Pressure Distribution At

k-2.06 (XZ Plane)
5

180.0

90.00

ar

t

0.000

a.

-90.00

-180.0

^i

U

^. 2.400

0 2.100

^{ 1.800

1.500

1.200

.9000

C

.6000

.3000

.0000

Copra

jnar

5

,,e

Lo

1

Figure 5.39 - Rectangular Cavity Pressure Distribution At

k-2.06 (ZY Plane)

a

-- -
137

Although the pressure distribution will not be

presented for the resonant peak at the wave number of 2.21

(758 Hz), some comments can be made regarding it. The next

expected mode shape should be the (0,0,2) mode which,

theoretically, is excited at a wava number of 2.09 (717

Hz). The numerically determined pressure distribution at

758 Hz was characteristic of the (0,0,2) mode but also

retained the nature of the (2,2,0) mode which occured at a

slightly lower frequency. Because the two resonances are

so closely related in terms of frequency, both the (2,2,0)

mode and the (0,0,2) mode are apparent in the distribution

at 758 Hz. A situation such as this is likely to happen

whenever the modal density is high in a particular frequency range cf interest.

A series of pressure distributions were gen^rated for

a non-resonant frequency of the rigid-walled rectangular

cavity. The distributions in figures 5.40 through 5.42 are

for a wave number of 1 (343 Hz) which is a frequency lower

than the first resonant frequency of 460 Hz. The pressure

distributions plainly show that at this non-resonant frequency, the applied source is the sole determinant of the

pressure field characteristics. This type of distribution

should be expected as all the resonant response will be

t

low at this frequency.
. U

. 0

00

00

0o

. 0

s

`Copra►

1nate

3

e% 1.600 c)

Q 1.400

X 1.200

1.000

.8000

.6000

-+ .4000

C

.2000

.0000

,,e

W-1 Gppt

Figure 5.40 - Rectangular Cavity Pressure Distribution At

k-1 (XY Plane)
180.0

90.00

4)

V7

L

0.000

a-

-90.00

-180.0

3

^. 1.600

q 1.400

X 1.200

1.000

.8000

D .6000

.4000 C

.2000

E .0000

3

u

4 0 atie

'Coordjn ooh d^

atP 2.G

Figure 5.41 - Rectangular Cavity Pressure Distribution At

k-1 (XZ Plane)
5

3

4

3

1.600

1.400

x 1.200

%..* 1.000

.8000

.6000

.4000

.2000

.0000

tie

l igure 5.42 - Rectangular Cavity Pressure Distribution At

k-I (ZY Plane)
.1

a:

141

5.4 Summary

The first section of this chapter presented numerical

results and, analytical results where uL-)plicable, of sound

fields in a spherical cavity. Excellent correlation was

obtained except in the low frequency regime. A hypothesis

for the low frequency difficulties was presented and partially supported through literature references. Using the

spherical cavity model, all three boundary condition+ types

(pressure, velocity and impedance) and the applied source

capabilities were verified. in addition, the sphere model

was used to briefly examine the acoustical radiation prediction capabilities of the program for a pulsating sphere.

The remainder of the chapter presented results

obtained for a 4x5x3 rectangular cavity. Three different

techniques for modeling corners and edges were presented.

The best of the three modeling techniques was to model the

corners and edges by leaving slight gaps. A model of the

cavity having gaps at the corners and edges was utilized

to develop a cavity response and pressure distributions

both on and off the resonant frequencies. The resonant

frequencies were located and identified with good accuracy

using the boundary element model.
142

s • _

V

 ..

^ rA

CHAPTER 6

CONCLUSIONS

The objective for this research was to develop a

numerical tool based on indirect boundary element theory

to predict the sound field within an acoustical cavity. In

addition, two additional features were desired. The first

additional feature was that the three boundary condition

types; pressure, velocity and impedance, be formulated and

implemented in the development. Secondly, the implementation of acoustic point sources was desired.

A code to implement these objectives was developed

and verified. In general, the conclusions which can be

drawn from this experience are;

1. The indirect boundary element method is capable of all

analyses for which direct boundary element methods have

been demonstrated in the literature. In addition, IBEM

N

techniques require only one boundary integral in most

3 
I

s

x

ry

143

cases. This will be more efficient than DBEM's in many

instances.

2. Forced response in three-dimensional cavities is a

straightforward application of boundary element

methods. Corners and edges cause difficulties as they

do with radiation problems. Further conclusions in this

regard are included in the next paragraph.

3. Implementation of internal source capabilities is reasonably straightforward and potentially quite powerful

in future applications such as active noise control.

Boundary element methods are decidedly better then finite element methods for representing internal point

sources since BEM's can accommodate singularities.

4. The application of impedance boundary conditions was

successfully demonstrated. It should be pointed out

that the formulation used in this development is not

limited to locally reacting impedance.

One suggestion can be made regarding the problem with

corners and edges evidenced in the rectangular box model.

The application of & 'noncompatible' element may resolve

these difficulties. A 'noncompatible' element is one in

which the nodes of the element are not located along the

edges of the element as for the isoparametric element used

in this research but located internal to the element. By
locating the nodes, or the locations of the boundary conditions, within the element, the boundary conditions will

never be situated at a boundary discontinuity. Thus, the

'free terms' in the velocity and impedance boundary condition equations will always maintain a value of t2 and,

more importantly, the difficulties with the boundary

discontinuities can be avoided.

The boundary element procedure developed for this

research can be applied to noise source identification.

Although this concept was not specifically pursued in the

case studies, the identification of acoustical generators

and optimizat-.on of the cavity acoustical characteristics

can be guided with the application of this boundary element procedure.

In summary, it appears from this work that the

indirect boundary element method should be the numerical

method of choice for a large precentage of investigators

interested in forced response problems in both interior

and exterior domains. This is especially true when only a

limited number of frequencies and response locations are

desired.

*WOO
LIST OF REFERENCES

i

F

i
LIST OF REFERENCES

[lJ A.K. Noor and W.D. Pilkey, editors, State-of-The-Art

Surveys 2a Finite Element Technology, ASME, 1983.

(2] R.J. Bernhard, "Using Structural Dynamics Code To Do

Acoustic Finite Element Analysis", Herrick Laboratory Report 85-5, Purdue University, 1985.

(3] C.-I J. Young and M.J. Crocker, "Prediction of

Transmission Loss in Mufflers by the Finite Element

Method", J. ACous. Soc. of Am., Vol. 57, No. 1,

1975, pp. 144-148.

(4] C.-I J. Young and M.J. Crocker, "Finite Element

Acoustical. Analysis of Complex Muffler Systems With

and Without Wall Vibrations", Noise Control

Engineering, Sept.-Oct. 1977, pp. 86-93.

[5J R.J. Astley and W. Eversman, "A Finite Element Formulation of the Eigenvalue Problem in Lined Ducts

with Flow", J. of Sound and Vim., Vol. 65, No. 1,

1979, pp. 61-74.

(6] A. Craggs, "The Application of Acoustic and Absorption Finite Elements to Sound Fields in Small Enclosures", Finite Element Applications in Acoustics,

ASME, 1981, pp. 1-19.

(7) L. Cederfeldt, "On the Use of the Finite Element

Method on Some Acoustical Problems", ASME paper 81-

WA/NCA-1, 1981.

(8] S-F. Ling, J.P. Hamilton and J.J. Allen, "The Use of

Two Dimensional Isoparametric Finite Elements for
146

Determining the Sound Field in Axisymmetric Lined

Ducts", Finite Element Agolicat.icng in Acoustics,

ASME, 1981, pp. 99-110.

[9] G.R. Doyle, Jr. and L . L. Faulkner, * ThreeDimensional Finite Element Acoustic Pressure

Analysis with Complex Impedance", Finite Element

Applications In Acoustics, ASME, 1981, pp. 57-72.

[10] A. Craggs, "An Acoustic Finite Element Approach for

Studying Boundary Flexibility and Sound Transmission

Between Irregular Enclosures", J. 2L sound and Vib.,

Vol. 30, No. 3, 1973, pp. 343-357.

[+1] S.H. Sung, " Automotive Applications of ThreeDimensional Acoustic Finite Elements", SAE paper no.

810397, 1981.

[12] P.K. Banerjee and R. Butterfield, Boundary Element

Methods ja Engineering Science, McGraw-Hill, London,

1981, p.4.

(13] L.H. Chen and D . G. Schweiker. t, "Sound Radiation from

an Arbitrary Body", J. Acous. Soc. 2f LkM., Vol. 35,

No. 10, 1963, pp. 1626-1632.

(14) G. Chertock, " Sound Radiation from Vibrating -Surfaces", J. Acous. Soc. of Am., Vol. 36, No. 7, 1964,

pp. 130571313.

(15) L.G. Copley, "Integral Equation Method for Radiation

from Vibrating Bodies", J. Lkcous. Soc. 21 AM., Vol.

41, No. 4, 1967, pp. 807-816.

(16) W. Williams, N.G. Parke, D.A. Moran and C.H. Sherman, "Acoustic Radiation from a Finite; Cylinder", J.

Acous. Soc. of Am., Vol. 36, No. 12, 1964, pp. 

2316-2322.

(17] L.G. Copley, "Fundamental Results Concerning

Integral Representations in Acoustic Radiation", J.

Acous. Soc. of Ai., Vol. 44, No. 1, 1968, pp. 28-32.

[18] A.J. Burton and G.P. Miller, "The Application of

Integral Equation Methods to the Numerical Solution

of Some Exterior Boundary -Value Problems", Froc.

Roy. Soc. Lond., A.323, 1971, pp. 201-210.

(19) H . A. Schenck, " Improved Integral Formulation for

E Acoustic Radiation Problems", J. Acous. S oc- of Am.,

6 Vol. 44, No. 1, 1968, pp. 41-58.

3
also" w -awmvW­­^

U ,

147

(20] C.M. Piaszczyk and J.M. Klooner, "Acoustic Radiation

from Vibrating Surfaces at Characteristic Frequenciec*, ff.. Acous. &M., Vol. 75, No. 2, 1984,

pp 363-375.

(21] R.R. Smith, J.T. Hunt and D. Barach, "Finite Element

Analysis of Acoustically Radiating Structures with

Applications to Sonar Transducers", Acous. lu2f, &M., Vol. 54, No. 5, 1973, pp. 1277-1288.

[221 J.J. En9blom and R.B. Nelson, *Consistent Formulation of Sound Radiation From Arbitrary Sturcture",

,1. 21 Apjal. HIgh., June 1975, pp. 295-300.

[23] W.L. Meyer, W.A. Bell, B.T. Zinn and M.P. Stallybrass, O Boundary Intagal Solutions of Three Dimensional Acoustic Radiation Probleme n , J. Sound And

Ylk., Vol. 59, No. 2, 1978, pp. 245-262.

(241 A.F. Seybert, B. Soenarko, F.J. Rizzo and D.J.

Shippy, "Application of the BIE Method to Sound

Radiation Problems Using An Isoparametric Element",

ASME paper no. 82-WA/NCA-1, 1982.

(25) R. Seznec, "Diffraction of Sound Around Barriers:

Use of the Boundary Elements Technique", J. Sound

and Vib., Vol. 73, No. 2, 1980, pp. 195-209.

(261 D.T.I. Francis, M.M. Sadek and S.A. Tobias, "Prediction of Acoustic Emission of Impact Forming Machines

from Design Drawings", ASME paper no. 81-WA/NCA-5,

1981.

(27) S. Vajpayee, M.M. Nigm and M.M. Sadek, "Acoustic

Diagnosis of a Hydraulic Hammer", ASME paper no.

81-WA/NCA-8, 1981.

(28) H. Benner and G.H. Koopmann, "Sound Power Prediction

Using the Helmholtz-Kirchhoff Integral Equation",

ASME paper no. 81-WA/NCA-4, 1981.

(29 1 G.H.Koopmann and H. Benner, *Method for Computing

the Sound Power of Machines Based on the Helmholtz

Integral", 9. Acous. Soc. 21 AM., Vol, 71, No. 1,

1982, pp. 78-89.

[30) G.H. Koopmann and J.C. Perrand, "Crack Location by

Means of Surface Acoustic Intensity Measurements",

ASME paper no. 81-WA/NCA-6, 1981.

(31] G.H. Koopmann amd M. Mezache, "Radiation Characteristics of a Slender Box-Type Structure", ASME
ILI:

148

paper no. 82 -WA/NCA- 9, 1982.

(32] G.R . C. Tai and R . P. Shaw, " Helmholtz-Equation Eigenvalues and Eigenmodes for Arbitrary Domains", J.

Acous. g,Q,^. Q,t ^., Vol. 56, No. 3, 1974, pp. 796-

804.

[33] T. Tanaka, T. Fujikawa, T. Abe and H. Utsuno, "A

Method for the Analytical Prediction of Insertion

Loss of a Two -Dimensional Muffler Model Based on the

Transfer Matrix Derived From the Boundary Element

Method", ASME paper no. 84 -WA/NCA-7, 1984.

(34] A. Sestieri, D. Del Vescovo and P. Lucibello,

"Structural-Acoustic Coupling in Complex Shaped Cavities", ,gound Ind Vj^., Vol. 96, No. 2, 1984, pp.

219-233.

(35] C.A. Brebbia and S. Walker, Boundary Element Tech_

ngges ja Angineerina, Newnes -Butterworths, London,

1980.

(36] P.K. Banerjee and R. Butterfield, Developments in

Boundary Fjement Methods z. j. , Applied Science Publishers, London, 1979.

(37) L.E. Kinsler, A.R. Frey, A.B. Coppens and J.V.

Sanders, Fundamentals gL Acoustics, Third Edition,

John Wiley & Sons, New York, 1982, p.113.

(38) W.A. Kaplan, Advanc2d Mathematics I_qL Engin.ers,

Addison-Wesley, Reading, Mass., 1981, p.475.

(39) Op. cit., Ref. 35, pp. 33-34.

[40] On. =.,, Ref . 12, p. 122.

(41] Op. cit., Ref. 35, pp. 30-31.

(42] Qp. cat., Ref. 12, p. 52-53.

(43] Qg. Oj^., Ref. 37, p. 111.

(44] E.B. Becker, G.P. Carey and J.T. Oden, Finite Elements An Introduction Volume I, Prentice-Hall,

Englewood Cliffs, New Jersey, 1981, p. 198.

(45] Qp. cit., Ref. 44, p. 93.

_ii

(46] 2g.rte. , Ref. 44, p. 199. {
149

, I

[47] Q2. 21&., Re.. 44, P. 190.

(48) Ref. 35, p. 68.

(49] H. Engels, Numerical Quadrature Aad Cub ►tuLl,

Academic Press, LorMon, 1980.

[50] D.J. Zimmerle and R.J. Bernhard, "Integration of

Certain Singular Boundary Element Integrals for

Applications in Linear Acoustics*, Herrick Laboratory Report 85-21, Purdue University, 1985.

[51] Q2. rte. , Ref. 37, p. 112-113.

(52] Q2. Ref. 37, p. 163.

[53] Q2. Ref. 37, p. 216.