N .
P(r) = P(r) +E af 2Strn) gs (2.23)
j=l Jas; 7

2.4 Impedance Matrix Formulation

Finding values of the Z matrix is the central part for solving the matrix equation.
Each matrix element Z; represents the contribution of the surface integration of j'* element at
the center of i element. The surface integration is numerically done, and its accuracy is very
critical to the overall accuracy of MOM. The seven-point Gaussian quadrature is used for the
numerical integration; an integration of any ‘smooth’ function within a triangle can be

accurately obtained by sampling seven points within the triangle.

1
{ f(r) ds=A > w; f(r) + Och) (2.24)
Iriangle i=l

where A is the area of the equilateral triangle and wj is the weighting value. The error of the

quadrature is on the sixth order of h, the radius of a circle circumscribing the triangle, so h
must be less than 1 for accuracy. A point P inside a triangle can be uniquely expressed in
terms of the areas of subtriangles defined by P and the vertices of the main triangle (see
Figure 2.4). The subtriangles L1, L2, and L3 are called area coordinates and they add up to
1. The area coordinates and weights for the seven points are listed in Table 2.2 [4]. The
coordinates for the point inside the triangle can be expressed as

rj= PLL1 + P2L2 + P3L3 (2.25)

where P1, P2, P3 are the three vertices coordinates of the triangle.

11



because the scattered fields on the surface of the reflectors have the same phase as the
incident fields on the rigid reflectors. Therefore when the source is positioned at 0.5
wavelength away from the center, the incident and reflected fields are also in-phase, adding
up the field amplitudes. The patterns show that the parabolic reflectors do not neccessarily
achieve high directivity. Large side lobes can be found depending on the positions of the
source. Figure 3.16 compares the diffraction patterns of a parabolic reflector, disk, and a
rectangular baffle with similar dimensions. The parabolic reflector achieves slightly higher

directivity, and it is not enough to justify for purchasing expensive parabolic reflector

microphones.

Focus

-—

A typical parabolic reflector

Focus

A disk-like parabolic reflector

Figure 3.10 The parabolic reflectors

32



y Incident planewave

+ Z-axis

p axis

Figure 3.4 Scattering from one sphere

Scattered pressure field amplitude in the azimuthal plane

0.14
0.12
0.1
0.08
0.06
0.04
0.02
0

ka=2

Field distance from the sphere: 5a

Figure 3.5 The BOR simulation and analytic scattering patterns of the sphere



Pressure, linear scale

Pressure, linear scale

Total pressure field amplitude on the surface of the disk
ka=1

i

0 1 i
0 0.2a 04a 0.6a 0.8a a

Distance from the center of the disk

Figure 3.18 The BOR simulation of the disk 1

Total pressure field amplitude on the surface of the disk
ka=2

0 0.2a 04a 0.6a 0.8a a
Distance from the center of the disk

Figure 3.19 The BOR simulation of the disk 2

37



The homogeneous Green’s function satisfies the Helmholtz equation, i.e., the

frequency domain wave equation:

(" +22) yor) = f(r) (2.3)
Cc

The speed of the sound is denoted as c and is assumed to be 340 m/sec in the air.

Substituting (2.1) in (2.3) yields,

v2 +22| Girrs) = - rr) (2.4)
(v? +a)

Multiplying (2.3) by G(r,rs) and (2.4) by y(r,rs) and rearranging terms in each equation,

the following relations are obtained:

Gorn) Vue) = Gerry [try -& yoo} 5)
yr) V°G(r,rs) = - Yn) (5r-9) + of Gir.) (2.6)

The Green's second identity states that for arbitrary scalars U and V,

[ler V VU) dQ= [eRe ds (2.7)

where Q is a volume and S is the surface enclosing Q. Replacing U with the pressure field,

V with the Green’s function and using relations (2.5) and (2.6),
| [yin deer) Gers) + Gersrs)] + Gere) [Ferg - wer] 40
lo 2 2

--| [ve 2S cere Mas (28)
IS

m on

where y(r') is the field on the surface and the partial derivatives were taken with respect to a
normal vector on the surface. The normal vector n points outward from the volume Q. Note

that the minus sign at the right hand side of (2.7) reflects that -n is taken instead of n (see



accuracy is the key in calculating the coefficients.

Once the coefficient ratios are determined accurately, the actual values of the
coefficients can be obtained in terms of an arbitrary coefficient value. Flammer uses a
normalization scheme [12] such that each spheroidal angle function reduces exactly to the
corresponding associated Legendre function when c becomes zero. This normalization
follows Chu and Stratton’s normalization scheme [12] except that the normalization is carried
out at 7) = 0. Thus using the normalization relations below, the coefficients are completely

determined.

SCIP rt gon = C1? a!

0 2t (ef! r= an ep! for neven (4.28)

= afte t! on CUE (at
x seal Fel) do “sell PH) for n odd (4.29)

Once the expansion coefficients are determined, No, from (4.2) is easily found as below

{ill

~ “e ‘0
Noa2 5) (aon?

20,1 (20 + 1)r! (4.30)

4.3 The Radial Functions
The radial functions satisfy the differential equation (4.9). The eigenvalue in (4.9) is
identical as in (4.8) of the angle function. The radial functions are found in terms of the
spherical Bessel, Neumann, and Hankel functions [8].
° et °°
RP Cic,i€) = LE, | x iF? gon Z(cE) (4.31)

where

Py) = 4/ -
WPoO=fE Jqu@) — forp=1

Z(2) = {= Yrs @) for p=2

(4.32)

47



n? (2n - 1)? (2n - 2)? (2n - 3)? . (n+ 1)? +22 M43 (+4?

7 (2n - 7) (2n - 5)? (2n- 3)3 (Qn- 1) (2n+1) (2n+1) (2n+3)* Qn+ 5)3 (2n +7)? (2n+9)

- (n + 1)4 (0 +2)4 - n4 (n- 1)4
(2n + 1)? (2n + 3)7 (2n +5)? (2n'- 3)? (2n- 1)? (2n + 1)?

- (n- 1)? n2 (n + 1)? (n +2)?
(2n - 3) (2n - 1)4 (2n + 1)? (2n + 3)* (2n +5)

(4.18)

In the numerical computation, truncating the summation up to k = 4 yields sufficiently

accurate eigenvalues for small c?. For a large argument, an asymptotic expansion is used

(11).

on(-ic) = - c2 + 2c (2V+ 1)- 2v(v41)-1+ Aon (4.19)

where v = 5 for n even and v= ah for n odd. Agn is defined as

Aon =X Bx ek
k=1

Bi'=-23 (q+), By" =-26[Sq* + 109? + 1}

BS” = - 2°9 q[33q4 + 114q? + 37] , Bo” = - 2-10 [6396 + 340g + 239q2+ 14] (4.20)

where q = n+ 1 forn even and q =n forn odd.

The power series and asymptotic expansion methods are not sufficient for
intermediate values of c2, and the eigenvalue should be refined by Boukamp’s method of
approximation. The eigenvalue is obtained from a transcendental equation U which is

originated from equations (4.15) and (4.16) [1 1]

U(Aon) = Ur (Aon) + Ua(Aon) = 0 (4.21)

where

45



gron(cic it) = © r+ 1) BYE" for neven
=0
‘ (4.47)

g'on(-ic,ié) = 5 2r BSP er for n odd
r=

Therefore the numerical means of calculating the radial functions are completely determined.

4.4 Analytic Solutions

The exact solution in (4.2) have been calculated in terms of the angle and radial oblate
spheroidal functions. Although (4.2) is a short equation, writing a FORTRAN code for it has
been challenging. There is a limited listing of the eigenvalues, coefficients, angle functions
and radial functions in [11] and [12]. A portion of those values listed in [11] and [12] have
been verified when debugging the code. In spite of the inadequate information, the code
calculates satisfactorily. Figures 4.2 and 4.3 compare the analytic solutions obtained by the
code and the corresponding BOR simulation results for thin disks. The exact solution in
(4.2) involves a summation from zero to infinity for n, but a summation from zero to 30
proves to be adequate for the code. Comparing the analytic solution has proved that the BOR

code is valid for a point source excitation as well.

51



REFERENCES

[1] C.A. Balanis, Advanced Engineering Electromagnetics, Wiley, New York, 1989.

[2] G.W. Swenson, Jr., E.R. Sandeen, L.L. Pater, and H.C. Zhuang, "The potential
for mitigation of gun blast noise through sheltering of the source," Technical Report
N-92/09, United States Army Corps of Engineers, Construction Engineering
Research Laboratory, April, 1992.

[3] N. Morita, N. Kumagai, and J. R. Mautz, Integral Equation Methods For
Electromagnetics, Artech House, Boston, 1990.

[4] J.M. Jin, The Finite Element Method In Electromagnetics, Wiley, New York,
1993.

[5] J. W. Benson, Y. L. Li, and G. W. Swenson, Jr., “A baffle-type directional
microphone,” J. Acoust. Soc. Am., vol. 95, pp. 2536-2538, 1994.

[6] A. F. Seybert, B. Soenarko, F.J. Rizzo, and D.J. Shippy, “A special integral equation
formulation for acoustic radiation and scattering for axisymmetric bodies and
boundary conditions,” J. Acoust. Soc. Am. , vol. 80, pp. 1241-1247, 1986.

(7] J.J. Bowman, T.B.A. Senior, P.L.E. Uslenghi, Electromagnetic and Acoustic
Scattering by Simple Shapes, Hemisphere Publishing Corp., New York, 1987.

[8]S. Wahlstrém, “The Parabolic Reflector as an Acoustic Amplifier,” J. Audio Eng. Soc.,
vol. 33, pp. 418-429, 1985.

[9] A. Leitner, “Diffraction of Sound by a Circular Disk,” J. Acoust. Soc. Am., vol. 21,
pp. 331-334, 1949.

[10] F. M. Wiener, “The Diffraction of Sound by Rigid Disks and Rigid Square Plates,” J.
Acoust. Soc. Am., vol. 21, pp. 334-347, 1949.

[11] C. Flammer, Spheroidal Wave Functions, Stanford University Press, Stanford,
1957.

[12] M. Abramowitz and I. A. Stegun, Handbook Of Mathematical Functions With
Formulas, Graphs, And Mathematical Tables, Dover Publications, Inc., New
York, 1970.

55



4.2 The Angle Functions

In the differential equation (4.9), Son is obtained from the expansions in the

associated Legendre functions of the first kind [1 1):

Son(-ic.n) = De dfn(-ic) P,(m) (4.10)

The prime over the summation sign indicates that the summation is over only even values of
r when n is even, and over only odd values of r when n is odd. Substituting (4.10) in (4.9)
and the use of the associated Legendre differential equation and the recursion formulas for

the associated Legendre functions yield the following recursion formula for the coefficients

9” (11):
(r+2)(rt+ 1)c? on
(r+ 3) (2r+5) “7
Le cet) doa tet et} aon PED 408 = 20) (4.11)

(2r - 1) (2r +3) “| (2r - 2) (2r-1)

Finding the solution to the above recursion formula is the key step in obtaining the angle
functions since the associated Legendre functions are relatively well-known functions. A
FORTRAN subroutine from Numerical Recipes was used for the Legendre functions (see
Appendix C).

For convenience, let

= 1.2 __1l___|
neree)-Leli+ampagy]: 62% (4.12)
12 (r- 1)* c4 >?
- cf lg 4.13
Br (ar - 1)? (2r - 3) (2 + 1) 22) @1)
N;= r(r-te?_ di" a9 (4.14)

(2r - 1) (2r + 1) a3

Substituting above notations in (4.11) yields,

43



| Total pressure field amplitude in the azimuthal plane

1.2
1

0.8
0.6
0.4

0.2

Tilt angle: 30°

Figure 2.8 MOM simulations of the tilted baffle 2

When one fires a gun inside a shed, the acoustic excitation may be represented as a
point source. A simple shed has been modeled as in Figure 2.9, and the diffraction patterns
have been studied. Figure 2.10 shows the diffraction pattern of a point source inside a shed.
Also a simple hill has been modeled to see the diffraction patterns of gun blast noises nearby
a hill (see Figures 2.11 and 2.12). The dimensions of the hill is small compared to a realistic

hill, and analyzing a larger hill requires more CPU time and memory.

15



Total pressure field amplitude in the azimuthal plane

1.2

08
0.6
04
0.2

<a Sounce at the surtacy

Source at 0.25 wayeleng Diameter: 2 wavelengths

Source at 0.5 wavelength
ee velengt Focus: 0.25 wavelength

Source at 1 wavelength

Figure 3.11 The BOR simulations of the parabolic reflector 1

Total pressure field amplitude in the azimuthal plane

1,2
1
0.8
0.6
0.4
0.2
0

Source at the surface
Source at 0.25 wavelength
Source at 0.5 wavelength Focus: 0.5 wavelength
Source at 1 wavelength

Diameter: 2 wavelengths

Figure 3.12 The BOR simulations of the parabolic reflector 2

33



_——————

Form Approved

| REPORT DOCUMENTATION PAGE OMB No. 0704-0188

timated to average 1 hour par response, including the time for reviewing instructions, searching existing data sources,
| gathering and maintaining the data needed, and completing and reviewing the collection of information. ‘Send comments regarding this burden estimate or any other aspect ofthis

Sollection of information, including suggestions for reducing this burden, to Washington Headquarters Services, Directorate for information Operations and Reports, 1215 Jefferson
Davis Highway, Suite 1204, Arington, VA 22202-4302, and to the Office of Management and Budget, Paperwork Reduction Project (0704-0188), Washington, DC 20503.

1. AGENCY USE ONLY (Leave Blank) | 2. REPORT DATE 3, REPORT TYPE AND DATES COVERED
May 1 Final

Public reporting burden for this collection of information is e:

997

4, TITLE AND SUBTITLE 5, FUNDING NUMBERS
BT25

Diffraction of Sound by Objects in the One-wavelength Size Range
TE6

6. AUTHOR(S)
Soon Young Lee

8. PERFORMING ORGANIZATION

7. PERFORMING ORGANIZATION NAME(S) AND ADDRESS(ES)
REPORT NUMBER

U.S. Amy Construction Engineering Research Laboratories (USACERL)
P.O. Box 9005 TTM 97/73
Champaign, IL 61826-9005

10. SPONSORING / MONITORING

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

US Army Center for Health Promotion and Preventive Medicine
ATTN: MCHB-DC-EEN

Aberdeen Proving Ground, MD 21010-5422

11, SUPPLEMENTARY NOTES,
Copies are available from the National Technical Information Service, 5285 Port Royal Road, Springfield, VA 22161.

12a, DISTRIBUTION / AVAILABILITY STATEMENT 12. DISTRIBUTION CODE

Approved for public release; distribution is unlimited.

13. ABSTRACT (Maximum 200 words)

14, SUBJECT TERMS 15, NUMBER OF PAGES
acoustics sound waves -- scattering 64
method of moments 16. PRICE CODE
sound -- measurement

17, SECURITY CLASSIFICATION 18, SECURITY CLASSIFICATION 19. SECURITY CLASSIFICATION 20. LIMITATION OF

OF REPORT OF THIS PAGE OF ABSTRACT ABSTRACT
Unclassified Unclassified Unclassified SAR
NSN 7540-01-280-5500 ‘Standard Form 298 (Rev. 2-89)
Prescribed by ANS! Std 239-18

298-102



Applying the above scheme to the line integral Equation (3.13) yields,
N 3 | ;
Cr) yor) = x x Vine K(r,r'(§)) No(E) PG) Im(G) d& + 4m wer) (3.17)
m=! a= at

where

3 1
Cr) =4n+ yi Ka(r,r'(&)) P(E) Im(E) 46 (3.18)
m=1 J-]

N is the total number of elements used to discretize the scatterer body, and J,,(&) is the

Jacobian of transformation for element m.
dp \2 2 1/2
In @)= (3 +(@P) =
ab) lo

[4 (93 - p1)& +(p1 +p3- 2 p27 +4(23-21)8 +(e1 +23 -2.207 |!

(3.19)

Therefore for any point P on the surface of the scatterer, whose global node number is j, the

integral equation can be expressed as follow:

3°53 N .
-L YL Vine a%j +l +z Cri} = yi" (3.20)
m= a=1 m=1
where
1
a-2[ K(P),8) No(E) p(E) Im(&) dé (3.21)
1
1
and Cinj | K2(Pj,6) P(E) Im(&) a& (3.22)
1

Equations (3.20) to (3.22) apply to each node, and those M equations make up a

matrix equation for BOR.

A vy [=| wire (3.23)

For integration on the axis, the Gaussian quadrature was used. Initially a 4-point Gaussian

quadrature formula was used to evaluate the integral. Comparing with the theoretical solution

26



Figure 2.4 Point P within a triangular element

Table 2.2 The area coordinates and weights

i Ll 12

1 1 1 1 270
3 3 3 1200

2 (9 + 215) (4-VT5) (7 - ¥T5) (155 - 15)
21 21 21 1200

3 (4-VT5) (9 + 2V15) (7 - V5) (155 - ¥I5)
21 21 21 1200

4 (4-115) (4-115) (13 + 2VT5) (155 - 15)
21 21 21 3200

5 (9 - 2V15) (4+ 115) (7 +115) (155 + ¥T5)
21 21 21 1200

6 (44115) (9 - 2V15) (7 +115) (155 +415)
21 21 21 1200

7 (44115) (4+115) (13 - 2¥75) (155 + ¥15)

2.5 Simulation Results
To verify the accuracy of the MOM program, a simulation of scattering from a

wooden baffle was run and compared to the physical measurements performed outdoors.

12



Total pressure field amplitude in the azimuthal plane

FAY
06 Bilin,

Frequency: 100 Hz

Figure 2.12 The MOM simulations of the hill

2.6 Problems In MOM

Although Alona Boag’s program accurately predicts the scattering pattern, there are
several difficulties in using her codes. First of all, the mesh generator she used is not user-
friendly. For example, to mesh the rectangular baffle, one has to enter four vertices
coordinates for each rectangular face of the baffle counter clockwise, looking from the
outside of the baffle. This ensures the correct calculation of outward normal vectors. It is
simple for the baffle with six faces, but when a geometry is complex and has many facets, it
is tedious to enter the coordinates manually. In addition, an accurate meshing of any curved
surface is impossible since the mesh generator assumes a flat surface. To create geometry
and mesh its surfaces with ease, a finite element analysis package called PATRAN is

explored. A few simple PATRAN commands create geometry, mesh surfaces, and output the

18



Total pressure field amplitude in the azimuthal plane

1.2
1

0.8
0.6
0.4
0.2
0

Frequency: 500 Hz

Figure 2.6 The MOM baffle simulation and measurements 2

Side view

+ Z-axis
oom |

Figure 2.7 The tilted baffle

14



revealed that the 4-point Gaussian quadrature introduced 20 - 30% error. 6- or 8-point
Gaussian quadrature must be used for an accurate numerical integration. In addition the
number of total nodes were varied, and the number of nodes did not contribute much to the
accurate analysis. For example, scattered patterns from 9, 17, and 33 node-bodies were
similar to one another. Thus beyond a reasonable number of nodes, the number of total

nodes will lengthen the computational time without any improvement on the accuracy.

3.3 Scattering From Spheres

Several cases were run to verify the validity of the code. To compare with the results
from the reference paper [6], identical test parameters were chosen. Seybert et al presented
the scattered pressure field patterns of spheres due to an incident plane wave [6]. The test
case results closely agree with the results from the paper. In addition, a theoretical scattering
pattern for Figure 3.4 was computed in Mathematica for a further comparison. Due to its
simple geometry, scattering by a sphere can be solved analytically. For a rigid sphere with a

plane wave incident in -z axis, the total field is [7],

yet = pine gy = Si)" (2 + 1) [jnke) - aly hOCkr)] Pa(cos®) (3.24)
n=0

where jn and h are spherical Bessel and Hankel functions, respectively. Also a'y = An(ka)
no'(ka)
where a is the radius of the sphere. P, is the associated Legendre function.Since the incident
field is exp(-ikr cos®), the scattered field is y*' - exp(-ikr cos@). In the actual calculation of
the analytic solution, the summation in (3.24) were added up to n = 30 (see Figure 3.5). The
test cases involving only one scattering body were compared to the analytic solutions, and
both my results and [6] closely followed the analytic solutions (see Figure 3.6). For test
cases involving two discrete spheres, approximate solutions, which neglected the interactions

between the two bodies, were compared to the numerical solutions (see Figures 3.7, 3.8,

and 3.9).

27



rr EE ————

| which solves the IE by dividing the scatterer surface into small patches called elements. The
elements can be any shape as long as the shape can be used tightly to mesh the surface,
leaving no “holes” or space on it. The simulation program written by Alona Boag uses
triangular elements to mesh the surface. For accurate simulation results, the lengths of three
sides of each triangular element must be approximately equal, i.e., isometric. Creating an
element with one extremely long leg compared to the others is not recommended. Also the
length of each triangle leg, denoted by grid length, should be less than one tenth of a wave
length.
Assuming the pressure fields do not very drastically within each element, P(r’) in
(2.15) can be regarded constant on each element. Usually the value of P(r’) on the center of
each element is picked; this is how the method of moment is named. The process of

assembling a matrix equation is described for a typical rectangular barrier (see Figure 2.3).

XKKX
KKK
OOO,

Let N be the total number of the triangular elements on the surface, and define pulse basis

functions

f,=1, on AS};

18
f,=0, onal other AS} “3 (2.16)

where AS’; is the j'® triangular element on the surface. Let r'j be the center coordinate of



CHAPTER 5
CONCLUSIONS

In this thesis numerical methods for calculating diffractions of acoustic waves have
been investigated. Diffraction patterns have been accurately calculated using the MOM code
for simple objects small compared to their wavelength. Chapter 2 introduces the MOM and
presents the simulation results. The MOM code is most useful for studying diffractions of
rectilinear objects, and it is not suitable for analyzing curved objects due to the unfriendly
meshing scheme. The BOR code is written to compensate for the shortcomings of the MOM
code. Diffraction patterns of axisymmetric objects with axisymmetric boundary conditions
can be efficiently calculated by the BOR code. The diffraction patterns of spheres, parabolic
reflectors, and thin disks have been successfully calculated, and the simulation results are
shown in Chapter 3. Diffractions due to a plane wave incidence are easily compared with
existing references. For a point source incidence, no graphical data are available, and only
the analytic solution is expressed in terms of the spheroidal wave functions. In Chapter 4 the
oblate spheroidal functions are discussed in detail. A code is written to calculate the analytic
solution, and the calculation results have verified the validity of the BOR code for a point
source incidence.

Note that the boundary conditions must also be axisymmetric for the BOR
simulations. All of the BOR simulation results in this thesis are done in the full space. For a
plane wave incidence, the direction of the propagation is parallel to the axis of symmetry, and
for a point source incidence, the source is located on the axis. Analyzing objects in a half
space, such as spheres above the flat ground, can also be done using the BOR code. Placing
acoustic images of the original objects in the full space is equivalent to the original objects in
the half space which loses the symmetry.

The BOR code can further be improved if it is modified to account for off-axis

incidence. The current code assumes that there is no field variation on a plane perpendicular

53



CHAPTER 1
INTRODUCTION

Loud noises such as gun blast noises from Army training facilities have brought
concerns from nearby communities. In response to this environmental issue, the Army has
built and implemented noise monitoring systems and noise barriers. For an accurate and
effective solution to noise problems, fundamental acoustic scattering and propagation should
be carefully studied.

Studying the acoustic scattering has been difficult and incomplete by theoretical
means. An efficient and accurate numerical method is necessary to investigate acoustic
scattering. A comprehensive investigation was done by Alona Boag under the direction of
George W. Swenson, Jr. at United States Army Construction Engineering Research
Laboratory (USACERL) Acoustics Team. Alona Boag concluded that method of moments
(MOM) was the best approach in solving acoustic scattering problems. MOM involves in
solving an integrodifferential equation in a matrix equation form. She has written an MOM
code which calculates the radiation patterns of an acoustic source scattered by an object. Her
simulation result closely matches experimental data taken from a scattering of a rectangular
baffle.

Her MOM program is an excellent way of solving scattering problems of small and
simple objects. However her geometry meshing scheme, explained in Chapter 2, becomes
tedious when creating curved or complex geometry. To effectively create 3-D objects, a
package called PATRAN has been explored. A few simple PATRAN commands create
geometry, mesh surfaces, and output node coordinates.

Another difficulty is in analyzing an object large compared to its wave length, 2. The
number of unknowns to be solved grows as an order of (d/A)*, and the central processing
unit (CPU) time grows as an order of (Ay, where d represents a dimension of the object.

For axisymmetric geometry, a method called, body of revolution (BOR) technique,



E-K
kk?

ez)

lS

(3.11)

2m
where E -| (1 - ’sin29)' 7a and k= (1 - RP). Substituting above expressions and
0

evaluating derivatives in (3.10),

k*R

Ko(rr) = s le -RE_.

cae 9 - EGo-2) nal (3.12)

Thus the surface integral in (3.1) is expressed as a line integral as below:

Cor) wr)

{ wr) K(r,r') p(r') dL(r') + 4x w"(r) (3.13)
IL.

where C(r) = 40 + [ K2(r,r’) p(r') dL(r').

3.2 Numerical Implementation
Equation (3.13) can be numerically evaluated along a contour L of an axisymmetric

body. The line L is divided into segments for integration (see Figure 3.2).

Z Element 1
Node 1
2
Element 2
3
4 >
p
5o/
6
7
9] 8
Element 3
Element 4

Figure 3.2 Discretization on the contour L

24



CHAPTER 3
BODY OF REVOLUTION

As discussed in Chapter 2, MOM has many limitations in analyzing large, curved
scatterers. A common example of such a scatterer is a parabolic reflector. A parabolic
reflector is axisymmetric, so a method called body of revolution (BOR) can be used to solve
such a case. The study of the BOR technique is based on a paper by Seybert et al. [6]. By
taking an advantage of the axisymmetric properties, the surface integral is reduced to a line
integral along the generator of the scatterer body and an integral over the angle of the
revolution. The integration over the angle is performed partly analytically in terms of elliptic
integrals and partly numerically using the Gaussian quadrature formula. A program in
FORTRAN is written, and the same test cases in the Seybert’s paper are run and compared to
the paper. My results agree well with both the paper and the theoretical solution. This chapter

describes BOR in detail.

3.1 Derivation Of BOR

The derivation of BOR is based on the Helmholtz integral formula,

C(r) vr) = | ver) sGer) ds + 4m yire(r) GB.1)
where
of 41 .
C(r) =4n+ { ———~| ds(r') (3.2
. aR i} )

Throughout the derivation of BOR, the time dependence of exp(iat) is used, and the Green’s
function is taken as G(r,r’) = exp(-ikR)/R to follow the notations of the Seybert paper. The
integral in (3.2) approaches -27 on the smooth surface of the body.

Figure 3.1 illustrates a simple axisymmetric scatterer and the nomenclature needed in

the derivation of BOR. When the scatterer, boundary conditions, and sources are

21



TABLE OF CONTENTS

CHAPTER PAGE
L INTRODUCTION oicccececeeseccceseesneeeeeenneeeeecnsaeeserraneeseeseeeeessenaeeseneaanes 1
2. METHOD OF MOMENTS (MOM) ......:c::cceeeeeeeeeeeeseessssseeeeeeeeneeeeeeeeeeeseeees 3

2.1 Mathematical Derivation Of The IE ...
2.2 Evaluation Of The Singularity .
2.3 Matrix Equation Assembly
2.4 Impedance Matrix Formulation ..
2.5 Simulation Results
2.6 Problems In MOM

3 BODY OF REVOLUTION (BOR) ...
3.1. Derivation Of BOR .

3.2 Numerical Implementation .. . 24
3.3 Scattering From Spheres .... . 27
3.4 Scattering From Parabolic Reflectors . . 31

36

3.5 Scattering From Thin Disks .

4 SPHEROIDAL WAVE FUNCTIONS
4.1. The Oblate Spheroidal Geometry .
4.2. The Angle Functions
4.3 The Radial Functions
4.4 The Analytic Solutions ...

5 CONCLUSIONS ...

REFERENCES

DISTRIBUTION ooo... eeecceeeeeeeecceneceseeeeestceteieeetreesenseeesseernaeenneetnness 56

(DTIC QUALITY INSPEITED 8



USACERL TM-97/73

Foreword

This study was conducted for USACHPPM under Project 4A161102BT25,
“Environmental Research—Corps of Engineers;” Work Unit IE6, “Terrain, Weather
and Source Direct Modeling for Blast Noise Propagation.” The technical monitor
was William Russell, USACHPPM.

The project was carried out under the U.S. Army Research Office Broad Agency
Announcement DAAH04-94-R-BAA3 dated January 1994, as part of the U.S. Army
Environmental Quality Basic Research Program. The point of contact was Larry L.
Pater, Planning and Mission Impact Division (LL-P) of the Land Management
Laboratory (LL), U.S. Army Construction Engineering Research Laboratories (USA
CERL). The author would like to take this opportunity to thank her advisor,
Professor George W. Swenson, Jr., for providing her with the opportunity to work
on this project. His encouragement and support have been invaluable to her
Master’s study at the University of Illinois. She would also like to thank Ms. Alona
Boag, who lay the foundation for this project, and Dr. Yuan-Liang Li, who taught
her avout the computational methods in acoustics. Dr. Harold E. Balbach is Chief,
CECER-LL-P and Dr. William D. Severinghaus is Operations Chief, CECER-LL.

Dr. Michael J. O’Connor is Director of USACERL.

iv



AS‘; . Also let the unknown pressure fields on the surface be expressed by,
N
Pr’) = > pf (2.17)
i=l
where p is a 1 x N unknown array for the surface pressure fields. Substituting (2.17) into

(2.15), and matching the field at the midpoint of each A’s, a matrix equation is obtained as

below:
pitt = z Zij Pj » 1=1,2,3,...N (2.18)

Expressing (2.18) in a matrix form,

p |= Z pinc (2.19)
where
dG(ri,r') - oe
Zij= —— ds 1F4j
4 I, on
(2.20)
Zy=5 »i=j

The derivative of the Green’s function with respect to a normal vector is

a ik exp(ik [rj - al ".
= (ik “im Ma) adr ri-r')-n (2.21)

The N x N matrix Z is called an impedance matrix, and it needs be inverted to find the

unknown surface pressures in p.

Zz

P= % Zi if pine (2.22)

i=l

To find a total pressure at any location r, sum the incident field and the scattered field found

by the above method;

10



3.18 The BOR simulation of the disk 1
3.19 The BOR simulation of the disk 2
3.20 The BOR simulation of the disk 3
3.21 The BOR simulation of the disk 4 ...
3.22 The BOR simulation of the disk 5 ...

4.1 The oblate spheroidal coordinate system .........s.sssseessseeeseesereeeseeeeeenees 41
4.2 The BOR simulation and analytic diffraction solution of the disk 1 ...
4.3. The BOR simulation and analytic diffraction solution of the disk 2 ...



Top view Side view

15m

7.0m

Front view

+ Z-axis A

+ X-axis

Figure 2.9 The shed

16



US Army Corps

of Engineers

Construction Engineering USACERL Technical Manuscript 97/73
Research Laboratories ly 1997

Diffraction of Sound by Objects in the
One-wavelength Size Range.

by
Soon Young Lee

19970528 026

The diffraction of sound by hard geometrical shapes the computational load. The validity of the resulting
has been studied by the Method of Moments (MOM), a procedure has been verified by comparing the MoM
boundary element technique for solving the scalar results for a flat circular disk illuminated by either a
Helmholtz integral equation. The scattered fields of plane or spherical wave with the results obtained by the
several rectilinear-shaped objects are analyzed. The classical method of separation of variables and
problem of scattering by an axially-symmetric obstacle expansion in oblate spheroidal wave functions.
has been simplified by using the symmetry to reduce

DTIC QUALITY KiSPEGTED 3

Approved for public release; distribution is unlimited.



Total pressure amplitude on the surface of the disk
ka=5

Pressure, linear scale

0 0.2a 04a 0.6a 0.8a a

Distance from the center of the disk

Figure 3.22 The BOR simulation of the disk’ 5

39



CHAPTER 2
METHOD OF MOMENTS

There are several methods for solving an acoustic scattering problem. A well known
method is geometrical optics (GO). GO yields an approximate solution for scattering from
bodies large compared to the acoustic wavelength [1]. Because it is an approximate solution,
only applicable to the analysis of high frequency waves, it is inappropriate for USACERL’s
scattering analysis. For example, a typical spectrum of an Army gun blast noise ranges from
5 Hz to 60 Hz [2]. In addition, accurate prediction of scattering from a complicated geometry
is almost impossible. The above reasons eliminate the further study of GO.

Another well known scattering mode is Rayleigh scattering. Rayleigh scattering is
typically used to analyze scattering from small particles such as air bubbles in liquids and
water molecules in the air. The dimensions of practical noise barriers are comparable to or
larger than the acoustic wavelength, so Rayleigh scattering is not an appropriate model.

Alona Boag has investigated two rigorous numerical methods in solving acoustic
scattering: a partial differential equation method (PDE) and an integral equation method (IE).
She presented her findings at a USACERL seminar, and Table 2.1 summarizes her findings.
Table 2.1 clearly indicates that the IE method is the best choice among available scattering
solutions. The IE method leads to MOM, and the remainder of this chapter is devoted to

MOM.

2.1 Mathematical Derivation of The TE
When an acoustic source radiates a wave in free space, the radiated field w at an

observation point r can be mathematically described as :
Wr) = - { f(rs) G(s) dv (2.1)

where f is the source at rs. The Green’s function G(r,rs) is an impulse response at r due to



eliminates such problems. BOR reduces the number of unknowns to an order of (d/A),
thereby reducing the CPU time. A BOR code written in FORTRAN is used to analyze
scattering from curved, axisymmetric bodies. The BOR code has been extensively tested for
its validity, including a comparison to an exact solution for scattering from oblate spheroidal
objects involving expansion of the field in series of spheroidal functions.

Chapters 2 and 3 of this thesis present the study of MOM and BOR in detail:
mathematical derivations, code implementation, and simulation results. In addition, an in
depth study of oblate spheroidal functions is given in chapter 4. Chapter 5 presents the

overall results and conclusion.



Total pressure field amplitude in the azimuthal plane

Diameter: 2 wavelengths

Point source on the surface

1.2
1

0.8
0.6
0.4
0.2

Figure 4.2 The BOR simulation and analytic diffraction solution of the disk 1

Total pressure field amplitude in the azimuthal plane

Point source on the surface

Figure 4.3 The BOR simulation and analytic diffraction solution of the disk 2

52



Total pressure field amplitude in the azimuthal plane

Diameter: 2 wavelengths

Focus: 100 wavelengths
“flat disk”

Source at 0.
Source at 0.5 wavelength
ath |

Total pressure field amplitude in the azimuthal plane

1.2
1

0.8
0.6
0.4
0.2
0

Plate
Disk
Parabolic reflector

Disk
|__--" Diameter: 1 wavelength

Parabolic Reflector Point source centered on the surface

Diameter: 1 wavelength
Focal length: 0.22 wavelength
Point source at focus

‘Square plate
Side: 1 wavelength
Point source centered on the surface

Figure 3.16 The BOR simulations of the parabolic reflector, disk, and baffle
35



The above equation can be decomposed into two parts,

an
_ exp(-ikR(r,r’) -1 :
Ky(rr') = { o{are Rier) ePCKRE) 1) goer) (3.5)
phi
and K.(nr) = f Stacey) doc’) (3.6)

As R approaches zero, the integrand of K;(r,r’) can be expressed in terms of a Taylor’s

series expansion,

exp(-ikR(rvr’) - 1 _ 1 - ikR(rr') + 0.5 (-ikR(rr'))? - 1 =-ik- k? R(r,r') BD
R(rr') R(r,r') 2 .
Let the normal component of R(r,r') be d, and the tangential component be x. Then the

derivative of (3.7) becomes,

es Rie)) 2 Rie) _ Kd (38)
on 2 2 on 24a

As x approaches zero, (3.8) will be -K, so Kj(r,r') is non-singular and can be evaluated

numerically at R(r,r') = 0. K2(r,r') can be expressed in terms of the elliptic integrals.

2n
9 on 4 dg _ 014
wen =h Re —1 _ aor 1-3 ne] ooo = jn [SA &] (3.9)

lo (1 -R sine)

where F(n/2, k) is a complete elliptic integral of the first kind with modulus Kk. Also,

R? =(pq + pr? +(zg+zpP and k= pape. The outward normal vector at point Q is
R

defined as f= npgp+n,q2, and the angle of revolution is 6 = @q - @p. Evaluating (3.9)
yields,

oF ok . OR Rez ok OR (3.10)

Karr’) = 4 [RS - Fo -FOR |,
arr) = ae ak 9Pq OPQ) Ik 0zq 3zq

The derivative of F is expressed in terms of the complete elliptic integral of the second kind,

23



CHAPTER 4
SPHEROIDAL WAVE FUNCTIONS

Chapter 3 describes BOR in detail and presents the BOR simulation results compared
with the experimental, analytic and numerical solutions from [5, 6, 7, 8, 9, and 10].
However those results are valid only for an incident plane wave, and no solution is available
for a point source excitation. In order to verify the BOR scattered patterns for a point source
excitation, an analytic solution for an axisymmetric geometry must be sought. Studying the
oblate spheroidal geometry is a good way to verify the BOR codes. A “fat” spheroid

represents a sphere, and a thin spheroid becomes a disk.

4.1 The Oblate Spheroidal Geometry
The oblate spheroidal coordinates (§,n,) shown in Figure 4.1 are related to the

rectangular Cartesian coordinates (x,y,z) by the transformation

xelav[?et (1-1?) cos
y=Lav(e?+i)(1-1) sino 4)

al
z=7dén

where 0 < § <e, -1 <1 1, and 0 << 2m. The z-axis is the axis of symmetry, and the
interfocal distance, minor axis, and major axis are d, dé and dV e +1 , respectively [6].
Because of the axisymmetry there is no 6 variation, so ¢ is assumed to be zero in the rest of
this chapter.

The exact solution for an acoustically hard oblate spheroid due to a point source

excitation on Z-axis (io = 0) is computed by Bowman et al. [7] as

40


axisymmetric, the pressure fields and the normal velocity fields on the surface of the scatterer
are also axisymmetric. Taking an advantage of this symmetric distribution of the fields, the

surface integral may be reduced to a line integral.

Figure 3.1 A simple axisymmetric scatterer

Because the axisymmetric body is independent of the angle of revolution 6, y(p,0,z)

depends only on z and p. The integral in (3.1) may be expressed as below:

2m .
. 0 (exp(-ikR(r,r') . , .
| volf nt Raa ja p(n’) dL(r')
L

dyer’) { (" expCiKROP) sore | oven at ce
| af Rar) d@(r')| p(r') dL(r') (3.3)
L
For convenience, let
20 .
yO 0 fexp(-ikR(r,r') '
K¢r,r') -{ S{nreee ee) d@(r') (3.4)

22



Figure

2.1
2.2
2.3
2.4
2.5
2.6
2.7
2.8
2.9
2.10
2.11
2.12

3.1
3.2
3.3
3.4
3.5
3.6
3.7
3.8
3.9
3.10
3.11
3.12
3.13
3.14
3.15
3.16
3.17

LIST OF FIGURES

A scatterer and a source in an unbounded space ..

Limiting procedure...........ccccceeeeeeeee
The meshed rectangular barrier and the coordinate system ....
Point P within a triangular element
The MOM baffle simulation and measurements
The MOM baffle simulation and measurements 2.. e
The tilted baffle... cece ce ccs cesceneceserseeseeeeeeceeeesseeesestenseeenae
MOM simulations of the tilted baffle ..
The shed

The MOM simulations of the hill...

A simple axisymmetric scatterer ...............0c00000cseseeeeeseeseseseesseseeeeeeereee 22
Discretization on the contour L

The element mapping onto the local coordinate system ..

Scattering from one sphere wo. ececeessesecssesseseeseesseseessesscsessecsscseee 28
The BOR simulation and analytic scattering patterns of the sphere ..
The BOR simulation of the sphere
Scattering from two spheres

The BOR simulations of the two spheres 1 ..
The BOR simulations of the two spheres 2 ..

The parabolic reflectors
The BOR simulations of the parabolic reflector 1 .... we 33
The BOR simulations of the parabolic reflector 2 ....

The BOR simulations of the parabolic reflector 3
The BOR simulations of the parabolic reflector 4 ....

The BOR simulations of the parabolic reflector 5 .
The BOR simulations of the parabolic reflector, disk, and baffle .. .
The thin disk... cece eesseesceeceeceseccsseessecessecssesessserseerstecseeseas



The baffle was 61 cm wide and 30.5 cm high, and a microphone was placed at the horizontal
center, about 7.6 cm above the ground, and 2.5 cm away from the face of the baffle. The
experimental data are taken from Benson et al. [5]. In the numerical analysis a point source is
| placed at the microphone position using the reciprocity theorem. The simulation results
| closely match the measurements (see Figures 2.5 and 2.6). Also diffraction patterns for the

same baffle but mounted at a tilt angle have been simulated (see Figures 2.7 and 2.8).

Total pressure field amplitude in the azimuthal plane

1.2
1
0.8
0.6
0.4
0.2
0 ar S
: |

ry

Frequency: 400 Hz

Figure 2.5 The MOM baffle simulation and measurements 1

13



(2r + 2) (2r + 3) B9",> + [(2r + 1) (2r + 2) - Aog(-ic)] BS® + 2 BY",

. (4.42)
=-2 Qb,(-ic) [cic DX cSR-ic) 2k (k+ t 1) for n even,
ksr+l
(2r + 1)(2r + 2) B9%,> + [2r (2r + 1) - Aga(-ic)] BS" + c2 BO, (4.43)

=~ 2 Qha(-ic) [ir OCI) E ecto) (2k + 1)(K)- Z BeCic) 2k(K=1)} for n odd.

ker+1

The initial coefficients BQ" are

The symbol (k) denotes the binomial coefficient « Ke TH

found in [11] as

O° = [c RP Cic,i0)]" - Qjy(-ic) RWCic,i0) for n even

; \ (4.44)
Be =[c RY (-ic,i0)}" for n odd
where
RY Ci -ic,i0) = = i gBMGic) for n even, (-ic,i0) = 0 for n odd
7 af"(-ic)
r=0 (4.45)
2c it! a?

ny:
2c dy'cic) for n odd, RS Gi ic,i0) = 0 for n even
3D dPrCic)

r=0

RW (-ic,i0) =

The remaining coefficients are completely determined from the recursion formulas in (4.42)
and (4.43).
The derivative of R®) is determined by simply taking a derivative of (4.37).

; Hy Dans
Ry CieaiB) = QhyCic) ROP ies) [tanh€) - 4] + Sou) Ronen) 3 CI) go(ieib) (4.46)
+1

where g'on is

50



SSANISNE WIdI4IO

128 ‘ON LINHAd $006-92819 SIONITH ‘NOIVdNVHO
WW NOWdWVHD $006 XO8 Od
aivd SH3SNIDNA 4O SddOO
aowised ae SSIMOLVHORV1 HOUVSSSY ONIUSANIONS NOILONYLSNOD
AWHY 3H1 4O INSWLuVdad



Scattered pressure field amplitude in the azimuthal plane

0.35
0.3
0.25
0.2
0.15
0.1
0.05

d=5RI

Field distance from sphere 1: 3 RI

Figure 3.8 The BOR simulations of the two spheres 1

Scattered pressure field amplitude in the azimuthal plane

0.25

d=5R1

Field distance from sphere 1: 3 R1

Figure 3.9 The BOR simulations of the two spheres 2

30



rT

3.4 Scattering From Parabolic Reflectors

Comparing the simulation results in Section 3.3 is not adequate for a complete
validation of the BOR code. Scattering from the spheres is axisymmetric regardless of the
propagation direction of the incident field. Also only the diffractions due to a plane wave
incidence have been compared and validated. For a complete testing of the BOR code,
scattering due to a point source incident on an object such as a parabolic reflector should be
tested.

Simulations on the diffractions of the parabolic reflectors can be useful since the
parabolic reflectors are used as acoustic amplifiers [8]. Recording bird sounds in the
wilderness is not a trivial task for ornithologists since the bird sounds may not be loud
enough to be distinguished in the presence of background noises. Many ornithologists use
parabolic microphones to record bird sounds. A parabolic microphone is a parabolic reflector
dish with a microphone attached to the dish. Typically the microphone is located at the focal
point of the parabolic reflector because the fields will arrive in equi-phase at the surface of the
reflector. Using the reciprocity theorem, a diffraction pattern due to a point source at a focal
point is identical to the reception pattern of the microphone located at the focal point.

Although the use of the parabolic microphones are common in bird recordings, exact
reception patterns are not supplied by their manufacturers. The receptivity of parabolic
reflectors is thought to be highly directional. Manufacturing parabolic microphones is not as
simple as building flat baffle-type directional microphones because of the curvature. Thus the
price of the parabolic microphones are quite expensive. The diffraction patterns from various
parabolic reflectors can be accurately calculated by the BOR code. Note that as the focal point
is place farther away from the reflector, the curvature of the reflector becomes more flat (see
Figure 3.10). Figures 3.11 through 3.15 show far field diffraction patterns for parabolic
reflectors whose diameters are twice the acoustic wavelength of the interest. The reflector
analyzed in Figure 3.15 is virtually a disk since the focal point is 100 wavelengths away
from the center of the reflector. Except for Figure 3.11, the highest directivity is achieved

when the point source is placed at 0.5 wavelength away from the center. This is plausible

31



Bn and Uj(Aon) = - Buea

Uj (Aon) = - + Aon - ————— —_
Yn-2- Aon - Nn-2 ‘Yn+2 - Aon - Nova

Note that Ny.2 is evaluated by the equation (4.16) and Ny+4 is evaluated by (4.15). Let ae

be the approximate eigenvalue calculated from either of the above two methods, and 5A, be

the difference between the actual eigenvalue and the approximate value.
0 = U(on) = U(Abn + Shon) = UAGR) + 3U (4.22)

by finding the first variation of U due to the variation 5Aon. If the variation 5Aoq is made in

the eigenvalue, the variation on (4.16) is

r,

8Ny42 = Son + Be nN, (4.23)
(N:P

Likewise, the variation on (4.15) is

(NP

Ir

BN, = "7 [8N42 - Bon] (4.24)

By iteration of (4.23) and (4.24) the variations of U are obtained

SU, = Hoy Bo, Bn Bao, Bn Bn-2 Bas | (425)
(Na? (Nn? (Nn.2? (Na? (Nn.2? (Nn-4?

8U2 = Ron] Ssh (Na? (Nns2? | (No)? (Nova)? (Nanos)? / (4.26)
Bn Bn Bn+2 Bn Bn+2 Basa

Substituting (4.25) and (4.26) in (4.22), 5Aon is found to be

= ui (ase)- usage
+ fn, BnBn-2 vo]o[ fast 4 (Nasal (Nea)? |

Bns2 Bn+2 Basa

(Ni? (NaF(Nn-2?

Bon = (4.27)

The new eigenvalues obtained by the above method yield exactly same eigenvalues listed in

[12]. Thus the Boukamp’s method leads to remarkably accurate eigenvalues, and this

46



coordinates. Creating and meshing curvatures are simple with PATRAN. However, there is
a drawback in PATRAN. To define outward normal vectors on surfaces of an object, the
coordinates should be listed in counter clockwise sequence. PATRAN meshes and lists
element coordinates arbitrarily, either clockwise or counter clockwise. Some normal vectors
point outward, and others inward. To remove this arbitrariness, a simple routine is written.
Assuming each element is approximately equal in size, the center coordinates of each element
are summed up and then divided by the number of elements, N. This yields an approximate
center coordinate of the scatterer, and a dot product of a vector from this center of the object
to acenter of any element and the corresponding normal vector should be a positive quantity.
If negative, then it indicates the normal vector is pointing toward the center of the object, and
the direction of the vector should be reversed.

Other difficulties are the memory and CPU time needed in computing. Since the
number of elements N is proportional to the surface area, analysis of an object large
compared to the wavelength is difficult, if not impossible. As N increases the required
memory for storing the impedance matrix increases as N2. In FORTRAN, reserving a large,
contiguous section of working memory becomes difficult when solving for a large number of
unknowns. For example, in the Michigan baffle simulation at 163 Hz, 7,308 triangle
elements are meshed, and the number of the matrix elements is 53,406,864. Each matrix
element occupies 16 bytes of memory for a complex double precision data type. Thus more
than 854 MBytes are needed in creating the matrix. The CPU time for a matrix inversion is
proportional to N>. Most of the acoustic scattering simulations have been solved in a Convex
C420 machine at the University of Illinois. The Convex C420 has four processors with 512
MBytes of memory. It is used for high speed, heavily vectorized computations. However
C420 has proved to be inadequate for the baffle simulation at 163 Hz. An access to a more
powerful machine is obtained from National Center for Supercomputing Applications
(NCSA): Power Challenge account. The Power Challenge Machine is built by Silicon
Graphics, Inc., and it has 16 shared memory multiprocessors. Its total memory is four giga

bytes, and it can perform 300 mega floating point operations per second (flops) per

19



Based on the above relations, Flammer has derived power series expansions for the oblate

radial functions as below [11].

ROGic, 18) = QaCcic) RorCic.i&) {tan 1(e)- 4 + gon(-ic,ié)

with

[x in ic) 0 jg) (28!

1!(2*(4)!]

QoaCic) = $"(-ic)

5 for n even

(1)
Qoa(-ic) mt ate sl o8"(-ic) (1 = 2r}! for n odd
r! [2 (1)?

. The coefficients c3f are

where a8 "(-ic) =
ia
p= ¥ (4 (c + }) ate for n even

o§p=—t_¥ (- th (+3 5) are for n odd

where (rk =r (r+ 1) (r+ 2)---(r+k- 1), (tp = 1.

The function gon in (4.37) satisfies the inhomogeneous radial equation [11]

(2 4 1) G- dog(cic) + ©2"] gon(-ic.i&) = - 2 Qa(-ic) -L REP -ic.iB)
d& dg dé

It is convenient to expand gon in the form

2re1
for n even

gon(-ic,i) = Borg
r=0

gon(-ieif) = © Bo £7" for nodd

(4.37)

(4.38)

(4.39)

(4.40)

(4.41)

since the power series expansions are used for the values of & which are close to zero [11].

After substituting (4.41) in (4.40) and rearranging, the following recursion formulas for BS?

are obtained [11].

49



i, General Comments. (Indicate what you think should be changed to make this report and future Teports
of this type more responsive to your needs, more usable, improve readability, etc.)

5. If you would like to be contacted by the personnel who prepared this report to raise specific questions or
discuss the topic, please fill in the following information.

Name:

Telephone Number:

Organization Address:

6. Please mail the completed form to:

Department of the Army

CONSTRUCTION ENGINEERING RESEARCH LABORATORIES
ATTN: CECER-TR-I

P.O. Box 9005

Champaign, IL 61826-9005



Pressure, linear scale

Pressure, linear scale

Total pressure field amplitude on the surface of the disk

0 0.2a 0.4a 0.6a 0.8a a

Distance from the center of the disk

Figure 3.20 The BOR simulation of the disk 3

Total pressure field amplitude on the surface of the disk
ka=4

0.2a 0.4a 0.6a 0.8a a
Distance from the center of the disk

Figure 3.21 The BOR simulation of the disk 4

38



3.5 Scattering From Thin Disks

Although the parabolic reflector simulations due to a point source could not be
compared to known results, the simulations can still be checked for the validity. When the
focal length of the parabolic reflector becomes large compared to the wavelength, the shape
of the reflector resembles a thin disk. Diffraction of a thin disk has been rigorously analyzed
previously by many scholars [7, 9, and 10]. Leitner presented the diffraction of sound by a
circular disk due to a plane wave incidence [7], and the exact theoretical values based on the
wave functions of the oblate spheroid were plotted in his paper. Also Wiener’s
measurements supported Leitner’s work [10].

In order to compare with the Leitner’s results, the same test cases were run using
BOR code (see Figure 3.17). Figures 3.18, 3.19, 3.20, 3.21, and 3.22 show the BOR
simulation results. The BOR simulation results agree with the Leitner’s exact theoretical
calculations, and they match closer to the experimental data than the Leitner’s results [9].

There is no study on the diffraction patterns due to a point source incidence although
an exact solution exists. The solution is expanded in terms of oblate spheroidal wave
functions which are not simple to calculate. Since the oblate spheroidal wave functions are
not available in any mathematical software package, a computer code is written to calculate

the functions which are discussed in Chapter 4.

Incident plane wave

Shadow a
side <<
—
Bright —+_———
side

<=

Figure 3.17 The thin disk

36



Figure 2.1). The surface includes the scatterer surface and the fictitious surface at infinity

(see Figure 2.1).

Figure 2.1 A scatterer and a source in an unbounded space

After substituting (2.5) and (2.6) on the left hand side of (2.8) and simplifying, (2.8)

becomes,

[we Or) Gary MM] as 29)
on on

Wr) - G(r,rs) f(rs) dQ= |
IQ

S

By equation (2.1), the integral on the left hand side of (2.9) is essentially equal to the
incident field Pi",

All fields in any unbounded medium must satisfy the Sommerfeld radiation condition:

[Z-#}v-oan-- (2.10)
oR

This condition ensures that the field y becomes negligible at a large distance from the


J and Y are Bessel and Neumann functions, respectively. Similar to Hankel functions, the

radial functions are related likewise [12]

ROG -ic,i€) = ROG -ic,i€) +i Oj ic,i§)
(4.33)

ROP Cie iB) = RPC ic,i&) - i RP(-ic,i8)
(4.33) makes defining the radial functions easier. The derivative of ROCic, i€) is easily
obtained by taking the first derivatives of RVC -ic,i§) and RO: -ic,i8).
Co -1 °°
R® (-ic,i&) | yr wo Li it ddr Zc) (4.34)
1=0,1 1=0,1
However calculating the radial function of the second kind for small values of & is not

suitable since the Neumann function converges poorly for small arguments. For small values

of €, a better numerical method should be used to calculate the radial function of the second

kind.
The angle functions and the radial functions are related by joining factors cP [12],
SP-ic,i&) = KP(-ic) RP-ici&) for p= 1, 2 (4.35)
where
n! » don
rane ic) = —=2__ forneven
On an of @ ia

© 3(n+ 1)! D' aon
Ko, (cic) = - ———_—_1= 1 for n odd
O° asl DT]

(4.36)
ic 29(BP! ap” =
(2) 73 :
Kon (ic) = aI x ae for n even
29n(n-1)\}(n+1)f}qon .,
ost

x? * 4On
 (-ic) = 3a 1)! »» d®" fornodd

48



to the symmetric axis. For example, when the plane wave is incident at an angle from the
axis of the symmetry, the incident fields on the surface of the object will vary accordingly.
The scattered fields will have the same variation due to the variation of the incident fields.

Finding out of a way to incorporate this variation in the BOR code will broaden its usage in

analyzing the diffractions.

534



At r=r', P(r’) is approximately equal to P(r) and may be assumed invariant within 5S.
Figure 2.2 indicates that -z is the normal vector, and taking a derivative with respect to z

yields,

2p’ dp'= -FPO)

laa 2 “al ep

Po
ro co Recewerd + (z-2')?

Figure 2.2 Limiting procedure |

As a and z approach zero, the integral at the singularity is reduced to - FPO). Therefore |

when the observation point is taken on the surface of the scatterer, (2.12) should be rewritten

as:

dG(r,r’)

OM as
on s

Pr(r) =4 P(r) +f P(r’) (2.15)

Ss

where f indicates the principal value of the integral. The above equation is solved to find the
pressure fields on the surface by assembling (2.15) into a matrix equation. Once the pressure
fields on the surface are known, then the real IE (2.12) can be solved for any observation

location.

2.3 Matrix Equation Assembly
Assembling a matrix equation from IE is the key step in MOM. MOM is a technique



USER EVALUATION OF REPORT

REFERENCE: USACERL Technical Manuscript 97/73, Diffraction of Sound by Objects in the One-
wavelength Size Range

Please take a few minutes to answer the questions below, tear out this sheet, and return it to USACERL. As
user of this report, your customer comments will provide USACERL with information essential for improving
future reports.

1. Does this report satisfy a need? (Comment on purpose, related project, or other area of interest for which
report will be used.)

2. How, specifically, is the report being used? (Information source, design data or procedure, management
procedure, source of ideas, etc.)

3. Has the information in this report led to any quantitative savings as far as manhours/contract dollars
saved, operating costs avoided, efficiencies achieved, etc.? If so, please elaborate.

4. What is your evaluation of this report in the following areas?

a. Presentation:

b. Completeness:

c. Easy to Understand:

d. Easy to Implement:

e. Adequate Reference Material:

f. Relates to Area of Interest:

g. Did the report meet your expectations?

h. Does the report raise unanswered questions?



Scattered pressure field amplitude in the azimuthal plane

0.35
0.3
0.25
0.2
0.15
0.1
0.05
0

ka=4

Field distance from the sphere: 5a

Figure 3.6 The BOR simulation of the sphere

Incident
planewave

4 + Z-axis

Sphere 1
- -P
p axis
|
d
Sphere 2

Figure 3.7 Scattering from two spheres

29



Total pressure field amplitude in the azimuthal plane

1.2
1

0.8
0.6
0.4
0.2

| B
i Frequency: 20 Hz

Source Location: (0, 10, 0.5) m

Figure 2.10 The MOM simulation of the shed

Top view
| + Y axis
|
i
9.0m 4.0m nd _—
6.0m
Front view
35m —| 4 + Z-axis
|
|
|
5.0m
|
|
l a

—>| + X-axis

T
o
}
3

Figure 2.11 The hill

17



The contents of this report are not to be used for advertising, publication,
or promotional purposes. Citation of trade names does not constitute an
official endorsement or approval of the use of such commercial products.
The findings of this report are not to be construed as an official
Department of the Army position, unless so designated by other authorized
documents.

DESTROY THIS REPORT WHEN IT IS NO LONGER NEEDED

DO NOT RETURN IT TO THE ORIGINATOR



_ RY ici . a .
i RM cied) - Ria RO ici RO (-ic,iE) Son(-ie.n) (4.2)
Non RO (-ic,i&1)

The above notations are consistent with the Flammer’s notations [11]. The coordinate of the
oblate spheroid surface, the source point, and the field point are denoted as (&1,72,0),
(&,0,0), and (&,7),0) respectively. Also &. and &, take the minimum and maximum values

between the source point and the field point coordinates respectively.

An=1 1 = 0.866

Figure 4.1 The oblate spheroidal coordinate system

The oblate spheroidal coordinate system is one of the coordinate systems in which the

scalar wave equation
(v2 +42) y=0 (43)

is separable [11]. To express this equation in spheroidal coordinates, the following relation

for metrical coefficients, he, hy, and hy are defined by

41



The line L is divided into four elements where each element is made up of three nodes. The

contour along the line L is assumed to be varying quadratically, and the variation is called a

shape function, Nq(€). A quadratic shape function is defined as below, and it provides a
more accurate description of the contour than does a linear shape function. & is called a local
coordinate and ranges from -1 to 1, and ois called a local node number and ranges from 1 to
3.

NaG=-t6+48  a=1

Nal) = 1-8 ,a=2 (3.14)

NoG=4+48  a=3

The actual coordinates can be represented in terms of the shape function and local coordinate,

3
p(&) = & Na(6) Pa
ol @.15)

3
26) = x Na(6) Za

The above equations map each element on the body onto a straight line in the local coordinate

g.

p2
pl g3 pl p2 93
——————- — © e >
l 0 1k

Element i Local coordinate

Figure 3.3 The element mapping onto the local coordinate system
Similarly the pressure fields y along L are approximated using the same shape function.

3
Wm() = X Nol) Wine
o=l (3.16)

where Wma is the value at node o on the element m. Because the fields distribution and the
contour are described using the same shape function, the formulation is said to be

isoparametric.

25



Total pressure field amplitude in the azimuthal plane

5 in Diameter: 2 wavelengths
ource at 0.5 wavelength
Focus: 1 wavelength

ree at 1 wavelength

Source at 2 wavelengths

Figure 3.13 The BOR simulations of the parabolic reflector 3

Total pressure field amplitude in the azimuthal plane

1.2
1

0.8
0.6
0.4
0.2

|
Sourve at the surface Diameter: 2 wavelengths
Source at 0.23 wavelength
Source at 0.5 wavelength
Source at 1 wavelength

Focus: 2 wavelegths

Figure 3.14 The BOR simulations of the parabolic reflector 4

34



processor, total of 4.8 giga flops.

In spite of the increased capacity, a more powerful machine is needed in analyzing
larger objects at high frequencies. The limited resources are unavoidable. One should look
for a clever method of solving the problem instead of looking for a faster machine. A
technique that overcomes all of the above problems, meshing, memory, and CPU time, is

sought and implemented in Chapter 3.

20



f(r). In a homogeneous medium,

G(r.) = 0k) (2.2)

where R =|r - rs]. The wave number k is equal to 2/X. The integration is taken over the

entire volume v of the source. Unless otherwise noted, time dependence exp(-iat) is

assumed.
Table 2.1 Alona Boag’s comparison between PDE and IE
Method Partial Differential Equation Integral Equation
Advanta; * applicable to bodies comparable to || ¢ IE is based on the Sommerfeld
Bes or smaller than the wavelength Radiation Condition
p : * the radiation condition
can analyze arbitrary media and . ;
geometry, including automatically satisfied
inhomogeneous bodies * lower number of unknowns
* sparse matrices to invert compared to PDE method
* high accuracy
Disadvantages ||° number of unknown N grows as _ |[* number of unknown N grows as
(d/Ay (d/n)
* requires artificial boundary + applicable only to piecewise
conditions to truncate the mesh homogeneous problems
« medium accuracy « CPU time grows as (any?

In the presence of an object in the medium, a portion of the radiated wave encounters
the object and re-radiates into different directions. The object acts as a scatterer of the original
wave. When a field is measured at r, it is a sum of the fields due to the original wave and the
re-radiated one. The field due to the original wave is called an incident field, and the other
field due to the re-radiated one is called a scattered field. The sum of the incident field and the
scattered is the total field. To distinguish each field, the incident, scattered, and total fields at
rare denoted as Pi"®, PS, and P respectively. Thus the field in (2.1) should be re-expressed
as the incident field Pi"°. The incident field is recalculated using the Green’s function
assuming that the presence of the scatterer does not influence the propagation of the original

wave.



USACERL DISTRIBUTION

Chief of Engineers
ATTN: CEHEC-IM-LH (2)
ATTN: CEHEC-IM-LP (2)
ATTN: CEMP
ATTN: CEMP-CE
ATTN: CEMP-EA (2)
ATTN: CEMP-ZM
ATTN: CERD-L
ATTN: CERD-M (2)
ATTN: CERD-ZA

ACS(IM) 22060
ATTN: DAIM-FDP

CECPW 22310-3862
ATTN: CECPW-E
ATTN: CECPW-FT
ATTN: CECPW-ZC

Norwegian Defense
ATTN: Chief Test & Dvipmt Section

HQ USAREUR & 7th Army
ATTN: AEAEN-EH
ATTN: Unit 29351

US Military Academy
ATTN: MAEN-A
ATTN: Civil Div Director
ATTN: Dept of Geo & Env Engr

US Army Combat Sys Test Activity
ATTN: STECS-SO-R 21010-5059

Commander FORSCOM
ATTN: FCEN-RDF 30330-6000

CEWES 39180
ATTN: Library

CECRL 03755
ATTN: Library

56

US Army ARDEC 07806
ATTN: SMCAR-ISE

Engr Societies Library 10017
ATTN: Acquisitions

US Army Environmental Center
ATTN: SFIM-AEC-NR 21010
ATTN: SFIM-AEC-CR 64152
ATTN: SFIM-AEC-SR 30335-6801
ATTN: AFIM-AEC-WR 80022-2108

National Guard Bureau 20310
ATTN: NGB-ARI

US Military Academy 10996
ATTN: MAEN-A
ATTN: Facilities Engineer

Naval Facilities Engr Command
ATTN: Facilities Engr Command
Code 20Y AZ (2)

US Army CHPPM
ATTN: MCHB-DC-EEN (5)

US Gov’t Printing Office 20401
ATTN: Rec Sec;Deposit Sec (2)

Nat’! Institute of Standards & Tech
ATTN; Library 20899

Defense Construction Supply Center
ATTN: DCSC-WI 43216-500

Defense Tech Info Center 22060-6218
ATTN: DTIC-O (2)

49 + 12
5/97



dx? + dy? + dz? = h dé? + m2 dn? + hy dg? (4.4)

where the scale factors are defined in [12] as,
2 1/72 2 1/2
ne oa [St np od (+
e=a[ 7] baa a
G+1

=4((1 -e)(e? + 1)!" (4.5)

With the use of the expression for the Laplacian operator in orthogonal curvilinear

coordinates, the following differential equation is obtained [12].

2 pPe2Paos Ste Fy oly ge
amt Wag *5e & ete TMi oe +a(P+)iv=0 a6

where c = 4 kd. By the usual procedure of the separation of variables, the solution of (4.6)
is found as

Winn = Smn(-ic, 1) Rmn(-ic,i&) exp(imd) (4.7)

The angle and radial functions, Sma(-ic,n) and Rma(-ic,i€), satisfy the ordinary differential

equations [11]

ale +1) 2 San(- seit) Aann(-ic) - 02? - M2] Raa (-ic,i8) = (48)
0& 2 +1

a 2

an [0- re) 2 Sp Stiles) | + | nic) + om a] Smn(-ic,n) = 0 (4.9)

The angle and radial functions are the key functions in computing (4.2). All of the variables
in (4.2) are evaluated at m = 0, and the rest of this chapter assumes m = 0. Flammer
describes the derivations of the angle and radial functions in detail [11], and the following
two sections briefly explain his derivations since understanding the derivations are crucial in

developing a computer code.

42



Nje——Pr__ (r22) (4.15)

Ye- Aon - Nevo
and reciprocally
Nes2 = hon te, (t 22) (4.16)
is

with Nz = Aon - Yo. N2 = Aon - Yi-

The recursion formula (4.11) must be convergent, and this condition leads to the
iterative computations of N, from (4.14) for large values of r. Conversely from the initial
values Nz and N3, N, is calculated from (4.15) for small values of r. When c is small, the
dominant coefficient for given n is d9" [11]. Thus it is convenient to obtain the ratios of the
coefficients d"/d9". It is observed that the series of the coefficient ratios obtained from
(4.15) is accurate only for large r, and the series from (4.14) is accurate only for small r.
Therefore for r less than n, the coefficient ratios from (4.14) are taken, and for r greater than
n, the coefficient ratios from (4.15) are taken.

The calculation of the coefficients are based on the accurate value of Agn. The

eigenvalue for a small c? is obtained from a series in powers of c? [11].

Donic) = YE (1K 198 c2k (4.17)
Par) 2k
where 13" =n (n+l)
Balas 1

2 (2n -1) (2n +3)

re -(n +1)? (n + 2)? + n?(n- 1)?
* "2 @n- 1) Qn +2) (2n+3)3(2n+5) 2(2n-3)(2n-1)> n+ 1)

r= n?(n- 1)? + (n +1)? (n+ 2)
$ (2n - 5) (2n - 3) (2n- 1)° (2n +1) (2n+3) (2n- 1) (2n + 2) (2n + 3)5 (2n +5) (2n +7)

On Lpeleoyl
Ig 2A+TEB+eC+ 5D

A= n2(n- 1)2 7 (n + 1)? (n+ 2)2
(2n - 5)? (2n- 3) (2n- 1)? Qn +1) (2n+3)? (2n- 1)? (Qn +1) n+ 3)? (Qn +5) (2n+ 7)

44



acoustic source. Applying the radiation condition to (2.10) eliminates the field contribution

from the fictitious surface. Thus the volume integral is reduced to a surface integral:

9G) Gee) PPO) as, (2.11)

P(r) = P(r) “| [r« ) nr’
ss on on

where y(r) and w(r') are re-expressed as the pressure field P(r) at r and P(r’) at r'
respectively.

P(r), P(r’), and the partial derivative of P(r’) are the three unknowns to be solved.
The surface is assumed to be perfectly rigid. Rigid surfaces do not move in the presence of
acoustic pressure, thus the normal velocity field on the surfaces is zero. The pressure field on
the surface is doubled, and its partial derivative normal to the surface is zero. Finally the IE

of interest is
IG) ag

nm

P(r) = PiP(r) + { P(r') (2.12)
Ss

(2.12) will be used to build a matrix equation that can be solved in a computer program.

2.2 Evaluation Of The Singularity

The equation (2.12) can be evaluated at any observation point r, and when r
approaches the surface of the scatterer, the integral in (2.12) becomes singular at r=r’.
Morita, et. al. cleverly evaluate the integral at the singularity [3].

Take an infinitesimally small surface around r' and call it 5S. For a convenient
evaluation of the singularity, shift the coordinate system such that the observation point r and
the integration point r' are (0,0,z) and (p',9',z'), respectively (see Figure 2.2) and set z'=0.

The integral is then approximated as,

exp(ikir-r'| 1
{ P(r’ py SO") as Ll set 4qtr-r| lise oe Mazel as (2.13)
Ss on on on
8s és
7