‘flfm‘rptfiwnma

I

fm
1 LI
“1}

‘1‘.

m

. .Ir
..

.2 .
fit.”
._ )- x
u, .
.4

:35? “a.
,‘ml‘iillck

an»)

F 9..

firearmumé. .mw.

. gfifinfifihfi E... .n

.2... .zfiyfi?

w... 1..
.-

 

 

.:. 2.
. 0.911

 

was?
33....
2 ‘2.

: 3.4.an z.
. LS.

. I! ‘-

1‘.‘
1! .v‘ch!‘

unrmmzw}
.
33:3;
\ 1U} I’lfi‘
«3.93994.th

:u v
. «isrtixz'ngl
«5:3 .

#3:? i :1...
t. t a: I". .
.I... Q; '3’?
333...!

2).!
m2 1. I.‘!‘\....k
. cit: v}: vi!-

‘7th\:%!
z . .fi
23.2.8. 412...!

3.135.. ...e.,$!ll|
'3. iii): ‘1

.11)....
.533 ‘2‘

 

 

LIBRARY
University

 

 

 

This is to certify that the
dissertation entitled

Rapid Numerical Evaluation of Ultrasound Pressure Integrals
and Potential Integrals

presented by

Duo Chen

has been accepted towards fulfillment
of the requirements for the

PhD. degree in Electrical Engineering

 

 

2&4

 

 

 

Dee. (9, Zoo?

 

Date

MSU is an Affirmative Action/Equal Opportunity Employer

 

 

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

 

DATE DUE

DATE DUE DATE DUE

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

5108 KzlProleoCGPresIClRCIDateDue.indd

 

_- .fi

RAPID NUMERICAL EVALUATION OF ULTRASOUND
PRESSURE INTEGRALS AND POTENTIAL INTEGRALS

By

Duo Chen

A DISSERTATION

Submitted to
Michigan State University
in partial fulfillment of the requirements
for the degree of

DOCTOR OF PHILOSOPHY

Electrical Engineering

2009

ABSTRACT

RAPID NUMERICAL EVALUATION OF ULTRASOUND
PRESSURE INTEGRALS AND POTENTIAL INTEGRALS

By
Duo Chen

Analytical expressions are derived for fast calculations of time-harmonic and tran-
sient near field pressures generated by triangular pistons. These fast expressions re-
move singularities from the impulse response, thereby reducing the computation time
and the peak numerical error with a general formula that describes the nearfield pres-
sure produced by any triangular piston geometry. The timedomain expressions are
further accelerated by a time-space decomposition approach that analytically sepa-
rates the spatial and temporal components of the numerically computed transient
pressure. Analytical 2D integral expressions are derived for fast calculations of time-
harmonic and transient nearfield pressures generated by apodized rectangular pistons.
The 2D expressions eliminate the numerical singularities that are otherwise present
in numerical models of pressure fields generated by apodized rectangular pistons. A
simplified time space decomposition method is also described, and this method fur-
ther reduces the computation time for transient pressure fields. The results, compared
with the Rayleigh-Sommerfeld integral, the Field II program and the impulse response
method, indicate that the FN M achieves smallest errors for the same computation
time among those methods. A 1D FN M for calculating the pressure generated by a
polynomial apodized rectangular piston is also obtained. The fast method is based
on the instantaneous impulse response. A trigonometric transform of the integrand is
performed and the order of integration is exchanged to obtain the 1D integral for the
apodized F NM for both apodization functions. The time and error comparisons are

performed among the 1D polynomial apodized FNM, the 2D apodized FNM and the

Rayleigh-Sommerfeld integral. The results show that the 1D polynomial apodized
FNM has the fastest convergence. Analytical expressions are derived for fast cal-
culations of potential integrals. These potential integrals inculde uniformly excited
volume potential integrals, polynomial apodized surface integrals and polynomial
apodized volume potential integrals. The derivation starts with the fast near—field
method (FNM), which originates from ultrasound pressure calculations generated by
polygonal pistons. For potential integrals evaluated over a volume source, the volume
source is first subdivided into subdomains about the observation points. The total po—
tential is the summation of the potential over each submain which can be reduced to
1D integrals. Those calculation methods remove the singularities from the Rayleigh-
Sommerfeld integral by subtracting sigularities in the integrand and thus can achieve
rapid convergence. Simulations results are compared with the Rayleigh-Sommerfeld
integral and the singularity cancellation method evaluated on a 3D grid. The results
indicate that the 1D FN M expressions reduces the computation time or the number
of sample point needed significantly than the Rayleigh-Sommerfeld integral and the

singularity cancellation method for a given number signifcant digits.

Dedicated to my wife, Jie Sun, and my son, Eric Yaoting Chen.

iv

ACKNOWLEDGMENTS

Firstly, I would like to thank my Ph.D. advisor Dr. Robert J. McGough for his
great support during the past five and a half years. I have learned a lot from his
guidance, his encouragement and his insightful discussions. Secondly, I would like
to thank my Ph.D. Committee members for providing helpful guidance and discus—
sions: Dr. Shanker Balasubramaniam, Dr. Edward Rothwell, Dr. Mei Zhuang,
and Dr. Neil Wright. Thirdly, I would like to thank my colleagues at the Biomed-
ical Ultrasonics and Eletromagnetics Lab (BUEL) at MSU for their friendship and
help during my stay at MSU: Ruihua Ding, Liyong Wu, Xiaozheng Zeng, Khawar
Khurshid, Don Chorman, Jason Biel, James F. Kelly, Josh Wong, Don VanderLaan,
Matthew Richard Jennings, and Christopher Johnson. Finally, I would like to express
my deepest gratitude and love to my dear wife, Jie Sun for her unconditional love

and support to me. Without her I would have never completed this achievement.

TABLE OF CONTENTS

LIST OF TABLES ................................

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

1 Introduction ..................................

1.1
1.2
1.3

Ultrasound pressure calculations .....................
Potential integrals .............................
Thesis Content ..............................

2 A Fast Nearfield Method for Calculations of Time-harmonic and
Transient Pressures Produced by Triangular Pistons ........

2.1

2.2

2.3

2.4

Time-harmonic and Transient Nearfield Pressure Calculations for Tii-
angular Sources ..............................
2.1.1 Impulse response calculations for a triangular source .....
2.1.2 The fast nearfield method for a triangular source .......
2.1.3 Superposition calculations with impulse response and F NM ex-

pressions ..............................
2.1.4 Transient input waveform ....................
2.1.5 Error calculations .........................
Results ...................................
2.2.1 Time-harmonic nearfield pressure calculations .........
2.2.2 Transient nearfield pressure calculations ............
Discussion .................................
2.3.1 Time and error calculations ...................
2.3.2 Advantages of the FNM for time-harmonic and transient cal-

culations ..............................
2.3.3 Field II ..............................
2.3.4 Smoothed impulse response ...................
Conclusion .................................

3 A 2D Fast Nearfield Method for Apodized Rectangular Pistons .

3.1

3.2

Existing calculation methods .......................
3.1.1 The Rayleigh-Sommerfeld integral ...............
3.1.2 The Field II program .......................
Fast nearfield method for apodized rectangular pistons ........
3.2.1 Steady-state Apodized F NM Expression ............
3.2.2 Transient Apodized FNM Expression ..............
3.2.3 Apodization function .......................

vi

X

xiv

“NH...

6
6
11

14
16
17
17
18
22
26
26

27
28
29
29

33
34
34
35

3.2.4 Input transient pulse ....................... 46

3.2.5 Time space decomposition .................... 46
3.2.6 Error Calculations ........................ 47
3.3 Results ................................... 50
3.3.1 Time-harmonic pressure calculations .............. 50
3.3.2 Transient field calculations .................... 55
3.4 Discussion ................................. 60
3.4.1 Large-scale computation ..................... 60
3.4.2 Time and error comparisons ................... 62
3.4.3 Apodization functions ...................... 62
3.4.4 Time Space Decomposition .................... 63

A 1D Fast Nearfield Method for Rectangular Pistons with Poly-

nomial Apodization ............................. 65
4.1 Polynomial apodization derivation .................... 66
4.1.1 Instantaneous Impulse Response ................. 66
4.1.2 Time-harmonic pressure calculations .............. 67
4.1.3 1D quadratic apodization .................... 69
4.1.4 2D Quadratic apodization .................... 72
4.1.5 The 2D apodized FNM ...................... 75
4.1.6 The Rayleigh-Sommerfeld integral ................ 76
4.1.7 Error Calculations ........................ 76
4.2 Results .................................. 77
4.2.1 1D quadratic apodization .................... 77
4.2.2 2D Apodization function ..................... 81
4.3 Discussion ................................. 86
4.3.1 Advantages and disadvantages .................. 86
4.3.2 Interpolation of the apodization function ............ 86
4.4 Conclusion ................................. 87

A Fast Nearfield Method for the Numerical Evaluation of 3D Po-

tential Integrals ............................... 88
5.1 F NM Calculations for Rectangular and Triangular Sources ...... 90
5.1.1 The potential integral ...................... 90
5.1.2 FNM calculations for a planar source .............. 90
5.1.3 FNM calculations for a volume source .............. 93
5.1.4 Error calculations ......................... 97
5.2 Results ................................... 98
5.2.1 Comparisons of potential evaluated at single points ...... 98
5.2.2 The potential evaluated on a 3D grid .............. 103
5.3 Discussion ................................. 107

vii

5.3.1 Other geometries ......................... 107

5.3.2 Sample point calculations .................... 108
5.3.3 Future work ............................ 108
5.4 Conclusion ................................. 109

6 A Fast Nearfield Method for Numerical Evaluation of Surface In-

tegrals with Polynomial Apodization .................. 110
6.1 Method .................................. 111
6.1.1 Uniformly excited source ..................... 112
6.1.2 Linear apodization ........................ 113
6.1.3 Quadratic apodization ...................... 114
6.1.4 Cubic apodization ........................ 115
6.1.5 Error calculations ......................... 117
6.2 Results ................................... 117
6.2.1 Linear apodization ........................ 118
6.2.2 Quadratic apodization ...................... 119
6.2.3 Cubic apodization ........................ 123
6.3 Discussions ................................ 124
6.3.1 Other Polygonal Sources ..................... 124
6.3.2 Higher order polynomials ..................... 126
6.4 Conclusion ................................. 127

7 A Fast Nearfield Method for Numerical Evaluation of Volume In-

tegrals with Polynomial Apodization .................. 128
7.1 Method .................................. 129
7.1.1 Uniformly excited planar source ................. 131
7.1.2 Linear apodization for a planar source ............. 133
7.1.3 Linear apodization for a subdomain ............... 135
7.1.4 Linear apodization in the y direction .............. 137
7.1.5 Linear apodization in the z direction .............. 139
7.1.6 The prism geometry ....................... 139
7.1.7 Global and local systems ..................... 140
7.1.8 Evaluating potential integrals generated by each subdomain . 145
7.1.9 Error calculations ......................... 146

7.2 Results ................................... 147
7.3 Discussion ................................. 148
7.3.1 Other volume sources ....................... 148
7.3.2 Higher order polynomials ..................... 149

7.4 Conclusion ................................. 149
8 Conclusion .................................. 151

viii

BIBLIOGRAPHY ................................ 153

ix

2.1

2.2

2.3

3.1

3.2

3.3

LIST OF TABLES

Basis functions for time-space decomposition with a Hamming—weighted
pulse. ...................................

Number of Gauss abscissas, computation times, and computation times
relative to the F NM that describe the reduction in the computation
time achieved with the fast nearfield method relative to the impulse
response and methods that approximate the impulse response for peak
errors of 10% and 1%. The FNM and exact impulse response results
are evaluated for time-harmonic calculations on a 81 by 101 point grid
located in the :1: = 0 plane, and the Field II and smoothed impulse
response results are evaluated on an 81 by 86 point grid in the a: = 0
plane that is slightly offset from the transducer face. (a) For a 10%
peak error and (b) for a 1% peak error. ................

Comparisons of computation times, input parameters, and computa-
tion times relative to the FNM that describe the reduction in the
computation time achieved with the FNM and time—space decompo-
sition relative to the exact and approximate impulse response for spec-
ified maximum errors of 10% and 1%. For FNM, impulse response,
Field 11 calculations with ‘use-triangles’, and Field II calculations with
’useJectangles’, these transient results are evaluated in an 81 by 101
spatial point by 85 time point grid, and for the smoothed impulse
response, the results are valued at the same temporal points in a re-
stricted 81 by 86 point spatial grid. (a) For a 10% peak error and (b)
for a 1% peak error. ...........................

Terms that define the time—space decomposition of the Hanning-
weighted pulse v(t — 'r) for transient apodized FNM calculations.

Terms that define the time-space decomposition of the derivative of
a Hanning-weighted pulse i)(t - T) for transient calculations with the
apodized Rayleigh-Sommerfeld integral ..................

Simulation parameters for time—harmonic calculations that achieve nor-
malized root mean square error (NRMSE) values of 0.1 and 0.01. Pa-
rameters listed include the number of Gauss abscissas or the corre-
sponding Field II parameters, the resulting computation time, and
computation time relative to the apodized F NM for the Rayleigh in-
tegral and the Field II program. (a) For a 0.1 NRMSE and (b) for a
0.01 NRMSE. ...............................

13

31

32

47

48

3.4

4.1

4.2

5.1
5.2

5.3

5.4

Simulation parameters for transient calculations that achieve normal-
ized root mean square error (N RMSE) values of 0.1 and 0.01. Parame-
ters listed include the number of Gauss abscissas or the corresponding
Field 11 parameters, the resulting computation time, and the compu-
tation time relative to the apodized F NM for the Rayleigh integral and

the Field II program. (a) For a 0.1 NRMSE and (b) for a 0.01 NRMSE. 61

Simulation parameters that achieve peak normalized error values of
10%, 1%, and 0.1% for the 1D quadratic apodization function u(x) =
u2 — 2:2. Parameters listed include the number of Gauss abscissas, the
computation time, and the computation time relative to the polynomial
apodized FNM for the Rayleigh-Sommerfeld integral. .........

Simulation parameters that achieve peak normalized error values of
10%, 1%, and 0.1% for the 2D apodization function u(:c,y) = (a:2 —-
r12)(y2 — a2), where a = 2 wavelengths. Parameters listed include the
number of Gauss abscissas, computation time, and the ratio of the
computation time relative to the polynomial apodized FNM. (a) 10%
peak normalized error, (b) 1% peak normalized error and (c) 0.1% peak
normalized error. .............................

The geometry of the prism. .......................

Reference potential fields are evaluated at three points (3:, y, z) =
(1/3m, 1/3m, d), where d = 0.5m, 1.0m and 1.25m. The reference
fields are computed using the singularity cancellation method with 30,
30, and 30 abscissas in the :r, y, and 2 directions, respectively.

Simulation parameters that achieve between 2 to 5 significant digits
with the FNM and the singularity cancellation method for the potential
evaluated on the observation points (x, y, z) = (1 /3, 1/3, d), where
d = 0.5, 1.0, and 1.25 [m]. Parameters listed include the number of
sample points for each observation point and the ratios of the number
of sample points required to achieve a specific accuracy relative to the
number required with the FNM. (a) Significant digits 2 and 3 and (b)
Significant digits 4 and 5. ........................

Reference potentials evaluated at three points (2:, y, z) =
(1/3m, 1/3m, d), where d = 0.1m, 0.01m, and 0.0001m. The refer-
ence potentials are computed using the singularity cancellation method

with 30, 30, and 30 abscissas in the 2:, y and 2 directions, respectively.

xi

77

82

98

100

102

103

5.5

5.6

6.1

6.2

6.3

7.1

Simulation parameters that achieve between 2 to 5 significant digits
with the uniformly excited 1D FNM expression for the volume integral
and the singularity cancellation method for the potential evaluated
at the observation points (:13, y, z) = (1/3m, 1/3m, d), where d =
0.1m, 0.01m, and 0.0001m. Parameters listed include the number of
sample points required at each observation point and the number of
sample points required to achieve a specified error relative to the FNM.
(a) Significant digits 2 and 3. (b) Significant digits 4 and 5 .......

Simulation parameters that achieve between 2 and 5 significant digits
with the uniformly excited 1D F NM expression for the volume integral
and the singularity cancellation method for the potential evaluated on
a 3D grid. Parameters listed include the number of sample points for
each observation point and the number of sample points required to
achieve a specified error relative to the uniformly excited 1D FNM
expression for the volume integral. (a) Significant digits 2 and 3. (b)
Significant digits 4 and 5. ........................

Simulation parameters that achieve 10%, 1% and 0.1% peak normalized
error in the potential obtained with the fast nearfield method and the
Rayleigh-Sommerfeld integral computed using linear apodization for
a triangular source. Parameters listed include the number of Gauss
abscissas, the computation time, and the ratios of the computation
times relative to the times obtained with the fast nearfield method.

Simulation parameters that achieve 10%, 1%, and 0.1% peak nor-
malized error in the potential obtained with the fast nearfield
nearfield method and the Rayleigh-Sommerfeld integral computed us-
ing quadratic apodization for a triangular source. Parameters listed
include the number of Gauss abscissas, the computation time, and the
ratios of the computation times relative to the times obtained with the
fast nearfield method ............................

Simulation parameters that achieve 10%, 1% and 0.1% peak normalized
error in the potential obtained with the fast nearfield method and the
Rayleigh-Sommerfeld integral computed using cubic apodization for
a triangular source. Parameters listed include the number of Gauss
abscissas, the computation time, and the ratios of the computation
times relative to the times obtained with the fast nearfield method.

Vertex locations defined for the prism source. .............

105

107

118

119

124

140

7.2 Simulation parameters that achieve between 2 and 5 significant digits
using the FNM and the singularity cancellation method to evaluate
a potential integral with linear apodization evaluated on a 3D grid.
Parameters listed include the number of significant digits achieved and
the ratio of the computation time relative to the time required for the
FNM. (a) Significant digits 2 and 3. (b) Significant digits 4 and 5. . . 147

xiii

2.1

2.2

2.3

2.4

2.5

LIST OF FIGURES

Triangular source geometries defined for nearfield pressure calculations.
The nearfield pressure is evaluated above the vertex A (indicated in
bold), and the shape of the triangle (right, acute, or obtuse) is defined
by the angle ABCA. The height of each triangle is indicated by l,
and the bases of the individual right triangles are indicated by s, 31,
and 32. The acute triangle in b) is represented by the sum of two
right triangles, and the obtuse triangle in c) is defined as the difference
between two right triangles. .......................

Superposition operations that calculate nearfield pressures generated
by an equilateral triangular source, where each side is four wavelengths
long. The vertex D (indicated in bold) is the projection of the observa-
tion point onto the source plane, which partitions the radiating source
into three triangles with sides (ai, bi: Ci). (a) The field point is located
inside of the equilateral triangular source. (b) The field point is located
outside of the equilateral triangular source. ..............

Simulated time-harmonic pressure field in the a: = 0 plane for an equi-
lateral triangular source with sides equal to 4 wavelengths. The refer-
ence field is generated by the impulse response method computed with
100,000—point Gauss quadrature. ....................

Peak normalized error for calculations of nearfield pressures generated
by the triangular source in Figure 2.2 plotted as a function of the com-
putation time. The results show that the FN M consistently achieves
smaller errors in less time than exact and approximate impulse response
calculations for time-harmonic excitations. ..............

Simulated transient pressure field in the a: = 0 plane for an equilateral
triangular source with sides equal to 4 wavelengths. For this calcula-
tion, the excitation is the Hanning—weighted pulse in Eq. 2.17, and the
transient pressure is evaluated at 85 time points in an 81 by 101 point
grid. The result is plotted at 1.8125ps after the initiation of the input
pulse. ...................................

xiv

2.6

3.1

3.2

3.3

3.4

3.5

3.6

The peak normalized error plotted as a function of the computation
time for the FNM/time—space decomposition method, the impulse re-
sponse method, and methods that approximate the impulse response.
These errors and times are evaluated for transient nearfield calculations
of an equilateral triangular source with sides equal to 4 wavelengths.
The excitation for these calculations is a Hanning-weighted pulse with
a center frequency of 2MHz. .......................

Orientation of the computational grid relative to the rectangular
source. The rectangular source, which has width 0. and height b, lies
in the z = 0 plane. The dashed lines define the extent of the compu-
tational grid in the :r = a/ 2 plane. The extent of the computational
grid is 2b by 0.99a2/4A in the y and 2: directions, respectively.

The decomposition of an apodized rectangular source into smaller rect-
angles, where each small rectangle is Au wide and AV high. The
apodization function f (p, v) is defined as constant over each small
rectangle. .................................

The apodization function f (p, U) = sin(mr/a) sin(mr/b) evaluated on
the face of a 4A by 4/\ square piston. The maximum value of the
apodization function is achieved when u = 2/\ and V = 2A. ......

Simulated reference pressure field generated by an apodized rectan-
gular source with each side equal to 4 wavelengths. The results are
evaluated in the a: = 2.0/\ plane for a time-harmonic excitation.

The normalized error distribution r}(:r, y, z; It) describes the difference
between the reference pressure field and the computed pressure field for
an apodized 4A by 4A source. The error distribution 77 is plotted for a)
the apodized F N M evaluated with 16-point Gauss quadrature in each
direction, b) the apodized Rayleigh integral evaluated with 16-point
Gauss quadrature in each direction. ..................

The normalized error distribution 11(23, y, z; It) describes the difference
between the reference pressure field and the computed pressure field
for an apodized 4A by 4A source. The error distribution 17 is plotted for
the Field 11 program evaluated with f3 = 48MHz and 30 subdivisions
in each direction. .............................

XV

23

36

51

52

3.7

3.8

3.9

4.1

4.2

4.3

4.4

Normalized root mean square error (NRMSE) plotted as a function of
the computation time for time-harmonic calculations with the apodized
FNM, the apodized Rayleigh-Sommerfeld integral, and the Field II
program. This figure demonstrates that the apodized FNM achieves
the smallest errors for a given computation time, and the apodized
F NM requires the smallest amount of time to achieve a given error
value. ...................................

Simulated reference transient field for an apodized square source ex~
cited by the Hanning-weighted pulse in Eq. (3.23) with f0 = 1.5 MHz
and W = 2.0).. The sides of the square source are equal to 4A. The
apodization function is given by Eq. (3.22). The transient reference
pressure, evaluated in the :r = 2.0A plane, is computed with 100,000
Gauss abscissas in each direction using the Rayleigh integral. Results
are plotted at a) t = 1.5625ps and b) t = 3.0625ps. ..........

Normalized root mean square error (NRMSE) plotted as a function of
the computation time for transient pressure calculations evaluated with
the apodized FNM, the apodized Rayleigh-Sommerfeld integral, and
the Field II program. For the same computation time, the apodized
F NM achieves the smallest errors, and for the same error, the apodized
F NM requires the least amount of time. ................

Normalized simulated time-harmonic reference pressure field in the
y = 0 plane for a rectangular source with each side equal to 4 wave-
lengths. The pressure field is computed with 10,000 Gauss abscissas
in each direction using the Rayleigh-Sommerfeld integral with the 1D

quadratic apodization u(a:) = 2:2 —— a2. .................

a) Maximum errors and b) computation times for the polynomial
apodized FN M and the Rayleigh-Sommerfeld integral. ........

Time vs. error comparison between the polynomial apodized FN M and
the Rayleigh-Sommerfeld integral. For the same computation time, the
polynomial apodized F N M achieves smaller errors, and for the same
error, the polynomial apodized FNM requires less time. .......

Absolute value of the simulated time—harmonic reference pressure field
in the y = 0 plane for a square source with each side equal to 4
wavelengths. The apodization function is a 2D function given by
u(a:,y) = (:1:2 — a2)(y2 — a2), where a =2 wavelengths. The pres-
sure field is computed with 10, 000 Gauss abscissas in each direction
using the Rayleigh-Sommerfeld integral. ................

xvi

54

58

59

78

79

80

4.5

4.6

5.1

5.2

5.3

5.4

5.5

5.6

5.7

6.1

a) Maximum errors and b) computation times plotted as a function
of the number of Gauss abscissas for a 2D apodization function with
the 1D polynomial apodized FNM, the 2D apodized FNM and the
Rayleigh-Sommerfeld integral. .....................

Time vs. error comparison between the 1D polynomial apodized FNM,
2D apodized FNM and the Rayleigh-Sommerfeld integral. The poly-
nomial apodized F N M converges to smallest errors with least time.

Parameters defined for FNM calculations with a triangular source
AABC. The lengths a,, bi: and c,- are defined for three different tri-
angles that share the vertex D in each subfigure. a). The projection
of the observation point D is located inside of the triangular source.
b). The projection of the observation point D is located outside of the
triangular source. ............................

Subdividing a prism about an observation point to form subtetrahe-
drons and subpyramids. (a). Subtetrahedron. (b). Subpyramid.

The geometric configuration indicating how the potential integral is
evaluated for a subtetrahedron. .....................

The prism geometry, where vertex B is coincident with the origin, and
vertices C, A, and E are located on the x, y, and z axes, respectively.

The number of significant digits in the computed potential for a prism
plotted as a function of the number of sample points for three observa-
tion points (x, y, z) = (1/3m, 1/3m, d), where d = 0.5m, 1.0m, and
1.25m. ..................................

The number of significant digits in the calculated potential achieved for
the prism in Figure 5.4 plotted as a function of the number of sample
points for three observation points (at, y, z) = (1 /3m, 1/3m, d), where
d = 0.1m, 0.01m, and 0.0001m. .....................

The number of significant digits achieved in calculations of the po—
tential of the prism shown in Figure 5.4 plotted as a function of the
number of sample points evaluated over a 3D grid. The results Show
that the FNM is accurate to a large number of significant digits than
the singularity cancellation method for the same number of sample
points. ..................................

a) Maximum errors and b) computation times obtained with the fast

nearfield method and the Rayleigh-Sommerfeld integral for a triangular
source with linear apodization. .....................

xvii

84

85

92

93

94

99

101

104

106

6.2

6.3

6.4

6.5

6.6

7.1

7.2

7.3

7.4
7.5
7.6
7.7

Numerical errors plotted as a function of the computation time for a
triangular source with linear apodization. The results show that the
fast nearfield method achieves much better convergence performance
than the Rayleigh-Sommerfeld integral. ................

a) Maximum errors and b) computation times obtained with the fast
nearfield method and the Rayleigh-Sommerfeld integral for a triangular
source with quadratic apodization. ...................

Numerical errors plotted as a function of the computation time for a
triangular source with quadratic apodization. The results show that
the fast nearfield method achieves much better convergence than the
Rayleigh-Sommerfeld integral. .....................

a) Maximum errors and b) computation times obtained with the fast
nearfield method and the Rayleigh-Sommerfeld integral for a triangular
source with cubic apodization. .....................

Numerical errors plotted as a function of the computation time for a
triangular source with cubic apodization. The results show that the
fast method achieves much better convergence performance than the
Rayleigh-Sommerfeld integral. .....................

Subdomains defined for a prism where the shared vertex is defined at
the observation point ............................

Geometric parameters defined for the potential calculations in a tetra-
hedral subdomain. ............................

The prism geometry, where vertex B is coincident with the origin, and
vertices C, A, and E are located on the x, y, and z axes, respectively.

The local coordinate system for the plane ADEB ............
The local coordinate system for the plane ADF C ............
The local coordinate system for the plane BEF C ............

The number of significant digits achieved in calculations 0d the poten-
tial integrals linear apodization over the prism plotted as a function
of computation time for evaluations on a 3D grid. The results show
that the FNM achieves a larger number of significant digits than the
singularity cancellation method for the same computation time.

xviii

121

122

123

125

126

130

132

141
142
143
145

148

CHAPTER 1

Introduction

1.1 Ultrasound pressure calculations

As increasingly complex transducer geometries are adopted for emerging applications
of ultrasound imaging and therapy, new methods are needed for rapid calculations
of pressure fields produced by these transducers. Fast numerical simulations are es-
pecially important for simulations of phased array structures containing hundreds
or thousands of tranducers that generate pressure in large computational domains.
Simulation methods that either directly evaluate the impulse response [1, 2, 3, 4]
or subdivide each transducer into smaller elements and then superpose the pressure
[5, 6] provide a convenient model for these calculations, but in the nearfield region,
the convergence of these methods is relatively slow [7] Furthermore, numerical im-
plementations of the impulse response for flat unfocused transducers encounter some
difficulties throughout the paraxial region [7, 8, 9].

Pressure fields from uniform planar sources can be calculated using several numer—
ical methods including the impulse response method [1, 2, 3], the Field 11 program
[10, 11], the Rayleigh-Sommerfeld integral [6], and the fast nearfield method (FN M)

[7, 8]. Among these methods, the FNM achieves the smallest error in the least time

[7, 8, 12]. The F NM for uniformly excited rectangular pistons [8] and circular pistons
[7] eliminates the 1/R singularity and therefore avoids problems with large errors in
the nearfield region. The F NM also converges rapidly in the nearfield, which has
been demonstrated for both time-harmonic and transient excitations [7, 8, 36]. For
uniform excitations, the FNM is ideal for nearfield pressure calculations and for ref-
erence calculations that evaluate the numerical errors associated with other methods
[13, 14].

Modeling the nearfield pressure generated by a spatially varying particle velocity
on the face of a rectangular piston is of practical importance for many acoustics
applications [15, 16], where accurate and fast computer simulations are required. Fast
ultrasound simulations are especially important for calculating pressure fields in large
computational domains, specifically steady-state HIFU simulations for ultrasound
therapy [17, 18] and transient calculations for ultrasound imaging [19, 20]. Individual
array elements are typically modeled as single baffled rectangular sources with uniform
surface particle velocities, but the particle velocity on the transducer face is in general
nonuniform. For example, Lin et. a1. [21] shows that the surface particle velocity of
a fluid-loaded piezoelectric element on a phased array is nonuniform, and Borges et.
al [22] demonstrates that the surface velocity of a single array element is apodized

and delayed when matching layers are used.

1.2 Potential integrals

Numerical calculation of potential integrals involved in integral equations are of
great importance in scattering problems. The polynomial function is a very pop-
ular function to approximate the electric and magnetic currents. For example, three—
dimensional polynomials are adopted to approximate the entire-domain normalized

current density by Moraros and Popovic [23] as applied to the optimization of volume

potential integrals involved in the moment-method analysis of 3D dielectric scatters.
The commonly used Rao-Wilton-Glisson basis that is also a linear basis function
which is used to approximate both the magnetic and electric currents [24]. P0-
tential integrals are often singular and direct evaluation of potential integrals may
encounter numerical difliculties. To improve the accuracy of the potential integrals,
singularity subtraction methods [23, 25, 26, 27, 28] or singularity cancellation meth-
ods [30, 31, 32, 33] are often used to acheive better performance. Typically, those
methods manipulate the integrands to eliminate singularties and thus the number
of dimensions over which the integration is performed remains the same for both
methods. This approach is reasonable for general potential integrals; however, more
efficient expressions can be acheived when polynomial apodized potential integrals
are considered.

Potentials generated by uniformly excited plannar source and uniformly excited
plannar source are computed with both the Rayleigh-Sommerfeld integral [6] and the
fast nearfield method (FN M) [7, 8, 12]. The existing FNM achieves much better accu-
racy by eliminating the 1/R singularity in the Rayleigh-Sommerfeld integral and by
simplifying the multiple integral into a single integral. Since the numerical evaluation
of potential integrals with polynomial apodization are routinely evaluated, a similar

fast nearfield method is needed for those integrals to improve the performance.

1 .3 Thesis Content

This dissertation investigates fast pressure calculation methods for planar pistons
and fast potential evaluation methods for potential integrals. Chapter 2 introduces
fast calculations for the time—harmonic and transient pressures generated by tran-
gular pistons. The transient calculations are further accelerated by the time-space

decomposition method. Analytically 2D integral expressions for fast calculations of

time-harmonic and transient nearfield pressures generated by apodized rectangular
pistons are derived in Chapter 3. As a special case of Chapter 3, fast caculations of
pressure generated by a polynomial apodized rectangular pistion are obtained based
on the instantaneous impulse response in Chapter 4. Chaper 5 introduces 1D fast
expressions for calculations of uniformly excited volume potential integrals that are
otherwise represented by a triple integral. Fast 1D calculation expressions for pres-
sures generated by surface integrals and volume integrals with polynomial apodization

are introduced in Chapters 6 and 7. Chapter 8 concludes the thesis.

CHAPTER 2

A Fast Nearfield Method for

Calculations of Time-harmonic and
Transient Pressures Produced by

Triangular Pistons

The substantial reduction in computation time demonstrated by the FN M for calcu-
lations of nearfield pressures generated by circular and rectangular pistons motivates
the derivation of similar integral expressions for triangular sources. After the impulse
response is obtained for right, acute, and obtuse triangular sources, general FNM
expressions for time-harmonic and transient inputs are then demonstrated for a tri-
angular source, and the time-space decomposition of the F N M integral is presented for
a transient excitation. Based on these expressions for the nearfield pressure generated
by a triangular source, computation times are evaluated for the same peak numerical
errors. For time-harmonic inputs applied to a triangular source, results show that
FNM calculations are several times faster than both exact and approximate impulse

response calculations, and for pulsed excitations, results demonstrate that FNM cal-

culations performed with time-space decomposition are also much faster than exact

and approximate impulse response calculations for triangular piston geometries.

2.1 Time-harmonic and Transient Nearfield Pres-

sure Calculations for Triangular Sources

2.1.1 Impulse response calculations for a triangular source

The geometry for a right triangular source with a right angle ABCA at vertex C is
depicted in Figure 2.1a. For this right triangle and the triangles in Figures 2.1b and
2.1c, the impulse response is evaluated at a point directly over the vertex A (indi-
cated in bold in Figure 2.1), where the the orthogonal projection of the observation
point onto the source plane is exactly coincident with the vertex A, and the distance
from the observation point to the source plane along this orthogonal projection is
represented by the variable 2. In Figure 2.1a, the acute angle [CAB = tan—1(s/l)
defines the angular extent of sector EAB with radius m, which has an im-
pulse response of c/(27r) tan-"1(s/l) for (23/0) S t S Wk. The impulse
response for the right triangle AABC contained within the sector EAB is obtained
by subtracting the impulse response of the region ECB between the curved outer
edge of the sector and the near edge of the right triangle so that only the contri-
bution from the right triangle AABC remains. The impulse response of the region
ECB is c/(27r)cos-1(l/m) for Wk S t S MR, and

therefore the impulse response at an orthogonal distance 2 above the vertex A is

c/(27r) tan—1(s/l) for t1 3 t 3 t2

h- (z;t)= c 27r {tan_1 s/l —cos—1—l———} fort gtgt ,
rzght /( l ( ) c2t2—z2 2 3

0 otherwise
(2.1)

   

 

\
N
[0
TI
CI)
H10

 

 

A

 

 

 

(c) Triangular source with obtuse ABCA

Figure 2.1. Triangular source geometries defined for nearfield pressure calculations.
The nearfield pressure is evaluated above the vertex A (indicated in bold), and the
shape of the triangle (right, acute, or obtuse) is defined by the angle ABCA. The
height of each triangle is indicated by l, and the bases of the individual right triangles
are indicated by s, 31, and 32. The acute triangle in b) is represented by the sum of
two right triangles, and the obtuse triangle in c) is defined as the difference between
two right triangles.

where the values of t1, t2, and 253 are z/c, Wk, and Wk, respec-
tively.

For other triangular sources, the impulse response is readily constructed from the
sum or difference between two right triangles. Figure 2.1b contains an example of
a triangular source with an acute angle ABCA at vertex C. The expression for
the impulse response evaluated at a point directly over the vertex A is obtained by
evaluating the sum of the contributions from the right triangles ACDA and ABDA,
each with a right angle at vertex D. The resulting impulse response above the vertex

A in Figure 2.1b is represented by

hsum(z;t) =

c/(27r) {tan—1(sl/l) + tan—1(32/z)} for t1 3 t g t2
2% {tan—1(31/l) + tan—1(52/l) — 2cos—1 (l/m)} for t2 3 t 5 t3
c/(27r) {tan—1(sl/l) — cos-1 (HA-22:77)} for t3 _<_ t 3 t4

0 otherwise

3

(2.2)
where the values of 31 and 32 are selected such that .91 2 32 and the values of t1, t2,
t3 and t4 are z/c, Wk, 2:2 +12 + sg/c, and z2 + l2 + Sig/c, respectively.

Similarly, Figure 2.1c contains an example of a triangular source with an obtuse
angle ABCA at vertex C, where the impulse response is again evaluated at a point
directly over the vertex A, but the impulse response is instead evaluated for the

difference between two right triangles. The impulse response for the triangle in Figure

2.1c is

c/(27r) {tan—1(sl/l) — tan—1(32/l)} for t1 S t _<_ t2
hdz'ff(3;t) = c/(27r) {tan—1(sl/l) — cos‘1(l/\/c2t2 — 22)} for t2 3 t 3 t3 1
0 otherwise
(2.3)

where the values of 31 and .32 are selected such that 31 > 32 and the values of t1, t2,

and t3 are z/c, 22 + l2 + sg/c, and 22 + l2 + s%/c, respectively.

Time-harmonic impulse response calculations

The time-harmonic pressure generated by these triangular source geometries is pro-
portional to the Fourier transform of the impulse response. Therefore, the formula

for the time-harmonic pressure generated by the right triangle in Figure 2.1a is

,jwt - -. . . 2 2 2
cup 18 J — 3 — z .7 _ S _ ,/
Fright(zik)=—9T{Etan 176 ‘7" —Etan l-l—e 3k 2 +l +3
,/ 2 2 2
2 +1 +3 . l
—'kfi —1
+//——— 6 J cos ———dfi . (2.4)
122-+12 ,82—22

The time-harmonic pressures produced the remaining triangles depicted in F ig-
ures 2.1b and 2.1c are obtained by adding and subtracting the contributions of two
right triangles, as for calculations of the impulse response in Eqs. 2.2 and 2.3, respec-

tively.

Transient impulse response calculations

Transient nearfield pressures are computed with the impulse response through the

convolution
pa; t) = 90130) a he; a (2.5)
where the time derivative of the particle velocity t)(t) is evaluated analytically from

the excitation pulse v(t), and the convolution ® is evaluated with the fast Fourier

transform (F F T). In particular, the discrete Fourier transforms of 13(t) and h(z, t) are
computed with the FFT, the results are multiplied, and the inverse FFT is applied to
the product. The forward and inverse FFT routines are computed with the Fastest

Fourier Dansform in the West (FFTW) library [34].

Field II

Field II is a software package [10] that computes the impulse response either by su-
perposing far field contributions from small rectangles or by evaluating expressions
similar to Eqs. 2.1, 2.2, and 2.3. With both approaches, Field II modifies the impulse
response according to the area under the impulse response curve [11]. This modifica-
tion allows Field 11 to reduce the temporal sampling for impulse response calculations,
which are directly applicable to transient and steady—state nearfield pressure compu-

tations.

Smoothed impulse response

The Fourier transform of Eq. 2.5 is

P(z; w) = iwp0V(w)H(z,w). (2.6)

Normally, the excitation V(w) is bandlimited, so the high-frequency components in
the Fourier transform H (z; w) of the impulse response are negligible. To exploit the
bandlimited characteristics of the excitation v(t), the formula for a smoothed impulse

response is given by [35]

_ O(c(t + At/2)) —— O(c(t —- list/2))
”smooth“ 0 : 27rtcAt ’

 

(2.7)

where 0(ct) is the area that formed by the intersection of the transducer and the

sphere with radius ct centered at the observation point, and At is the length of the

10

rectangular pulse that smooths the analytical impulse response. The constraint

1
fEmax < E (2.8)

insures proper smoothing, where f Emaa: is the highest-frequency component of the
excitation pulse, and At = 0.02113 in the simulations that follow. The result obtained
from Eq. 2.7 is then directly applied to calculations of the nearfield pressure for time-

harmonic and transient inputs.

2.1.2 The fast nearfield method for a triangular source

Integral expressions that describe the fast nearfield method (FNM) for a triangular
source excited by a time-harmonic input are obtained by replacing the inverse cosine
term with the integral form of the inverse tangent and then exchanging the order of
integration in the impulse response expressions for right, acute, and obtuse triangles.

The procedure is illustrated by:

,/ 2 2 2

Z +l +8 . l

/\/—2—2 e—lkflcos_1—————dfi
2 +1 32 —z2

2 2 2 2 2 2
\/z +l +3 ., ,6 —z —l
2/ e—Jkfltan_1\/ dfi
\/22+l2 l

Z2 2 ,/ 2

+l +82 / 5 —z2 2e—jkfia2

\/ 2+12

/3 Vz2+12+s2 jkfi [3 l ,6 ( )
e" d dad . 2.9

0 v 02+z2+l2 02 +12

After defining a new variable of integration and subtracting the singularity at z = 0

 

 

leadfi

 

 

from each integrand, the resulting FNM expression for a right triangle (Figure 2.1a),

the sum of two right triangles (Figure 2.1b) that share a common side of length l, and

11

the diflerence between two right triangles (Figure 2.1c) that share a common side of

length l is

 

_. jLUt I l ' ~ 2 2 2 ' 7
P(z;k) — —pC—”—e—/ B (49’ka" +3 +1 —e_-7k”) do, (2.10)

_ 27;- :1: C 02 + [2
where l represents the height of the triangle, and :1: B and 2:0 represent the x-
coordinates of B and C, respectively. In Figures 2.1a, 2.1b, and 2.10, the values
of (23,320) are (3,0), (31,—32), and (31,32), respectively. Thus, a single FNM ex-
pression represents all three triangle geometries in Figure 2.1, whereas the impulse

response requires a separate expression for each triangle in Figure 2.1.

'IYansient FNM calculations

The inverse Fourier transform of Eq. 2.10 generates the F NM expression for the
transient response. The transient pressure generated by a triangular source above

the vertex A is represented by

zzt =—— —— v t—— 2: +0 +l —vt—z c do, 2.11
p( .) 27r $0 [2+02 C ( /) ( l

where the transient excitation is represented by v(t). By retaining the v(t —- z / 6) term

within the integral and subtracting the singularity, Eq. 2.11 maintains the rapid rate

of convergence achieved for time-harmonic calculations with Eq. 2.10.

Time-Space Decomposition

Transient FNM computations are accelerated by decoupling the temporal and spatial
dependence of Eq. 2.11. The time-space decomposition approach, demonstrated pre-
viously for a circular source [36], expands the delayed input pulse v(t - r) in terms of

temporal weighting functions 971(t) and spatially-dependent terms fn(r) that depend

12

Table 2.1. Basis functions for time-space decomposition with a Banning-weighted

 

 

pulse.

temporal basis functions 971(t) 1 spatial basis functions fn('r)

91m = %sin(27rf0t) m7) = cos(27rfOT)

g2(t) = —% cos(27rf0t) f2(7') = sin(27rfOT)
g3(t) = —% cos (2V?) sin(27rf0t) f3(T) = cos 2%; COS(27Tf0T)
g4(t) = $005 2g;- cos(27rf0t) f4(’r) = cos 2a;- sin(27rfOT)
g5(t) = %sin 2a;- sin(27rf0t) f5(T) = sin 2&7": COS(27Tf0T)
96(t) = ésin 2&5!- cos(27rf0t) f6(T) = srn 2V7?- sin(27rf07)

 

only on the coordinates of the observation point and the variable of integration 0

through 7' = %\/ 22 + 02 + [2. The decoupled input pulse is thus represented by

N
m — T) = rect (ET—V1); fn(T)gn(t), (2.12)

where the time duration of the pulse is indicated by the parameter W. The decom-
posed pulse in Eq. 2.12 is then inserted into Eq. 2.11, and then time-dependent terms
are factored out of the integral. The result consists of N edge wave terms specified
by

_ pc IB M) t-7
E-n— 27ft IC 02+I2T€Ct W (10' (2.13)

and a direct wave term given by

at
D = —Ev(t — z/c)l/ B ——2——1—-—2-d0. (2.14)
27f $0 (7 +1

The temporal dependence of the edge wave integrand in Eq. 2.13 is eliminated when
the effect of the rect function is instead shifted to the limits of integration. This oper-
ation, which restricts the edge wave contributions by only considering those that have
reached the observation point without completely passing the observation point, com-
pletely removes all temporal variables from the integrand. As a result, calculations of
transient pressure fields are converted into the numerical evaluation and subsequent

superposition of N spatial integrals that are weighted by analytical timedependent

13

terms. Further reduction in the computation time is achieved by storing redundant

edge wave calculations from Eq. 2.13 in the matrix

n(ij)= Z ,Umf_7_2___z(‘rl:m2l)_ (2.15)

In Eq. 2.15, tum represents the weights and am represents the abscissas computed

for Gauss quadrature, the value of T[0ml is obtained from the relation T[O’m] =

 

%(/22 + 072” + Z2, and the indices 2' and j indicate the shortest and longest times
that correspond to the limits of integration. The values in Kn(z', j) are initialized
within the computation procedure only for the points that are needed, and then the
time-space decomposition calculations superpose the numerically computed results
of the spatial integrals with analytical time-dependent weighting factors to achieve
a significant reduction in computation time for transient pressure calculations in the

nearfield region.

2.1.3 Superposition calculations with impulse response and

FN M expressions

At observation points away from the normal that passes through a vertex of the tri-
angular source, impulse response and F N M calculations project the observation point
onto the source plane and then superpose the contributions from two or three triangles
as in Figure 2.2. The contributions from three triangles are either added, as shown
in Figure 2.2a for an observation point within the lateral extent of the source triangle
AABC, or added and subtracted as demonstrated in Figure 2.2b for an observation
point outside of the lateral extent of the source triangle. Whether a contribution is
added or subtracted depends on the location of the projected observation point in
the source plane relative to each side of the triangular source.

The FNM admits some additional simplifications for nearfield calculations of pres-

sures generated by the triangle AABC' in Figure 2.2. If the three lines that are co—

14

 

 

 

 

 

 

 

(b)

Figure 2.2. Superposition operations that calculate nearfield pressures generated by

an equilateral triangular source, where each side is four wavelengths long. The vertex
D (indicated in bold) is the projection of the observation point onto the source plane,
which partitions the radiating source into three triangles with sides (ai,bz-,cz-). (a)
The field point is located inside of the equilateral triangular source. (b) The field
point is located outside of the equilateral triangular source.

15

incident with the three sides of the source triangle AABC are defined in the general
form E51: + Fig + G,- = 0, then the distances from the projected observation point to
each of the three sides are represented by lz- = lEix+Fiy+Gi|/ E? + F 22. Likewise,
the sign of each contribution is defined as S,- = (El-1: + Fiy +Gi)/|Ez-x + Fig + Gil for
coefficients Ei, F11 and G,- chosen such that 32' is positive within the lateral extent
of the source AABC. Furthermore, the lower and upper limits of integration are
defined as ((122 — b? - €22)/(2ci) and (az2 + c? — b?)/(2cz-), respectively. The resulting

nearfield pressure generated by AABC in Figure 2.2 is therefore represented by

_pwejwt

P(:v,y,2; k)= 27,

 

Z3: Eix+Fiy+Gi

i=1 ‘/Ez-2+Fz.2
2 2 2

a-+c——b- _--/———2 2 2 .
_1_2_CZ_T__L8 3k 0 +2 +lz'_e—sz

x 2 do 2.16
1.2-12.2 .2 +13 < >

Ci

 

Calculations with Eq. 5.4 compute the values of Ci, E2" F21 and C,- only once for each
edge of AABC, whereas the values of a,- and b,- are calculated once for each (x,y)
pair. Unlike the expressions for the impulse response that change depending on the
spatial coordinate, Eq. 5.4 is a general formula that computes the nearfield pressure

with a single expression that is valid at all points in space.

2.1.4 Transient input waveform

Evaluations of the impulse response and the FNM with time—space decomposition are

performed for the Hanning-weighted pulse specified by
1
v(t) = 2 [1 — cos(27rt/W)] sin(27rf0t)'rect(t/W), (2.17)

where rect(t) = 1 if t 6 [0,1] and rect(t) = 0 otherwise. In the simulations that
follow, the input is a Hanning-weighted pulse with a center frequency f0 = 2MHz

and a pulse duration W = 1.5113. Time-space decomposition performed on this pulse

16

with N = 6 yields the entries in Table 2.1, where the spatial edge wave integral
in Eq. 2.13 is evaluated once for each row entry applied to each edge of the source
triangle AABC in Figure 2.2, and then the results are weighted by the temporal basis

functions in Table 2.1.

2.1.5 Error calculations

For time-harmonic nearfield pressure calculations, the numerical error 17(23, y, z) is
defined as the normalized difference between the reference field and the computed

field according to

|P($vya Z) _' Pref(x1y7 Z)l
maxIPTef(a:,y, z)|

 

71(1“, 11, 2) = (2.18)

where Pref“, y, z) is the reference time-harmonic nearfield pressure. For transient
nearfield pressure calculations, the numerical error n(;r, y, 2) between the computed

transient field and the reference transient field is defined by

”PCT, 31, Z; t) ‘Preflx, y, 2; t)”
maxx,y,z ”Prefix, y, 2; 15)“

 

77(1) y) 2:) = , (2.19)

where H - M denotes the energy norm used with respect to time, and 127.8 f(:1:, y, z; t)
is the reference transient pressure field as a function of time. The maximum error is
defined as 77mm: = max/33y; 77(513, y, z), and this value is computed for both time-

harmonic and transient excitations.

2.2 Results

All simulation programs are written in C, then compiled and executed within a
Matlab—C language MEX interface. The simulations are performed on an eMachines

T3958 personal computer with a 2.93MHZ Celeron D processor. The operating sys-

17

I
[’1 I I
I 1,1,”,
’I,”// ’o ' .
. 1,, II, fife/0,6,
’1 ”1’0; ’0 "0’0,
/ II, II] I], ’I/ ’I, In’
Illa/000,, ’ll,”I/,”I/," 1,, ,,,
I ’I, II I/ II] II] II] II I/ r
I’ll]I’llll’ll’lfil/Il,’Ill/’07I’I;I”ZI”Z'"I"
’I ’I, I, I] I] I] I] l/ l/ I, I/
I, I, I], I, l] I] I], l/, ’1, 'I, I: '
9,, ’I,’ I,, 1,, I], I], I], /,l l], 1,] I101,
’ I” I], 1,, I,, I], , .
9,000,, 0,,//, II, , II, I a,
’I ’II ’II /I II ’II I; '
phi/I”,1,,”//,,//,,,//II’I, 1,,
‘ I
1,,

I

’II I ”

] I] I
”Ill/loll” [”1”

.

III,

I ’I

,o,I],I,,",,"/,,1,,III,IIII,IIII,III,,,/II,,II/I/,ll///,,///I//,/I’I/,,l’l

1 o:o:o,'o,;o,0,,7,,,I;/,,;,,,;;/,,;;II/,,,;;/x/,,”Ix/,7],
"4’1, ’0 1,, ”I1 0”,,

normalized pressure
.0
01

{7’00“
I I, I
oft/I, I, I]
"""':":"I/’I/ ///
41/90, I,,

  

0.4

 
    

0.2 ,
2 (units of aZ/x) 0 3 y (units of half length a)

Figure 2.3. Simulated time-harmonic pressure field in the x = 0 plane for an equilat-
eral triangular source with sides equal to 4 wavelengths. The reference field is gener-
ated by the impulse response method computed with 100,000-point Gauss quadrature.

tem on this computer is Fedora Core 3 Linux. All simulations are run sequentially

under similar operating conditions.

2.2.1 Time-harmonic nearfield pressure calculations
Reference pressure distribution

The reference pressure field is computed in Figure 2.3 for an equilateral triangular
source with sides equal to 4 wavelengths. In Figure 2.3, the acoustic field is evaluated
in the a: = 0 plane defined in Figure 2.2. The reference nearfield pressure distributions

in Figure 2.3 are obtained when the impulse response is calculated for all triangles

18

with 100,000-point Gauss quadrature. This pressure distribution is selected as the
reference because nearfield pressures computed with 100,000 abscissas produce nor-
malized errors that converge to 15 significant digits throughout the nearfield region,

which represents the smallest error achievable with double precision arithmetic.

FNM and impulse response calculations

The numerical errors and computation times for the fast nearfield method and the
impulse response method are shown in Figure 2.4. For the FNM, the exact impulse
response, and Field II with ‘use_triangles,’ nearfield pressures are evaluated in an 81 by
101 point grid in the :c = 0 plane as shown in Figure 2.3. Field II with ‘useJectangles’
and the smoothed impulse response require an offset due to a singularity on the piston
face and are therefore evaluated on a smaller 81 by 86 point grid. The FN M and the
exact impulse response are evaluated with Gauss quadrature, and all three integrals
corresponding to the three sides of the source triangles are evaluated with the same
number of abscissas. The remaining methods that approximate the uniformly sampled
impulse response (i.e., Field II and the smoothed impulse response) are evaluated
with the midpoint rule as described in the user’s guide on the Field II web site
(http://www.es.oersted.dtu.dk/staff/jaj/field/). Figure 2.4 shows that the error for
a given computation time is consistently smaller with the FNM, where smaller errors
are located nearer to the horizontal axis on the bottom of this log-log plot. Likewise,
the time required to achieve a given error is consistently smaller with the FNM, since
the FNM plot is consistently located to the left of the impulse response plot.
Comparisons between the impulse response and the FNM evaluated for the same
peak error are summarized in Table 2.2. For a 10% peak error, the FNM is 4.39 times
faster than the impulse response, and for a 1% peak error, the FNM is 3.44 times
faster than the impulse response for this grid and piston geometry. Even greater im-

provements are observed for smaller peak error values due to the rapid convergence

19

O

 

—L
O
1

 

 

 

 

 

 

 

x .

(U - - -- =."‘Ifl: -

E skid, L7 17::
g 9 - -o- - -o- - o - -o
5 5 ' °
t 10 0°" *
(D
.8 FNM

—o—
.L‘.‘
_(5 ---o-- impulse response
g 1 0‘1 O '-I-- Field ll (use_triangles) _
8 -A- Field ll (use_rectangles)
f5 - 0 - smoothed impulse response
03
Q.
-15 J -

 

1o 3

10‘1 10° 101 1o2 10
computation time (seconds)

Figure 2.4. Peak normalized error for calculations of nearfield pressures generated
by the triangular source in Figure 2.2 plotted as a function of the computation time.
The results show that the F NM consistently achieves smaller errors in less time than
exact and approximate impulse response calculations for time-harmonic excitations.

of the FNM. Although these values change somewhat for different source and grid ge—
ometries, the FNM is consistently faster than the exact and the approximate impulse

response for nearfield calculations of time—harmonic pressures.

Field II calculations

The Field II simulation program [10] includes the ‘use_triangles’ option for calcula-
tions that model rectangular and triangular pistons as the superposition of triangular
sources. For calculations of time-harmonic pressures with the ‘use-triangles’ option

applied to the source geometry in Figure 2.2, Field II requires a temporal sampling

20

frequency of f3 = 16MHz to achieve a peak error of 10%. The computation time
for Field II with ‘use.triangles’ is 29.02 times slower than the F NM evaluated on the
same grid. For time-harmonic calculations, Field 11 with the ‘use_triangles’ option
requires a temporal sampling frequency of f3 = 32MHz to achieve a peak error of
1%. This results in a computation time that is 26.43 times longer than that required
for the F N M evaluated on the same grid.

Field II also provides a ‘useJectangles’ option that introduces a numerical singu-
larity on the piston surface, so the pressure is evaluated on a smaller 81 by 86 point
spatial grid that is offset from the piston face. Field II with ‘useJectangles’ evaluated
on this reduced grid produces a 10% peak error in 0.8908 seconds, which is 9.5 times
slower than the FN M on the full 81 by 101 point grid. For a 1% peak error, Field II
with ‘useJectangles’ computes the result on the restricted grid in 218.2492 seconds,

which is 1697 times slower than the FNM on the full grid.

Smoothed impulse response calculations

Time-harmonic calculations with the smoothed impulse response [35] evaluate the
pressure on a smaller 81 by 86 point spatial grid that is offset from the piston face. The
offset is required for smoothed impulse response calculations so that the singularity in
Eq. 2.7 at the piston face is avoided. For calculations of the time-harmonic pressure
generated by the triangular source depicted in Figure 2.2 and evaluated within an 81
by 86 point subset of the grid shown in Figure 2.3, the smoothed impulse calculation
converges to a peak error of 10% with a temporal sampling rate of f3 = 32MHz.
This computation is completed in 1.0105 seconds, which is 10.74 times longer than
the time required for the corresponding FNM calculation evaluated on a larger 81
by 101 point spatial grid. Time-harmonic calculations with the smoothed impulse
response achieve a peak error of 1% for a temporal sampling rate of f3 = 128MHz.

This computation is completed in 3.9983 seconds, which is 31.26 times longer than

21

     
 
 
 
 
  
  
  
  

   
   
  
 

 
  
    
  

  
 
 

 

g
a 0.5. .. ; (all
a) ,v. 9,1"
'5- :i‘il‘ l
"5‘." 33‘ 5
'O 0 ' ‘11“ (l
I 3(111‘(),’,‘1?‘;::,(‘,“

a) J)L:;.’o,t:..¢:¢‘1\\l‘lo’,
.u 0...“ ‘““..$.“ ‘\‘\““:‘\‘: :’/'/ ’0’ 49’2”
— _ , . cit. “3.?“ .| “\““’."z::‘| $2,017 3,1912, ’17,’
m -0.5 \ ' ' ' i I . l : 1":7-3 , 3‘39. ’3‘f’3”'»,."’ off/1’,” ,” ,I,’ ..
h
8

'1 s -

1

2
2 (units of 320,) 0 3 y (units of half length a)

Figure 2.5. Simulated transient pressure field in the a: = 0 plane for an equilateral tri-
angular source with sides equal to 4 wavelengths. For this calculation, the excitation
is the Hanning-weighted pulse in Eq. 2.17, and the transient pressure is evaluated at
85 time points in an 81 by 101 point grid. The result is plotted at 1.8125/18 after the
initiation of the input pulse.

time required to obtain the FNM result with 1% peak error in Table 2.2.

2.2.2 'D'ansient nearfield pressure calculations
Reference pressure distribution

The reference nearfield pressure distribution for transient excitations is calculated
with impulse response waveforms that are sampled at f3 = 524.288GHz, zero padded,
and convolved with FFTs. The resulting temporal variations in the nearfield pres-

sure, which are evaluated for an equilateral triangular piston with 4 wavelengths on

22

 

I
—-L

.3
O
.4
.4
i

. . . . .F a 1 ,
x P “ ~‘A5 ‘ h
\ " .~ ~
m \ 3 T l ‘A‘ 4
E \ °~! i ‘ x
\ In,‘ \ x
- o -
10 2, s I. lb'u b-wu 7“
E \\ El. 0. "s “
l- ‘- ' '
: x J H" . l"s
- ~3‘0---o---o-4=g,-o--'la~
. ’ '''''' 0" +
-3 ., 'I

_x
O

 

L ,' o

E

. —t—time-space decomposition
~-o---impulse response

 

I
h

_x
O

.—u-- Field ll (use__triangles) 1
E - A- Field ll (use_rectangles) 3

peak normalized error 11

- 0 - smoothed impulse response

 

 

 

 

(J'l

 

10 . 1Q 10 10
computation time (seconds)

 

—-L
0I

_1‘ 3

_x
C

Figure 2.6. The peak normalized error plotted as a function of the computation time
for the F NM/ time-space decomposition method, the impulse response method, and
methods that approximate the impulse response. These errors and times are evaluated
for transient nearfield calculations of an equilateral triangular source with sides equal
to 4 wavelengths. The excitation for these calculations is a Hanning-weighted pulse
with a center frequency of 2MHz.

each side, are then downsampled and stored at f3 2: 16MHz. The reference field is
calculated for a sound speed of c = 1.5mm/us on an 81 by 101 point spatial grid
evaluated at 85 time points, and the result at time t = 1.8125113 is shown in Figure
2.5. This error reference is accurate to 5 significant digits for calculations in the a: = 0

plane.

23

FNM and impulse response calculations

Figure 2.6 shows the numerical error plotted as a function of the computation time
for the FNM with time-space decomposition and calculations based on the impulse
response method. All methods, except for Field II with ‘useJectangles’ and the
smoothed impulse response, are evaluated relative to an 81 by 101 spatial point by
85 time point reference transient pressure distribution. Field II with ‘useJectangles’
and the smoothed impulse response are singular at the piston face, so a smaller 81
by 86 point spatial grid that incorporates an offset from the piston face is again
required for transient field computations. The input for the reference is generated
by a Hanning—weighted pulse with a center frequency of f0 = 2MHz. The transient
nearfield pressures are compared for f5 = 16MHz, which is the original sampling
rate for the FNM calculations and the resulting rate after downsampling for impulse
response calculations. Figure 2.6 shows that the FNM with time-space decomposition
is consistently faster than the impulse response and the methods that approximate
the impulse response. Similarly, Figure 2.6 indicates that the FNM with time-space
decomposition achieves much smaller numerical errors than the impulse response and
approximations to the impulse response.

Table 2.3 shows that the FNM with time-space decomposition achieves a 10% peak
error with 5 Gauss abscissas in 0.4867 seconds. To achieve a 1% peak error, the FNM
with time-space decomposition needs 9 Gauss abscissas and the computation time
is 0.6160 seconds. In contrast, the impulse response method achieves a peak error
of 10% with a sampling frequency of f3 = 128MHz in 1.8911 seconds. To achieve
a peak error of 1%, the impulse response method requires a sampling frequency of
f3 = IGHz and a computation time of 23.5241 seconds. Thus, the reduction in the
computation time with time-space decomposition applied to the FNM relative to the
impulse response is a factor of 3.89 for a peak error of 10% and a factor of 38.19 for

a peak error of 1%.

24

Field II calculations

The Field 11 result obtained with the ‘use_triangles’ option for the transient excitation
in Eq. 2.17 requires a sampling frequency of f3 = 16MHz to achieve a peak error of
10%, and the computation time for this combination of parameters is 5.3317 seconds.
For a peak error of 1%, Field 11 with the ‘use_triangles’ option requires a sampling
frequency of f3 = 64MHz, and the computation time is 6.8078 seconds. Therefore,
the F NM with time-space decomposition is 10.95 times faster than Field 11 with
‘use_triangles’ for a peak error of 10% and 11.05 times faster for a peak error of 1%.

Transient Field 11 calculations that subdivide the aperture into small rectangular
sources with ‘useJectangles’ reach a peak error of 10% with a temporal sampling
frequency of f3 = 32MHz in 3.9926 seconds. Field II with ‘useJectangles’ achieves
a peak error of 1% with a temporal sampling frequency of f3 = 48MHz in 221.6569
seconds. Therefore, the FN M with time-space decomposition is 8.2 times faster than
Field II evaluated with subdivided rectangular sources for a 10% peak error and
359.81 times faster than Field II evaluated with subdivided rectangular sources for a

1% peak error.

Smoothed impulse response calculations

Transient calculations with the smoothed impulse response in Eq. 2.7 evaluate the
pressure at 85 time points on a smaller 81 by 86 point spatial grid that is offset from
the piston face. The offset is required in order to avoid the singularity in Eq. 2.7
on the piston face. For calculations of the time-harmonic pressure generated by the
triangular source depicted in Figure 2.2 and evaluated within an 81 by 86 point subset
of the grid shown in Figure 2.3, the smooth impulse calculation converges to a peak
error of 10% with a temporal sampling rate of f3 = 32MHz. This computation is
completed in 1.0122 seconds, which is 2.08 times longer than the time required for the

corresponding FNM calculation evaluated on a larger 81 by 101 point spatial grid.

25

Time-harmonic calculations with the smoothed impulse response achieve a peak error
of 1% for a temporal sampling rate of f3 = 128MHz. This computation is completed
in 4.1261 seconds, which is 6.7 times longer than time required to obtain the FN M

result with 1% peak error in Table 2.2.

2.3 Discussion

2.3.1 Time and error calculations

W'hile computer processor speed and memory has increased substantially in recent
decades, the size and complexity of ultrasound therapy and imaging simulations has
grown accordingly. Simulations of large ultrasound therapy arrays are now applied to
thousands of transducer elements and computational volumes spanning hundreds of
wavelengths in three dimensions, and simulations of diagnostic imaging arrays have
demonstrated a corresponding increase in the number of active elements and the
number of scatterers. As a result, large simulations of ultrasound phased arrays can
require 24 hours or longer on modern computers.

For these large simulations, the evaluation of computational time and numerical
error is essential. The computation time remains the primary bottleneck in these
time-consuming calculations, but fair comparisons of computation time also require
calculations of the numerical error. In recent years, evaluations of the numerical error
have been neglected due to the slow convergence of the impulse response and methods
that approximate the impulse response. Figures 2.4 and 2.6 demonstrate this slow
convergence, which is further emphasized by the time-harmonic reference field that
requires 100,000 Gauss abscissas for convergence to 15 significant digits and by the
transient reference field that requires a sampling frequency of f3 =524.288GHz for
convergence to 5 significant digits.

The rapid convergence of the FNM demonstrated in Figures 2.4 and 2.6 suggests

26

that the FNM is ideal for calculating nearfield pressure reference fields. In Figure 2.4,
time-harmonic FNM calculations converge within 15 significant digits in less than
one-third of the time that the impulse response requires for convergence within 5
significant digits. Likewise, in Figure 2.6, transient FNM calculations with time-space
decomposition converges within 5 significant digits in less than one-fifth of the time
that the impulse response requires for convergence within 2 significant digits. In these
simulations of a triangular piston source excited by a pulse with a center frequency
of 2MHz, impulse response calculations require a sampling rate of lGHz to achieve
only 2 significant digits of accuracy, whereas the FNM with time-space decomposition
requires only 9 Gauss abscissas applied to each integral and a sampling rate of 16MHz

to achieve 2 significant digits of accuracy throughout the nearfield region.

2.3.2 Advantages of the FNM for time-harmonic and tran-

sient calculations

The computational advantages of the FN M are obtained from several sources. First,
the FNM replaces time-consuming calculations of inverse trigonometric functions with
a ratio of polynomials in the integrand. This reduces the computation time without
increasing the numerical error. Second, the FNM reduces the numerical error by
subtracting a singularity in the integrand. This step, which reduces the numerical
error without significantly increasing the computation time, is particularly effective in
eliminating numerical problems that occur along the edge of the source and through-
out the paraxial region. Third, the FNM defines a single analytical expression that
describes the pressure throughout the nearfield region, whereas the impulse response
requires multiple expressions to define the field generated by a single source. Thus,
relative to calculations that employ exact or approximate calculations of the impulse
response, convergence is faster with the FNM, and the FNM expressions are easier to

evaluate.

27

The advantage of the F NM with time—space decomposition is that an integral
expression with temporal and spatial dependencies is replaced with an equivalent
expression that instead evaluates N spatial integrals for each edge of the triangular
source and weights the result of each integral with an analytical temporal term. This
results in greatly reduced overhead for transient nearfield calculations, considering
that the impulse response requires sampling rates of f3 = 128MHz for a peak error
of 10% and f3 = 1GHz for a peak error of 1% for the source geometry in Figure 2.2.
The F N M eliminates these high sampling rates, which therefore facilitates much more

efficient utilization of computer memory.

2.3.3 Field II

The Field II calculations with ‘use_triang1es’ are evaluated within the same 81 by
101 spatial grid defined previously for these nearfield calculations, whereas the same
calculations with subdivided rectangular subapertures (i.e., ‘use.rectangles’) are eval-
uated in an 81 by 86 spatial grid that includes an offset from the piston face. The
offset is required for these nearfield calculations, otherwise the error grows exces-
sively large on the piston face, which translates into much longer computation times
for 10% and 1% peak errors. This occurs because subdividing the aperture introduces
a numerical singularity on the piston face. Although Field II reduces the sampling
frequency relative to other impulse response calculations, the exact impulse response
consistently outperforms Field 11 for these time-harmonic nearfield calculations, and
the FN M evaluated with Gauss quadrature outperforms both of these by a wide mar-
gin. Furthermore, the FNM with time-space decomposition is also considerably faster
than Field II for transient nearfield calculations, and the FNM with time-space de-
composition, unlike Field II with ‘useJectangles,’ allows the computational grid to

extend up to the piston face.

28

2.3.4 Smoothed impulse response

Unlike the FNM and the exact impulse response, the smoothed impulse response
requires an offset from the piston face for nearfield calculations. This offset is re-
quired because the denominator in Eq. 2.7 produces a numerical singularity on the
piston face. Despite evaluating the nearfield pressure on a smaller 81 by 86 point
spatial grid, the smoothed impulse response is slower than the FNM and the exact
impulse response for time-harmonic calculations, as demonstrated in Figure 2.4 and
Table 2.2. The exact impulse response is faster than the smoothed impulse response
for these time-harmonic calculations because the exact impulse response is evaluated
with Gauss quadrature, and Gauss quadrature generally converges much faster than
other numerical integration methods that uniformly sample the integrand. For tran-
sient calculations with both exact and approximate impulse response expressions,
uniform sampling is required for convolutions with the FFT. In these transient calcu—
lations, the smoothed impulse response gains some advantage over the exact impulse
response by evaluating the pressure at a smaller number of spatial grid points and by
reducing the problems with aliasing at higher frequencies. Nevertheless, as demon-
strated in Figures 2.4 and 4.3, the smoothed impulse response converges more slowly

than the FNM for time-harmonic and transient nearfield calculations.

2.4 Conclusion

A fast nearfield method is presented for numerical calculations of the pressure gen-
erated by a triangular source. For time-harmonic nearfield computations, the FNM
expression in Eq. 2.10 achieves smaller peak errors in less time than the exact im—
pulse response, the smoothed impulse response, and the Field II program. The results
show that the F NM is 4.39 times faster than the exact impulse response for a 10%

peak error, and the FNM is 3.44 times faster than the exact impulse response for

29

a 1% peak error. The F NM is at least an order of magnitude faster than Field II
and the smoothed impulse response for time-harmonic calculations compared at 10%
and 1% peak error values. In transient nearfield computations, the F NM in Eq. 2.11
combined with time-space decomposition achieves a substantial reduction in the com-
putation time relative to exact and approximate impulse response calculations for a
given peak error value. Transient nearfield pressures are evaluated with a Hanning-
weighted broadband pulse, and the resulting transient calculation is transformed into
the superposition of six spatial integrals. The results demonstrate that the FN M with
time-space decomposition is 3.89 and 38.19 times faster than the impulse response for
peak errors of 10% and 1%, respectively, evaluated on an 81 by 101 spatial grid at 85
time points. Comparisons between smoothed impulse response results evaluated on
the smaller 81 by 86 point offset spatial grid and the FNM with time-space decompo-
sition evaluated on the larger 81 by 101 point spatial grid indicate that the FNM with
time-space decomposition is 2.08 times faster than the smoothed impulse response for
a 10% peak error and the FNM with time-space decomposition is 6.7 times faster for
a 1% peak error. Compared to the Field II program, the FNM is at least an order of
magnitude faster for 10% and 1% peak error values. The results also suggest that the
F NM, which eliminates the numerical problems that are encountered in exact and
approximate impulse response calculations, provides a superior reference for nearfield

pressure calculations evaluated with time—harmonic and transient inputs.

30

Table 2.2. Number of Gauss abscissas, computation times, and computation times
relative to the FNM that describe the reduction in the computation time achieved
with the fast nearfield method relative to the impulse response and methods that
approximate the impulse response for peak errors of 10% and 1%. The FNM and
exact impulse response results are evaluated for time-harmonic calculations on a 81
by 101 point grid located in the a: = 0 plane, and the Field II and smoothed impulse
response results are evaluated on an 81 by 86 point grid in the :1: = 0 plane that is
slightly offset from the transducer face. (a) For a 10% peak error and (b) for a 1%
peak error.
(a)

Time-Harmonic Nearfield Computations

 

 

10% peak error
impulse Field II smooth Field II

FNM response ’use_triangles’ imp. resp. ’useJectangles’

Parameters N =8 N=11 f3 = 16MHz f3 = 32 f3 = 32MHz

 

 

 

 

 

MHz N=16x 16
Time 0.09385 0.41125 2.72125 1.00515 0.8908
Computation 1x 4.39x 29.02x 10.72x 9.5x
time relative
to the FNM

 

 

 

 

 

 

 

(b)

Time—Harmonic Nearfield Computations

 

 

1% peak error
impulse Field II smooth Field II

FNM response ’use_triangles’ imp. resp. ’useJectangles’

Parameters N211 N212 f3 = 32MHz f3 = 128 f3 = 48MHz

 

 

 

 

 

MHz N=256x 256
Time 0.12865 0.4419s 3.39935 4.01305 218.24925
Computation 1x 3.44 x 26.43x 31.21 x 1697.00x
time relative
to the FNM

 

 

 

 

 

 

 

 

31

Table 2.3. Comparisons of computation times, input parameters, and computation
times relative to the F NM that describe the reduction in the computation time
achieved with the F NM and time-space decomposition relative to the exact and ap-
proximate impulse response for specified maximum errors of 10% and 1%. For FNM,
impulse response, Field II calculations with ‘use_triangles’, and Field II calculations
with ’useJectangles’, these transient results are evaluated in an 81 by 101 spatial
point by 85 time point grid, and for the smoothed impulse response, the results are
valued at the same temporal points in a restricted 81 by 86 point spatial grid. (a)
For a 10% peak error and (b) for a 1% peak error.

(a)

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Transient Nearfield Computations
10% peak error
impulse Field II smoothed Field II
FNM response ’use_triangles’ imp. resp. ’useJectangles’
Parameters N: f3 = 128 f3 = 16MHz f5 = 32 f3 = 32MHz
MHz MHz N=16 x 16
Time 0.48675 1.89115 5.33175 1.01225 3.99265
Computation 1x 3.89x 10.95x 2.08x 8.2x
time relative
to the FNM
0))
Transient Nearfield Computations
1% peak error
impulse Field II smoothed Field II
FNM response ’use_triangles’ imp. resp. ’useJectangles’
Parameters N =9 f3 = 1 f3 = 64MHz f5 = 128 f3 = 48MHz
GHz MHz N=256 x 256
Time 0.61605 23.52415 6.80785 4.12615 221.65695
Computation 1x 38.19x 11.05x 6.70x 359.81x
time relative
to the FNM

 

32

 

 

CHAPTER 3

A 2D Fast Nearfield Method for

Apodized Rectangular Pistons

Although several methods, including the Rayleigh-Sommerfeld integral [6] and the
Field II program [10, 11], calculate the pressures generated by apodized rectangular
pistons, the numerical performance of these methods suffers in the nearfield region.
The numerical evaluation of the Rayleigh-Sommerfeld integral converges slowly in
the nearfield region because of the singularity introduced by the 1 / R term, which
approaches infinity when R approaches zero. The Field II program subdivides a
rectangular piston into smaller rectangles and computes the pressure using the far
field approximation for the impulse response of the velocity potential, which also
contains a l/R term. Thus, the Field II program generates relatively large errors
and converges slowly in the nearfield region, especially near the piston face. To
address this problem for circular pistons, the FNM has been recently extended to
include axisymmetric particle velocity distributions, and the resulting 2D integral
also demonstrates rapid convergence [37]. However, this apodized F NM expression is
specific to pistons with circular or cylindrical symmetry, and methods for modeling

apodized rectangular pistons are still needed.

33

To improve the performance of nearfield calculations for apodized rectangular pis-
tons, an apodized FNM expression is derived from the FNM expression for uniformly
excited rectangular pistons. The derivation of the apodized F NM for rectangular pis-
tons begins by subdividing the piston into small uniformly excited subelements. The
total pressure is obtained by superposing by the pressure produced by all of the subele—
ments. After performing a summation variable exchange and integrating by parts,
the apodized FN M expression for rectangular pistons is obtained. Next, the apodized
FNM for transient pressure calculations is obtained by inverse Fourier transforming
the time-harmonic apodized FNM expression. The apodized FNM expression, the
Rayleigh-Sommerfeld integral, and the Field II program are then evaluated in the
nearfield of a square piston that extends 4 wavelengths in both directions. The re-
sults of time-harmonic and transient computations indicate that, when compared to
calculations performed with the Rayleigh-Sommerfeld integral and the Field II pro—
gram, 1) the apodized F NM achieves the smallest errors for a given amount of time,

and 2) the apodized FNM requires the least time to achieve a given error.

3.1 Existing calculation methods

3.1.1 The Rayleigh-Sommerfeld integral

The time-harmonic pressure generated by an apodized rectangular source is also com-

puted with the Rayleigh—Sommerfeld integral [6] via

:jw_P___v0€jwt a. b )3;ij
0 0

where w is the excitation frequency in radians per second, p is the density of
the medium, c is the speed of sound, v0 is the constant normal particle veloc-

ity evaluated on the surface of the rectangular source, k is the wavenumber, and

34

 

R = \/(2: — p)2 + (y — u)2 + .22 is the distance between the source point (12,12, 0)
and the observation point (3:, y, z). The transient pressure generated by an apodized
rectangular source with a temporal excitation component v(t) is given by the inverse

Fourier transform of Eq. (3.1), which yields

R/ C)

pRayleigh($ y,z t)— " “(:b/flflfl ———-—-d udu, (3.2)

where 1'2(t) is the time derivative of the input excitation pulse v(t).

3.1.2 The Field II program

The Field II program [10] is a software package that computes transient and steady-
state pressures generated by phased arrays and individual ultrasound transducers.
The Field II program with the userrectangles option divides each piston source into
small rectangular elements and applies the far field approximation of the spatial
impulse response to each small rectangular element [11], where Field II specifias the
apodization at the center point of each small rectangular element. The accuracy of
the Field II program is dependent on two factors, namely the sampling frequency and
the number of small rectangular elements. With an increase in the sampling frequency
or the number of small rectangular elements, Field II achieves smaller errors, but the

computation time increases accordingly.

3.2 Fast nearfield method for apodized rectangu-
lar pistons

In the derivation that follows, each observation point is denoted by (13,31,21), and each
source point is denoted by (,1, 11,0). Figure 3.1 shows the coordinate system used in

the derivation and subsequent evaluations. The rectangular source is located in the

35

Yil

 

2b

 

 

 

0.99112 /4/\ "

 

Figure 3.1. Orientation of the computational grid relative to the rectangular source.
The rectangular source, which has width a and height b, lies in the z = 0 plane.
The dashed lines define the extent of the computational grid in the 1‘ = a/ 2 plane.
The extent of the computational grid is 2b by 0.99112 /4)\ in the y and 2 directions,
respectively.

36

z = 0 plane, and the origin of the coordinate system 0 coincides with the lower left
corner of the rectangular source. The apodization function of the source is given by
f (p, V), which is zero outside of the rectangular source. The width of the rectangular
source is a and the height is b.

The FNM expression for a rectangular piston that is excited uniformly is given by
McGough [8]. A more general unapodized FNM expression for a rectangular piston
is obtained from the expression in Chen and McGough [12] for a triangular piston.
Here, the FNM expression for the uniformly excited rectangular piston is denoted
by p0[p, u](x, y, z; k), where the subscript ’0’ indicates the uniform excitation and p
and V represent the width and height of the rectangular source, respectively. The

nearfield pressure for the uniformly excited piston is given by

. m- _' /2 2 2 .
pcvoejwti/Z e ]k a +7: +hi_e-sz

h-
27T z 0'2 + h?

 

pol/1,140.11. 2; k) = do, (3.3)

1:17;,-
where or is an integration variable as defined in [8, 12]. The values of mi and ”2'1
which are functions of p and u, represent the upper limits and lower limits of the
integral, respectively. In Eq. (3.3), the values of (mi, 71,-) are (m1,n1) = (m2,n2) =
(,u—x, —:c) and (m3,~n3) = (7714,714) = (u —y, —y), and the values of h,- are h1 = g,
11.2 = u — y, h3 = :c, and I14 = ,u — :1: for i = 1,2, 3, and 4, respectively. According to

Eq. (3.3), there are two special cases, namely

pol0. 1/](Jc,y. z; k) = 0 and Polfl, 01(1“. 'y, z; k) = 0, (3-4)

where the subelement has zero width or height, respectively. These special cases are
utilized in the following derivation of the apodized FNM expression for a rectangular

source.

37

yl

 

-————--
—-——¢--i

 

 

 

 

(j — 1)AV """"""""""""""

 

 

 

 

p.
(0, 0) (i - 1)AM my 9?

Figure 3.2. The decomposition of an apodized rectangular source into smaller rectan-
gles, where each small rectangle is Au wide and AV high. The apodization function
f (n, V) is defined as constant over each small rectangle.

3.2.1 Steady-state Apodized FNM Expression

The pressure field pap0d(:c, y, z; k) is obtained by subdividing the rectangular source
into N x N small rectangles, where the subscript apod indicates that the pressure
is computed with the apodization function f (,u, V). One of these small rectangles is
depicted inside of the rectangular source in Figure 3.2. The values Au = a/N and
AV = b/N denote the width and the height of each small rectangle in the x and y
directions, respectively, and S [2 j] represents the rectangle at the ith column and
jth row of the subdivision. The four coordinates of the vertices of S [1, j] are given
by ((i-1)A#. (J' -1)A#), (mu, (3' -1)Au), (73AM, J'Au), and ((i-1)Au, J'Alu)
for the lower left, lower right, upper right, and upper left coordinates, respectively.
The expression p[i, j] = p[2’A]r, j AV] (1:, y, 2) represents the pressure produced by the
uniformly excited rectangular source having width iAu and height jAu, where the

lower left corner is located at (0,0), and the aperture function over the rectangle

38

S[i,j] is defined by q[i,j] = f(z'A,u, jAu). The total pressure pap0d(x,y,z;k) is

approximated by

pap0d($,y,2;k) 2 Z Z qli.J'l(p[i.J'l+p[i-1,J'- 1] -p[i,J'-ll-pli-1.J'l). (35)
i=1j=1

According to Eq. (3.4), p[l,0] = 0, and p[0, l] = 0, l = 1, , N. By utilizing these

restrictions and rearranging Eq. (3.5),

papodtca ya 2; k) 2

N- 1N— . .
qli+1 J'+ll-(1[i+1 J'l)-(q{i J+ll- (11ml)
;1 glpfi] j](( AVA/J, AVAp)

 

+2195 N]( q[’i+1,1:;]#- (Il'i,Nl)A#

 

+ j: 11p[N,j](—Q[ N’j +1] Q[N’j])Au+p[N,N]q[N,N] (3.6)

All

is obtained. Letting N —> 00 such that Au —> 0 and AV —* 0, Eq. (3.6) becomes

papodmv y’ 2:; k) = pap0d1($ay1 Z; k) + pap0d2(xa 3’: Z; k)

8
where papodfi :r 31,2 k) =fa fb iffy—>155. Vl(1‘.y, 2; k)dltdv,

a3
pap0d2($7yizik) : — f0 P—féZ—ZPOl/J’abl(xayrzrk)dfly

3 .
pap0d3(x.y,2;k) = - 6’ Jégfllpolawl($.y,2;k)dm
and pap0d4(:r,y,z;k) = f(a,b)p0[a,b](:r,y,z;k). Within these expressions,

39

p0[p,u](:1:,y,z;k) is the single integral from Eq. (3.3), so pap0d1(:r,y,z;k) is ac-
tually a triple integral, pap0d2(:r,y, z; k) and pap0d3(:r,y,z; k) are double integrals
that admit further simplification, and pap0d4(:c, y, z; k) is already simplified.

After substituting the uniformly excited fast nearfield method expression

POlI‘aVKxai/a Z; ’9) from Eq. (3.3), pap0d1(:r,y, z; k) in Eq. (3.7) becomes

jwtz 62(
, , _p__2__cvo: f (1W!)
.W) :43] f /—

 

 

i=1nZ-O
e—jk,/02+z2+h2_ e—jkz
x h. dadpdu. (3.8)
2 02 + I122
The derivation of the first two triple integrals
u—x a b
.k _pC’Uo:jwt 82““
pap0d1,i:1’2(xvyizi ): T—aau)
—:r 0 0
_' / 2 2 2 .
h e ]k 0' +2 +hi *e__szd d d (39)
x - a ,u u .
z 02 + h?
is outlined here, where pap0d1(x,y,z;k) = pap0d1,2-:1(:r,y,z;k) +

papod1,z'=2($? y, Z; k) + 19111901115239?)2 .z; k) + 19111902112243?!) Z; k)- Let 9i=1,2(#)

denote the integral

— 2
( ) p/x, e—Jk,/02+22+hz. —e‘jkzd (3 0)
'2 u = 1' a. '1
g, 1’2 z 02+h22

—$

 

where h1 = y and h2 = V - y as indicated above. The function 9,21,2(u) is defined
in terms of a variable in the upper limit of the integral, so the derivative of 9,21,2(11)

with respect to ,u is

 

_-.. _ 2 2 2 .
I e J’W/(l‘ 33) +z +hz' _e—jkz
9i:1,2(l1) : hi

 

(,1 — 2:)2 + 1.22 (3'11)

40

After integrating by parts with respect to the variable y, papodl i=1 2(rr,y, z; k) is

rewritten as

. b
pcvoejwt / [(3701.11)

pap0d1,i:1,2(xayv Z; k) = _ 27F all 9121,2(#)I8
0
“an > ,
um
—/—8—V——gi:1’2(p)dfl du. (3.12)
0

I
According to Eq. (3.10), 9,21,2(0) = 0. Substituting the expression for 91:1 2(a)
into Eq. (3.12) and performing an exchange of variables yields the following analyti—

cally equivalent expression for the first triple integral,

 

b
, Jwt a a
_ ___P____Ct0€ f(aV)_f_(__u,1/)
pap0d1,i=1,2(xiy’z’k): O/a/ (—;— 69—1/—
0
2
he—jk)/(fl_) —:1: )2+z2+hi —e_]kz

 

dydu. (3.13)

 

. u—y a b
- ( 'k) __ pcheJWt //32f(/1,V)
pap0d1,z=3,4 1r,y,~, _ 27,. @1611
—y 0 0
—jlc‘/02+zQ+hz2 _ 63—sz
xh- dad/Adz! 3-14
z 02 + h2 ( )

in Eq. (3.8) is outlined in the following. Let 9,23,4(11) denote the integral

u— _- /2 2 2 .

9': V = h-
2 3,4( ) z 02+h2
_y z

 

do, (3.15)

41

where h3 = a: and h4 = p — :1: as indicated above. The function giz3,4(V) is defined
in terms of a variable in the upper limit of the integral, so the derivative of gi:3A(V)

with respect to V is

 

,- 2
I ( ) h e—]k\/(V—y)2+z2+hz. —e—jkz
51-: V = -
z 3’4 z (V-y)2+hz2

After integrating by parts with respect to the variable V, papod1 2:3 4(x,y, z; k) is

 

(3.16)

rewritten as

- a
pcvoejwt / [Bfm V)

 

. . _ b
pap0d1,i:3,4(1,y,z,k) — -‘ 2” a“ 9i=3,4(V)l0
0
b0f< ) ’
IM/
_/ Bu gi=172(l/)d1/ d,u. (3.17)
0

I
According to Eq. (3.10), 9,23,4(0) = 0. Substituting the expression for 92-:3 4(V)
into Eq. (3.12) and performing an exchange of variables yields the following analyti-

cally equivalent expression for the first triple integral,

 

b
t a
_, _p___cv0e9'w an (rub) _Bfmm)
pap0d1,i=3,4(x’y’zik): //(5,u (9p
0 0

,- . 2
e‘Jk\/(V“y)2+z2+hi _ 8.3-]...
(V - y)2 + h?

By substituting Pol/1: b](:r:, y, z; k) and p0[a, V] (:13, y, z; k) from Eq. (3.3) into

 

 

x h,- dudV. (3.18)

papod2($7y’ z; k) and pap0d3(x,y,z; k) and by using integration by parts,
papod2(x’ y, z; k) and papod3($’ y, z; k) are converted into the sum of two single in-
tegrals and two double integrals, while pap0d4 (:12, y, z; k) is already a single integral.
After the four terms in Eq. (3.7) are added and common terms are canceled, the

complete F NM expression for an apodized rectangular piston is:

42

 

 

 

 

a b — k +z2+hfi -.
J _
-k — pet—Li“; T h \/"12 e Md d
papodcrayaz) )‘—_ 12 1282 ’12,“. V
10 0 “12‘+ 12'
_ 2 .
_p_2__w0:jwt4 ejkf22+z +h2’b- e_]kz
Z [Tm-122,6 2 do (3.19)
. h2 .
izlnz- 022+ 22

where the values of Tli are T11 = T12 = _fléfigfll and T13 = T14 = —6f(i:’u ,
the values of ah- are 0‘11 = 0112 = p -— a: and a13 = 014 = V — y, and the values
of flu are hll = y, h12 = V - y, [113 = :12, and flu = u — :r. The values of T2,- are
T21 = T22 = f(a, b) and T23 = T24 = f(a,a), the values of 022- are 0121 = (122 =
a - :1: and 023 = (124 = a —— y, the values of h2z‘ are h21 = y, h22 = b - y, h23 2: 2:,
and h24 = a — :L', and the values of (mimi) are (m1,n1) = (m2,n2) = (a, 0) and
(m3,n3) = (m4,n4) = (b, O) for 2' = 1,2, 3, 4.

The apodized F NM expression in Eq. (3.19), which contains the summation of
four double integrals and four single integrals, describes the pressure generated by an
apodized piston for any boundary condition. The apodized FNM expression in Eq.

(3.19) admits further simplification if the apodization function is equal to zero on the

piston edge, where

f(0, V) = O, f(a,V) = 0, f(p,0) = 0, and f(p, b) = O. (3.20)

The boundary conditions given by Eq. (3.20) are equivalent to setting all of the terms
ng- equal to 0, so the single integrals with respect to a in Eq. (3.19) disappear. The
resulting apodized F NM expression is the summation of the four double integrals with
respect to p and V, so only the first line of Eq. (3.19) is needed when the boundary
values are all zero.)

The FNM expression for the uniformly excited rectangular piston is also a special

case of the apodized FNM expression in Eq. (3.19). For the uniform case, the

43

apodization function is f (,u,V) = rect(p/a)rect(V/b) where rect(t) = 1 if t E (0, 1)
and rect(t) = 0 otherwise. The weak derivative of f 01, V) is given by 8 f (p, V) / 8;» =
601) — 6(p — a) and 8f(,u, V)/8V = 6(V) — 6(V -— b). When this apodization function
is substituted into Eq. (3.19), all terms T2,- are equal to O, and the double integrals

reduce to a single integral, which is the same as po [(1, b] (:r, y, z; k) in Eq. (3.3).

3.2.2 Transient Apodized FNM Expression

The transient response for an apodized rectangular piston is obtained from the inverse
Fourier transform [37] of Eq. (3.19). Defining v(t) as the temporal component of the
transient normal particle velocity v(t) f (V, V), the transient pressure generated by an

apodized rectangular piston is represented by

 

 

 

 

a b
t _ pc 4 T -h .v(t— a¥i+z2+h%i/c)-v(t—z/c)d d
papodcrvyiz’ )— —'—7r 2 12 12 2 +h2 ,U. V
i=10 0 0‘12“ 12'
mi 2 2 2

_E Z ['13-th“ — \/a22. +2 + 1%./c) —v(t — z/c)d0 (3 21)

Z Z 2 2 '

2“ 1:17;, “22' + ’2'

where the values of T123 012': h”, T2,, 022': and h22- are listed immediately after Eq.
(3.19) in the previous section. The expressions for the Rayleigh-Sommerfeld integral
in Eqs. (3.1) and (3.2) are analytically equivalent to the apodized FNM expressions
in Eqs. (3.19) and (3.21), but the numerical properties of the two methods differ as

demonstrated in the results shown below.

3.2.3 Apodization function

The apodization function selected for comparisons between the apodized F NM, the
Rayleigh-Sommerfeld integral, and the Field II program is the product of sinusoidal

functions given by

44

 

 

J

 

  
     
     
   
   
     

    
 
  
 
 

 
   
  
   
    

 

    
    

             
        

   
   
 

> o 5 r I, I, 6'
. 1 a \
.1, 1, 1.
, u 1 u
- 11,,"1, '1/1 1 “mm.
w, 11,, q, I,,'.,' “mmum
v' 1,, 11,011, 1,, I, o, , “uuuuu
V .rm 11,, I, 1;, I, 1, I \u “mum.
g... .. 0,001,, 1,, II, 1, 1, v 1 “m “nun“ u ,
1 ,, 11,, 1,, 1,, I, 1, I, “N “mum“.
n, 11, ’I, 1, 1, I, \ m‘ “\unuu‘
n, '1, 1, I, 1, 1, I, \‘ “t“ \\\\\\\\u ‘
u ’1 ’I I, 1 I I \‘ \\\‘ nun
1 1 1 I 1 1 v \\\I
11,,”1,,’1,,,1,’I,,'1,'1 , “u“
, I I
1,011,004,110” '
1],]; 1 ”[1,/G
I c

   
  

 

1 I
11 ”1,", I

  
 

15°

    
 
   
   
 

  
 

\ \\
‘\\\ “\u “m u‘
\“\\\‘“\\\.\
u“ \\\““\‘\““\‘l““‘
\\ \\‘ “N
‘u

 

\

\\ “l “\
\‘ ‘\ \\

‘ “‘\\“‘\\““\

\

‘ \

 
 
 

\
‘\

|‘ \‘

“\“‘\““I \
‘ N u‘

“ ‘u“‘u\“‘\\
“.“‘\\““

\\

    

     

v (units of A) o 0 ll (units of 1)

Figure 3.3. The apodization function f (n, V) = sin(p7r/a) sin(V7r/b) evaluated on the
face of a 4A by 4A square piston. The maximum value of the apodization function is

achieved when p = 2A and V = 2A.

f(u, V) = sin(mr/a) sin(V7r/b). (3.22)

This function corresponds to the lowest order vibration mode of a rectangular mem-
brane with fixed edges [14, 38]. Eq. (3.22) is plotted for an apodized rectangular
source with each side equal to 4 wavelengths in Figure 3.3.

The apodized FNM equations of Eq. (3.19) and Eq. (3.21) admit additional sim-
plification when applied to the apodization function in Eq. (3.22). The apodization
function is the product of two sinusoidal functions, so the apodization function is
separable with respect to the variables y and V. After the derivative of the apodiza-

tion function with respect to p and V is substituted into Eqs (3.19) and (3.21), the

45

first and third double integrals contain [(1)3 cos(V7r/b)dV and f6 cos(,uTr/a)dp terms,
respectively. These two integrals are exactly equal to 0, so the first and third double
integrals in Eq. (3.19) and Eq. (3.21) are equal to 0 for the apodization function
in Eq. (3.22) at any observation point. Thus, only the second and fourth double
integrals are needed. Furthermore, in Eq. (3.19) and Eq. (3.21), the second and
the fourth double integrals share several terms. Shared terms in the apodized FNM

equations are always computed once and stored for use in repeated calculations.

3.2.4 Input transient pulse

For transient calculations, the excitation pulse v(t) is specified by the Hanning-

weighted pulse

v(t) = a1 — cos(27rt/W)] sin(27rft)rect(t/W), (3.23)

where rect(t) = 1 if t E (0, 1) and rect(t) = 0 otherwise. In the transient simulations
that follow, the center frequency f0 and the pulse duration W are f0 = 1.5MHz and

W = 2.0115, respectively, for the Hanning-weighted pulse.

3.2.5 Time space decomposition

For transient calculations, most of the computation time is expended while evaluating
v(t - T) in Eq. (3.21) and 1)(t — T) in Eq. (3.2). These terms are calculated at each
time t. The variable T is a function of the observation coordinates only, so T represents
the contribution from the spatial variable. The function v(t — T) can be separated
according to the time space decomposition approach in Kelly and McGough [37]. This

expression is decoupled as

M
v(t — T) = rect (t L—VT) fm(T)gm(t), (3.24)

m=1

 

46

where M =6 for transient F NM calculations with the Hanning—weighted pulse in Eq.
(3.23), W is the length of the pulse, and the fm(T) and gm(t) terms are given in
Table 3.1. The transient input to Eq. (3.2) is decoupled in the same manner, and
the corresponding fm(T) and gm(t) terms for i2(t — T) are given in Table 3.2 where
M = 10 for transient Rayleigh-Sommerfeld calculations with the Hanning-weighted
pulse. A simplified version of the time space decomposition algorithm is outlined
below.

1. Pre-compute and store gm(t) in advance for all values of t.

2. Evaluate T once for each spatial coordinate.

3. Compute the values of each fm (T) term immediately after T is calculated.

4. Calculate the value of v(t — T) according to Eq. (3.24).
By exploiting repeated calculations, this approach dramatically reduces the compu-
tation time without increasing the numerical error.

Table 3.]. Terms that define the time—space decomposition of the Banning-weighted
pulse v(t — T) for transient apodized FN M calculations.

temporal basis functions gm(t)

spatial basis functions fm(T)

 

 

910:) = %sin(27rf0t)

92(t) = —% cos(27rf0t)
g3(t) = -% cos ($7?) sin(27rf0t)
g4(t) = %cos ($69 cos(27rf0t)
g5(t) = —% sin (2%?) sin(27rf0t)
96(t) = ésin (21551:) cos(21rf0t)

3.2.6 Error Calculations

 

f1 (T) = cos(21rf0T)
f2(T) = sin(27rfOT)

f3(T) = cos 2a;- cos(27rf0T)
f4(T) = cos 2a;- sin(27rf0T)
f5(T) = sin 21377- cos(27rf0T)
f6(T) = sin EZW- sin(27rf0T)

Two error metrics are used in this paper. One is the normalized error distribution
n(:r, y, z; k), which describes the absolute value of the pressure difference at each spa.

tial point for time—harmonic calculations. The other is the normalized root mean

47

Table 3.2. Terms that define the time-space decomposition of the derivative of a
Hanning—weighted pulse iv(t —T) for transient calculations with the apodized Rayleigh-
Sommerfeld integral.

 

 

 

temporal basis functions gm(t) spatial basis functions fm(T)
91 (t) = 7Tf0 cos(27rf0t) f1 (T) = cos(27rf0T)
g2(t) = Trfo sin(27rf0t) f2 (T) = sin(27rf0T)
g3(t) = —7rf0 cos $6!- cos(27rf0t) f3(T) = cos 2131f; cos(27rf0T)
g4(t) = —1rf0 cos 2V? sin(27rf0t) f4(T) = cos 2VWVT' sin(27rf0T)
95(t) = —7rf0 sin 2&9!- cos(27rf0t) f5(T) = sin 26f;- cos(27rf0T)
96(t) = —7rf0 sin ‘21? sin(27rf0t) f6(T) = sin 217i;- sin(27rf0T)

V

g7(t = {£7 sin (%7{7£) sin(27rf0t) f7(T) 2: cos vafi: cos(27rf0T)
98(t) = -1711} sin (2V7?) cos(27rf0t) f8(7') = cos 2‘5?- sin(27rf0T)
27rT

99(t) = —% cos (2??) sin(27rf0t) f9(T) = sin W) cos(27rf0T)
910(t) = {£7 cos (Zr/76;) cos(27rf0t) f10(T) = sin (2V7?) sin(27rf0T)

 

square error (NRMSE), which describes the overall error performance with a sin-
gle value. For computed pressure field p(z,y, z; k) and the reference pressure field
prefix, y, z; k), the normalized error distribution n(:r, y, z; k) for each spatial point in

time-harmonic calculations is given by

77(16, :1. z; k) = |P(:r, y, z; k) - Prefixn ,z; k)|/lpref(rr, y, z; k)|max(x,y,z)- (3-25)

The value of n(:r, y, z; k) is shown in each depiction of the error mesh. The normal-
ized root mean square error (NRMSE) across all spatial points for time-harmonic

calculations is

 

 

NRMSE = Z Ip(1‘,y,2,k)—pref($ay,zik)l2/ Z lpref(z7yizlk)|27

x,y,z 1:,y,z

(3.26)

48

where 21mg, 2 denotes summation over all of the spatial grid points. The N RMSE

for transient calculations is

 

 

NRMSE: Z 1111.31.21;t>-p,.ef(x.y,z;t)l2/ Z Ipref(:c,y,z;t)|2,
:c,y,z,t $,y,2,t

(3.27)
where me’zi denotes summation over all spatial and temporal grid points. The val-
ues of the NRMSE are tabulated for each method evaluated in the results section 3.3.
The reference fields prefer, y, z; k) and pref”, y, z; t) for time-harmonic calculations
and transient calculations, respectively, are calculated using the Rayleigh-Sommerfeld
integral.

For nearfield pressure calculations with the apodized FNM expressions and the
Rayleigh-Sommerfeld integral, the number of abscissas in the p and V directions are
the same for a square source. The NRMSE values for both methods are computed as
the number of abscissas ranges from 2 to 100. The number of abscissas for a given
NRMSE is determined by the smallest number of abscissas that has an NRMSE
smaller than the desired NRMSE value. For nearfield pressure calculations with the
Field 11 program, the sampling frequency f3 and the number of small rectangular
elements are the two factors that determine the value of the N RMSE. In these calcu-
lations, the sampling frequency is varied from 16MHz to 160MHz with a step size of
16MHz, and the number of small rectangular elements in each direction ranges from
10 to 60 with a step size of 5. The NRMSE for each combination of the sampling
frequency and the number of small rectangular elements is then computed. The sam-
pling frequency and the number of small rectangular elements for a given N RMSE

are determined by the combination that has the smallest computation time.

49

3.3 Results

3.3.1 Time-harmonic pressure calculations
Reference pressure field

The reference pressure field in Figure 3.4 is computed for the apodization function
shown in Figure 3.3. In Figure 3.4, the acoustic field is evaluated in the a: = 2.0A
plane, where the grid in the :1: = 2.0A plane extends from —2/\ to 6A in the y direction,
and the grid spacing in the y-direction is 0.1A. The grid extends from 0.01a2/4A to
1.0a2/4A in the z direction with a spacing of 0.01a2/4A. The computational grid
in the z direction is shifted by 0.0102/4A relative to the source in the z = 0 plane
for all three methods in order to avoid the most severe singularities in this location.
Although the apodized FNM eliminates the worst singularities in the z = 0 plane, the
grid is nevertheless shifted slightly to reduce the problems that the Rayleigh integral
and the Field II program encounter on the piston face. The reference pressure field
shown in Figure 3.4 is obtained using the Rayleigh-Sommerfeld integral evaluated
with 100,000 Gauss abscissas in each direction, and the results are computed on an

81 x 100 spatial grid.

Error distributions

Figures 3.5 and 3.6 shows the normalized error distribution 17(x, y, z; k) between the
reference pressure field in Figure 3.4 and the simulated pressure field computed using
the apodized FNM expression, the apodized Rayleigh-Sommerfeld integral, and the
Field II program with the apodization function specified by Eq. 3.22. The region
near the piston face contains the largest errors for all three methods. The NRMSE
for the apodized F NM evaluated with N = 16 Gauss abscissas in each direction is
0.0005. The apodized Rayleigh-Sommerfeld integral evaluated with N = 16 Gauss

abscissas in each direction has an N RMSE of 0.0450. The Field II program computed

50

   
    
     
  
 
     
 

I ~~""“i§i\u<m
(8?“ \{lllliliihuuuuxu

   

 

   
 
   
   
  
  

      
   

  

  
    
  
     
   

  

 

    
   

        
   

 

 

   

 

          
 
 

 

  
 
 
   
   

    

          

 

”(536183 \\\\\\\\\\\\\\
g e: ‘willlllll\\\\\\\‘§‘lllli\\\\\\\\\\\\\\\
g 0.8\ . _ -.
' - \\ mu \\
o "‘\\\“““‘“~~.. \\\\\\\‘\ll\\\\\\\\\\\\ \\\\\\\\\\\ \\
b [,111.0‘\\\‘ \\\\\\\\\\8\\\\\\ \\\ \\\\“\\ \\\ \\\\\\\\ \
a 0.6“ I ml’;”lallmlmm‘\llllalilililll\\\\\\i“\\\ \\\\\\\\\\\\\\\\\\ll\\\\\\\\\\\\\\\\ \
ID Mfg/77333;}: “ll“
“’ III I x
N 0.4\ . I?”
= [I \ r
g 0 2 . Ill/Illllll/Illlllllilll1\\ \\\ “\\\ \\\“ “ \\““\x\“ M .I \.
.. . . z lull/”1155555555;at”.‘ ‘“‘“\“\\\“\““““‘““‘lllllll“‘“ll“ill““lllli\‘llllllttllist33333§§r=§§tr _
2 o ,lllltélmar» \\\\\\\\\\\\llllll\llil\lll\\xig}(llllllxgiltllgatlliilllliatllliatlliitiiiittggggg..wfi\t‘
’ - \\\\\\\ \\\\ \\\ x\\\xiit“rainiirsaggu3g . .
6 \\\\\\“\\\\\\\\\\l\\\\“xu‘xlll:gllS ‘
\\\\\“‘ \\\\“\\\““tttl\:::t: “
4 \\lllllllrtttW 1
2 " o 6 0.8
0 o 2 0.4
y (units of 7t) '2 0 ' ' 2
z (unlts of a I49»)

Figure 3.4. Simulated reference pressure field generated by an apodized rectangular

source with each side equal to 4 wavelengths. The results are evaluated in the a: = 2.0A
plane for a time-harmonic excitation.

with a 48MHz sampling frequency and 30 subdivisions in each direction achieves an

NRMSE of 0.0808.

Time vs. error comparisons for a timebhmmonic input

Figure 3.7 shows the computation times and the error values for nearfield pressure
calculations with a time-harmonic input applied to the apodized FNM, the Rayleigh-
Sommerfeld integral, and the Field II program. In Figure 3.7, the apodized FNM
achieves the smallast errors for a given computation time, and the apodized FNM uses

the least amount of time for a given error value. Comparisons between the apodized

51

normalized error 1]

 

  
 
 
 

0.8

   

0.6

 

0.4

. _ 0.2
Y (“""5 °f 7*) 2 o 2 (units of a2I4A.)

a) Apodized FN M.

P
'5‘

  
 
  
 

normalized error 11
P
5"

 

‘1

05°

0.8

     

o 2 0.4 0'6
y (units Of A) -2 0 . 2 (units of azl4x)

b) Rayleigh integral.

Figure 3.5. The normalized error distribution 1)(:1:,y, z; k) describes the difference
between the reference pressure field and the computed pressure field for an apodized
4A by 4) source. The error distribution 7) is plotted for a) the apodized FN M evaluated
with 16-point Gauss quadrature in each direction, b) the apodized Rayleigh integral
evaluated with 16—point Gauss quadrature in each direction.

52

A
I

normalized error 11
o
01
I

 

O
L

  
  

........ 1

z i V A 0.8
' 0.4 0.6

     

0.2

. _2
y (units Of f“) 0 2 (units of azl47t)

Figure 3.6. The normalized error distribution r)(z,y,z; It) describes the difference
between the reference pressure field and the computed pressure field for an apodized
4A by 4). source. The error distribution 17 is plotted for the Field II program evaluated
with f3 = 48MHz and 30 subdivisions in each direction.

F N M, the Rayleigh integral, and the Field II program are summarized in Table 3.3. To
achieve an NRMSE of 0.1, the apodized FNM needs 6 abscissas in each direction, and
the computation time is 0.07 seconds. The apodized FNM achieves an NRMSE of 0.01
in 0.18 seconds with 10 abscissas in each direction. Likewise, to achieve an NRMSE
of 0.1, the Rayleigh-Sommerfeld integral needs 12 abscissas in each direction, and
the computation time is 0.29 seconds. The Rayleigh-Sommerfeld integral achieves
an NRMSE of 0.01 in 2.25 seconds with 34 abscissas in each direction. Instead
of evaluating an integral with a certain number of abscissas, the Field II program

subdivides the aperture into small rectangular sources and calculates the impulse

53

 

10

-2

1o -----

NRMSE

    

 

--Apodized FNM g g
1o" '-°-'Ray|ei9h
'AIFie-d." . 3 55225:;
0.01 0.1 1 10 100
computation time (seconds)

I

 

 

 

 

 

 

Figure 3.7. Normalized root mean square error (NRMSE) plotted as a function of
the computation time for time-harmonic calculations with the apodized F NM, the
apodized Rayleigli-Sonunerfeld integral, and the Field 11 program. This figure demon-
strates that the apodized FNM achieves the smallest errors for a given computation
time, and the apodized FNM requires the smallest amount of time to achieve a given
error value.

54

response at a specified sampling frequency. To achieve an NRMSE of 0.1, the Field
II program requires a temporal sampling frequency of 64MHz and 20 subdivisions
in each direction, and the result is obtained in 4.61 seconds. The Field II program
achieves an NRMSE of 0.01 with a temporal sampling frequency of 160MHz and
50 subdivisions in each direction. This calculation is completed in 27.91 seconds.
Thus, for an NRMSE of 0.1, the apodized FNM is 4.14 times faster than the Rayleigh
integral, and the apodized FNM is 59.43 times faster than the Field II program. For an
NRMSE of 0.01, the apodized FNM is 12.50 times faster than the Rayleigh integral,
and the apodized FNM is 155.06 times faster than the Field II program. These
results are specific to the grid and piston geometry in Figure 3.3 and the apodization
function in Eq. (3.22). Although the results vary somewhat for different parameter
combinations, the apodized FNM consistently achieves the best performance in time
vs. error comparisons. The apodized FNM is also an ideal reference for time-harmonic

calculations due to the rapid convergence of Eq. (3.19).

3.3.2 Transient field calculations
Reference pressure field

The reference transient field is computed using the transient Rayleigh-Sommerfeld
integral in Eq. (3.2), where the numerical integration is performed with 100,000
Gauss abscissas applied in both the height and width directions. The input pulse is
the Hanning—weighted pulse in Eq. (3.23) with center frequency f0 = 1.5MHz and
pulse length W = 2.0113. The speed of sound is c = 1.5mm / us, and the wavelength at
the center frequency is /\ = 1.0mm. The sampling frequency is 16MHz. The transient
field is evaluated in the :1: = 2.0A plane, where the computational grid extends from
—2)\ to 6A in the y direction and from 0.01a2/4/\ to 1.0a2/4A in the z direction as
shown in Figure 3.1. The reference transient field is calculated on an 81 x 100 spatial

grid evaluated at 192 time points. Figure 3.8 shows two transient plots of the reference

55

Table 3.3. Simulation parameters for time—harmonic calculations that achieve nor-
malized root mean square error (N RMSE) values of 0.1 and 0.01. Parameters listed
include the number of Gauss abscissas or the corresponding Field II parameters, the
resulting computation time, and computation time relative to the apodized F NM for

the Rayleigh integral and the Field II program. (a) For a 0.1 NRMSE and (b) for a
0.01 NRMSE.

(a)

Time-Harmonic Calculations

 

 

 

 

 

 

 

0.1 NRMSE
Apodized FNM Rayleigh Field II
Simulation 6 x 6 12 x 12 f3 = 64MHz
Parameters abscissas abscissas 20 x 20 subdivisions
Time 0.073 0.293 4.613

Computation Time 1x 4.14 x 59.43 x

Relative to
Apodized FN M

 

 

 

(b)

Time-Harmonic Calculations

 

 

 

 

 

 

 

 

 

0.01 NRMSE
Apodized FNM Rayleigh Field II
Simulation 10x10 34x34 f3 = 160MHz
Parameters abscissas abscissas 50x50 subdivisions
Time 0.183 2.253 27.91s
Computation Time 1x 12.50x 155.06x
Relative to
Apodized FNM

 

 

 

 

56

pressure evaluated at times t = 1.5625115 and t = 3.0625115.

Time vs. error comparisons for a Hanning-weighted input pulse

Figure 3.9 shows the computation times and the numerical errors for the transient
pressure field generated by the Hanning-weighted input pulse evaluated with the tran-
sient apodized FNM, the transient Rayleigh-Sommerfeld integral, and the Field II pro-
gram. For both the transient apodized FNM and the transient Rayleigh-Sommerfeld
integral, the simplified time space decomposition approach is applied as described in
Section 3.2.5. In Figure 3.9, the transient apodized FNM achieves the smallest errors
for a given computation time, and the transient apodized FNM requires the smallest
amount of time to achieve a given error. Comparisons between the transient apodized
F NM, the transient Rayleigh integral, and the Field II program are summarized in
Table 3.4. To achieve an NRMSE of 0.1 with the Hanning weighted pulse, the tran-
sient apodized FNM needs 8 abscissas in each direction, and the computation time
is 0.64 seconds. The transient apodized FNM achieves an NRMSE of 0.01 in 0.72
seconds with 10 abscissas applied in each direction. To achieve an NRMSE of 0.1 with
the Hanning weighted pulse, the transient Rayleigh-Sommerfeld integral needs 14 ab-
scissas in each direction, and the computation time is 1.48 seconds. The transient
Rayleigh-Sommerfeld integral achieves an NRMSE of 0.01 in 8.57 seconds with 34
abscissas in each direction. To achieve an NRMSE of 0.1 with the Hanning weighted
pulse, the Field II program requires a temporal sampling frequency of 64MHz and 20
subdivisions in each source direction, and the result is obtained in 4.68 seconds. The
Field II program achieves an NRMSE of 0.01 with a temporal sampling frequency
of 128MHz and 50 subdivisions of the source in each direction. This calculation is
completed in 24.04 seconds. Thus, for an N RMSE of 0.1, the transient apodized FN M
is 2.31 times faster than the transient Rayleigh integral, and the transient apodized

FNM is 4.66 times faster than the Field II program. For an NRMSE of 0.01, the

57

normalized pressure
0

     
  
  
 
 

0.8

   

0.6
0.4

0.2

- -2 0

Y (“""s °f 7‘) 2 (units of azm.)

a) t = 1.5625113.

0
L-
: 1
m
(D
d)
h
D.
'U
0
I!
a 0.5
E
h
0
C

’ 0.6 0.8

   

- -2
y (units Of M 0 2 (units of a2I4k)
b) t = 3.0625113.

Figure 3.8. Simulated reference transient field for an apodized square source excited
by the Hanning-weighted pulse in Eq. (3.23) with f0 = 1.5 MHz and W = 2.0/\. The
sides of the square source are equal to 4/\. The apodization function is given by Eq.
(3.22). The transient reference pressure, evaluated in the x = 2.0) plane, is computed
with 100,000 Gauss abscissas in each direction using the Rayleigh integral. Results
are plotted at a) t = 1.562511s and b) t = 3.0625113.

58

 

NRMSE

   

 

”‘APOdized FNM 3 2:51:35 ; '35::
mime-Rayleigh :2
-A-Fieldll “ E 3 33.51 3 £3;.23§

 

 

 

 

 

i

 

 

10'2 10'1 10° 101 102
computation time (seconds)

Figure 3.9. Normalized root mean square error (NRMSE) plotted as a function of
the computation time for transient pressure calculations evaluated with the apodized
FN M, the apodized Rayleigh-Sommerfeld integral, and the Field II program. For the
same computation time, the apodized FNM achieves the smallest errors, and for the
same error, the apodized FNM requires the least amount of time.

59

transient apodized F NM is 11.90 times faster than the transient Rayleigh integral,
and the transient apodized F NM is 24.04 times faster than the Field 11 program.
Once again, these results are specific to the grid, apodization function, and piston
geometry utilized in this paper. Overall, the transient apodized FNM achieves the

best performance in time vs. error comparisons with the Hanning-weighted pulse.

3.4 Discussion

3.4.1 Large-scale computation

Despite the dramatic increase in processor speeds and computer memory in re-
cent decades, nearfield pressure calculations that numerically evaluate the Rayleigh-
Sommerfeld integral are still very slow. This is especially evident in calculations
performed on large grids, which take an hour or more on the fastest modern desk-
top computers. For example, more than one hour was needed in a recent paper to
compute the reference pressure field in the Fresnel zone for a single element using the
Rayleigh-Sommerfeld approach [14]. More recently, nearfield pressure calculations
with the Rayleigh-Sommerfeld approach for a single element required 117 hours [39]
per simulation to achieve a NRMSE of 10’7. Thus, even for simulations of single ele-
ments performed in the nearfield region, traditional computational methods perform
poorly, and faster methods such as the FN M are needed.

Computation times with the Rayleigh-Sommerfeld integral can also be excessively
long for phased array simulations. Extending the results shown in Tables 3.3 and
3.4 and Figures 3.7 and 3.9 to a linear phased array with 1024 rectangular elements
evaluated on a computational grid with 100 x 100 x 100 points in the x, y, and 2
directions, the Rayleigh-Sommerfeld integral will require about 10 hours for a 0.1
NRMSE and about 80 hours for a 0.01 NRMSE. For the same array evaluated on the
same grid, the FNM will require only 2.5 hours for a 0.1 NRMSE and 6.5 hours for

60

Table 3.4. Simulation parameters for transient calculations that achieve normalized
root mean square error (NRMSE) values of 0.1 and 0.01. Parameters listed include
the number of Gauss abscissas or the corresponding Field 11 parameters, the resulting
computation time, and the computation time relative to the apodized FNM for the
Rayleigh integral and the Field II program. (a) For a 0.1 NRMSE and (b) for a 0.01
NRMSE.

(a)

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Transient Calculations
0.1 NRMSE
Apodized FNM Rayleigh Field II
Simulation 8x8 14x14 f3 = 64MHz
Parameters abscissas abscissas 20x20 subdivisions
Time 0.643 1.483 2.983
Computation Time 1 x 2.31 x 4.66x
Relative to the
Apodized FNM
(b)
Transient Calculations
0.01 NRMSE
Apodized F NM Rayleigh Field 11
Simulation 10X 10 34x34 f3 = 128MHz
Parameters abscissas abscissas 50x50 subdivisions
Time 0.723 8.573 17.263
Computation Time 1x 11.90x 24.04x
Relative to the
Apodized FNM

 

 

 

61

a 0.01 NRMSE. Thus, for large—scale computations, the time savings achieved by the
apodized F N M is substantial. The apodized FNM is also an ideal reference for both

time-harmonic and transient calculations.

3.4.2 Time and error comparisons

In addition to the results presented shown in Tables 3.3 and 3.4 and Figures 3.7
and 3.9 that compare the Rayleigh-Sommerfeld integral and the apodized FN M, the
individual computation times and error values were also determined as a function of
the number of abscissas. Those results show that, for a given number of abscissas, the
apodized FNM achieves much smaller errors than the Rayleigh~Sommerfeld integral,
and the computation time of the apodized FNM is slightly less than that of the
Rayleigh integral for the same number of abscissas. Therefore, the rapid convergence
of the apodized FNM is responsible for the enhancement in performance relative to
other methods. However, fair comparisons of these methods require the evaluation
of both times and errors as shown in Tables 3.3 and 3.4 and Figures 3.7 and 3.9.
Otherwise, a method that generates large errors in a similar amount of time or a
method that requires much more time to achieve a similar error could be incorrectly
identified as having comparable performance. The method with superior performance

should consistently achieve the smallest error in the shortest time.

3.4.3 Apodization functions

Eq. (4.24) and Eq. (3.21) can be applied to a large class of apodization functions. The
only requirement for the apodization functions is that 8 f (n, V) /3p and 6 f (a, V) /8u
exist. Other approaches have different requirements. For example, Fresnel approxi—
mations have been obtained for several apodization functions including exponential,
sinusoidal, Gaussian, or error-function apodization functions [14]. One restriction

for these Fresnel approximations is that the surface velocity distribution should be

62

separable such that f (pm) = fflp) fy(1/). The advantages of the apodized FNM
expressions in Eq. (4.24) and Eq. (3.21) are that: 1. separability of the apodization
function is not required, and 2. accurate results are obtained very close to the piston
face.

Pressures generated by a Gaussian apodization function [40, 41, 42], given by
f(p, V) = e-(”_u0)2/20"2‘e_(V_V0)2/203 Where #0, V0, 0“, and 0;; are constants,
were also evaluated with Eq. (4.24). When the Gaussian apodization function is
applied, the entire expressions in Eqs. (4.24) and (3.21) are needed because the
T2,- terms are nonzero in this case. Furthermore, the integration limits of the single
integrals are either (0, a) or (0, b), so the abscissas and weights computed for the
single integrals can also be used for the double integrals. The computation time
needed for all terms in Eqs. (4.24) and (3.21) is only slightly more than that needed
to evaluate only the double integrals, especially when the number of abscissas is larger

than 10 in one direction.

3.4.4 Time Space Decomposition

The expressions for the transient nearfield pressures are computed with the simpli-
fied time space decomposition approach outlined in section 3.2.5 for both the tran-
sient apodized FN M expression and the transient Rayleigh-Sommerfeld integral in
Eq. (3.21) and Eq. (3.2), respectively. The simplified time space decomposition
approach provides an important advantage in that, unlike the original time space de-
composition method described by Kelly and McGough [36], one dimensional and two
dimensional nearfield diffraction integrals are readily computed. The original time
space decomposition method of Kelly and McGough [36] stores a 3D K matrix that
captures the arrival and departure information for a finite duration pulse, where the
transient pressure propagation is described by a single integral F N M expression. Ex-

tending the original time space decomposition to handle transient calculations with

63

the 2D apodized FNM expressions in Eq. (3.21) or the Rayleigh-Sommerfeld integral
in Eq. (3.2) involves replacing the three dimensional K matrix with a complicated
five dimensional K matrix. With the simplified time space decomposition approach,
no matrix storage is needed. Thus, the algorithm complexity and the memory us-
age are greatly reduced. The simplified time space decomposition in section 3.2.5 is
preferred for transient calculations with 2D integrals, where a significant reduction
in the computation time is achieved relative to other approaches without increasing
the numerical error. For transient calculations with 1D integrals, there is a trade-off
between the simplified time space decomposition and the original time space decom-
position method. The original time space decomposition requires more memory, is
more complicated to program, and is slightly faster, while the simplified time space

decomposition is much easier to program.

CHAPTER 4

A 1D Fast Nearfield Method for
Rectangular Pistons with

Polynomial Apodization

Pressure fields generated from apodized rectangular sources are readily calculated us-
ing the Rayleigh-Sommerfeld integral [6] and the analytically equivalent 2D apodized
FNM expression [43]. Both of these methods evaluate double integrals, which handle
various types of apodization functions. However, the numerical performance of those
methods are quite different. Numerical calculations of the Rayleigh-Sommerfeld in-
tegral converge very slowly in the nearfield region due to a 1 / R singularity. The
performance of the 2D apodized FNM is much better than the Rayleigh-Sommerfeld
integral in the nearfield since the singularities have been eliminated. Thus, the 2D
apodized FNM is a better choice if a general 2D apodization function is given. How-
ever, the computation time can be reduced further if the 2D FNM integral can be
converted into a 1D FNM integral. This is achieved with polynomial apodization
functions.

To improve the numerical performance of nearfield calculations with polynomial

65

apodized rectangular pistons, a 1D polynomial apodized F NM expression is derived
from the instantaneous impulse response proposed by Scarano et al [44]. The deriva-
tion of the apodized F N M for rectangular pistons starts by convolving the instanta-
neous impulse response with the polynomial apodized particle velocity. Then, after
obtaining a double integral and specifying the particular polynomial apodization func-
tions, the inner integral can be integrated analytically, and the polynomial apodized
F NM is obtained. Two polynomial apodization functions are considered, where one
is a 1D quadratic polynomial and the other is the product of two 1D quadratic poly-
nomials. The results show that the convergence of the polynomial apodized FNM is
much better than the 2D apodized FNM and the Rayleigh-Sommerfeld integral for

both of these apodization functions.

4.1- Polynomial apodization derivation

4.1.1 Instantaneous Impulse Response

The instantaneous impulse response [44] is obtained at a fixed time t on a specific
plane where the value of z is a constant. The spatial particle velocity distribution

3(33, y) on the piston face is given by

s(:r, y) = u(:r)rect(:c/a)w(y)rect(y/b) (4.1)

where u(a:) and w(y) are the normal particle velocity distributions on the vibrating
surface along the x and y directions, rect(x) = 1 for [2:] S 1 and rect(x) = 0 for
|:c| > 1 and a and b are the width and height of the rectangular piston in the :1: and
y directions, respectively. The separable particle velocity distribution can also be

expressed by means of two-dimensional convolution as

8(17, y) = (11(I)I‘ect(1‘/ a')<5(y)) * *(w(y)rect(y/b)5(l*)) (4-2)

66

The impulse response h(z; t) is given by

6(1‘ — 7‘0)

h(z; t) = %u(x)rect(x/a)6(y) * * w(y)rect(y/b)6(z) * a: r0

6
= 27141? ereCWB/a) (y )* *f(I 31) (43)
where r = V152 + 312, and To = V c2t2 — 2:2. The instantaneous impulse response for

y > b is given by [44]

0 r0<y—-b

C9K,€h%K—£—) y-bSr Sy+b
h(z;t)= E 1 2 (Hal 0 , (4-4)

7r 29([6 {1, £2) )(rect (aT—il)

 

 

—,(g3 gummy—$12)] r0 2 y + b
_Owh/TQ-rz-y)
where 9((61,€2)::_f€12u<003-062)”? d6, $1 2 «7‘3 — (y — b)2,

 

$2 = \/""(2) - (y + (1)2151: max(:1: ‘ a, -131), 63 = max(:z: — a, —:1:2),
52 = min(:z: + (2,131), {4 = min(:1: + a, 1:2).

The apodized pressure based on these expressions is readily extended to all values

of 3].

4.1.2 Time-harmonic pressure calculations

The time-harmonic pressure generated by a polynomial apodized rectangular piston

is given by

+oo
P(:1:,y,z;k)=jwpejwt / e_JWTh(z;’r)dT (4.5)

—00

67

where w is the excitation frequency in radians per second, p is the density of the
medium, h(z; t) is the instantaneous impulse response in Eq. (4.4) which accounts for
the apodization, k is the wavenumber, and t is the elapsed time. Distance parameters
are defined as (11: \/(x—a)2+(y-b)2+22, d2 = \/(x—a)2+(y+b)2+z2,
d3 = \/(:c+a)2+(y—b)2+z2, and d4 = \/(x+a)2+(y+b)2+z2, and then

 

 

 

 

h(z;t) in Eq. (4.4) is substituted into the pressure calculation expression in Eq.

(4.5). For :1: 2 a, the expression for the nearfield pressure is obtained,

6 3 "12' —jkfi "22'
P(:1:,y, z; k) = ? Z (112' e f(€7fi)d€d!87 (4-6)
7r ’ m . m .
2:1 12 22
where the values of qz- are qll = q12 = 1 and q13 = —1, the values of (mlz‘fllli)

are (m11,n11) = (d3,d4), (m12,n12) = (d1,d3), and (m13,n13) = (d2,d4), the

values of ("12237120 are (7n21,n21) = (:1: — a,:r + a), (m22,n22) = (:1: — man), and

(m23m23) = (I - 0,12), and f(€,fi) = Uf-T - {WM/T3 - $2 - IDA/32 - 22 - £2-
The inner integrals £7,322 f ((5 , fi)d§ in Eq. (4.6) can be represented as the sum of

 

two integrals,

712.,- 7122: 0
/ mmdg = / mm: + f(€,fi)d€- (4.7)
"122' 0 "122'
Using Eq. (4.7), the pressure expression for 2: 2 a and y 2 b in Eq. (4.6) forms six

double integrals which are given by

3 711' '. 712'
P($,y.z;k) = 56;; Z 011/ z fjwf Z f(€,5)d€dfi
i=1 77212“ 0
3 711' . 0
5; q1z‘ / if?“ f(£,fi)d€d/3 (4.8)
i=1 "‘12” "122'

And by combining the limits in the outer integral, Eq. (4.6) can be simplified as

713' _-,. 714'
P<x.y,z;k> = £21127: / 2e 43/0 WWW/3, (49)

23:1 m3i

68

where the values of qz- are gm 2 q24 = —1 and q22 = q23 = 1, the values of
(manner) are (m3iin31) = ((11,612), (m32vn32) = (613,614), (m33m33) = ((11,013)
and (m34,n34) = (d2,d4), and the values of 17.4,- are 7241 = :1: - a, n42 = a: + a,
1143 = 1:1, and 1144 2 51:2, and f(§,fi) = u(:1: — £)w( r8 -— x2 — y)/\/fl2 — 22 — £2.

 

4.1.3 1D quadratic apodization

In the derivation of the FNM expressions for 1D quadratic apodization, the apodiza-

2

tion function u(:c) is a quadratic function with u(x) = p + qx + r1: , where p,

q, and r are constant coefficients and w(y) = 1. Then, u(:r - 5) is written as
u(:1:—-§) = A+B€+C§2, where A =p+qx+rx2, B = —(q+2rx), and C: r and

w( r8 - 2:2 — y)=1. The inner integral in Eq. (4.9) is computed as

£2 dg
/ “(5’3 ’ 0““— : Iterml + Iterm2 + Ite'rm3’ (4'10)
0 7.2 _ £2
0
where Iterml = (A + 2CT3) tan—1 (62/ V T8 _ 6%), 11:8er 2

— (B + %C€2) [Hg — €2,1term3 = By/fl2 — .22, and r0 = V52 - 22.

The first term

After substituting the first term from Eq. (4.10) into Eq. (4.9), the tan—1(-) term
is converted into an equivalent integral expression using the identity tan‘1(:z:) +

tan"1(1/:1:) = 7r/2. The tan—1(-) term then becomes

:r—a

1
\/fi2-22-($-a)2

tan

 

 

 

 

2—22— ,H, 2 _ a
g- /0\/fl ( ) a: do. (4.11)

(a: —a)2+02

69

By substituting Eq. (4.11) into Eq. (4.9) and exchanging the order of integration,

the pressure calculation result for Iterml in Eq. (4.10) becomes

 

Pterm1($,y,2;k)2:2(131/"3/29w:2?sz , (4.12)

where the values of the q,- are q31 = q33 = 1 and q32 = q34 = —1, the values of 0,- are
01: J02+(x—a)2+22, 02 = \/02+(x+a)2+22, 03 = \/02+(y—b)2+22,

and a4 = \/02 + (y +b)2 + 22, the values of (mt-,ni) are (m1,n1) = (m2,n2) =

 

 

 

 

(y — b,y + b) and (m3,n3) = (m4,n4) = (a: — a,a: + a), and the values of h,- are
hl =x-a, h2=x+a, h3=y—b, and h4=y+b.

The four integrals in Eq. (4.12) are double integrals, but the inner integral with
respect to ,6 can be evaluated analytically. The resulting 1D integral for Eq. (4.12)

is given by

 

(-Ak2 C)h- _-k. _-k
PteTm1(:c,y,z;k) 2:2q3i/:i|:jk3(h2+02)2(6 .7 Uz—e J z)

 

+jC C }- _-- . _- ..
+2—19hi63-Jj2kaz-[pl2lz 2(01'8 Jkal—Ze 1164)] d0. (4.13)
z- +0
2

The second term

Substituting Iterm2 into each double integral in Eq. (4.9) and performing the change

 

the variable 0 = \/132 — 2:2 —— 6% produces the analytical expression given by

—jk0‘- 2

e 20

Pterm2<T y z; k) _ ‘78 i§:1:q2z/:z(_7—) d0
2

+11) (6-3-m3 _ e—jkdl) _ y_fb (e—jkca _ (M2)

70

4 n- —jk0’z‘ 2
~22 Q2'ihz' / Z (:4) da, (414)

. m. 0'
2:1 1

z
where the values of the (122' are gm 2 q24 = —1 and (122 = q23 = 1.

The third term

After substituting Iterm3 into Eq. (4.9), the integrands of all four double inte-

grals are the same. The first and fourth double integrals can be combined to form

 

a new double integral with the integration limits from \/(:c — a)2 + (y - b)2 + 22

 

to \/(a: + (1)2 + (y + b)2 + 22 with the integrand —e‘jkflB £32 - 22, and the sec-

ond and third double integrals are combined into a double integral with the inte-

 

 

gration limits from \/(x — (1)2 + (y — b)2 + 22 to \/(x + a)2 + (y + b)2 + 22 for the

integrand e’jkfiBVfi2 — Z2, thus

Pterm3($v y, z; k) = 0. (4.15)

The final equation for the 1D quadratic apodization

The total pressure in Eq. (4.9) is the sum of three terms given by Eq. (4.13), Eq.
(4.14) and Eq. (4.15), where P(;r, y, z; k) = Pterm1($a y, z; k) + Pterm2(33v'ya z; k) +

Pterm3(xv y, z; k), which is

 

4 n- . 2
c z 3(Ak -C' h- _- . _-
P($,y’z;k) : 57.; Z q3'i/I [ ) 2 (e JkO'z _ 6 JkZ)

 

(32' —'k- C h' —'k- —'k
_ . J 0 _ l . .7 0 _ .7 Z
+2khze 2+k2h22+02 (ale 2 ze ) d0
2 n- —jk0" 2
c z e 20
i=1 "‘2' Z

71

 

+§7cl3 (e—jkd3 _ e—jkdl) _ yjjfkb (e—jkd4 _ e—jkd2)]

4 n- —jk0’z' 2
——£ (12,12,- / z (g___0_) do. (4.16)
47r , m-

z=1 2 of
4.1.4 2D Quadratic apodization

In the derivation of the FN M expression for 2D quadratic apodization, the apodization

functions u(a:), w(y) are quadratic functions with u(:1:) = p1 +qla: +r1$2 and w(y) =

P2 + (1211 + 7‘2y2

resulting expressions for u(:r — 5) and w(‘ (7(2) — {2 — y) are

, where p1, ql, 7‘1, 192, qg, and Q are constant coefficients. The

11(11: — 5) = 141+ Blé + 0162

rah/1‘3 — £2 — y) = A2 + Bm/rg —§2 +C2(rg —§2)

Where A1 = P1 + (1193 + 7‘12, Bl = -(41 + 27130.01 = T1, A2 = 102 - (121/ + r2312,

32 = (p2 - q2y), 02 = 7‘2, and r0 = (M32 — 22.

The inner integral in Eq. (4.9) is computed as

52 urn/r3 - a?) - M
f “(33 - E) 2 2 = Iterml + Iterm2 + Iterm3
0 ‘/TO _g

where Iterml = §(C1C'2T3 + 401242712) + 4A1C2Tg + 8141/12) tan (7%) ,
r —5
0

1
Iterm2 = \(7'8 — {3862(4A1C2—CIC2T3+20102§%—4A2C'1)-VT?) —- 62(311‘12-1-
:1;3102(T(2) - £2» and Iterm3 = 1130182§3 + %BIBQ{% + A1326 .

 

(4.17)

The first term

The Iterml is rewritten as

72

T(,8) = t1fi4 + be? + (£3
where t1 = §0102, t2 = gay/12 + A102) — $010222 and t3 = $010254 —
$42142 + A102)22 + .4142.
Substituting the first term Iterml from Eq. (4.17) into Eq. (4.9) and transforming
tan—1(-) according to Eq. (4.11), the pressure calculation result for the first term

Ite'rml in Eq. (4.17) becomes

 

_' h
Pterm1—2D($ 31,2 k) - 2132/ ‘/(:T(fi ejkfldfi 2 2 d0 . (4.18)

2
hi + 0
where q31— — q33— —- 1 and q32 = q34- _—
The inner integral in Eq. (4.18) is computed analytically

o . o .
/ e—Jkflrmw = e—Jkflulfi‘l + :25? + t3)dfi
02' ‘72"

1 _ - .
—_— fie 9W (T1 44 + T253 + T362 + T4fi + T5)|32.

where T, ___ 411:4, T2 : 4jk3t1, T3 = 12k2t1, T4 = 2jk3t2 — 243441, and T5 =
2t2k2 — t3k4 — 2441.

The first integral in Eq. (4.18) is given by

Pterml—2D($v y, Z; k) =

 

712‘
. 5 Z 932‘] [Tlhz' (02 + (12)8—jk01 + h (2le2 +T202' +T3)e— jkaz’
—]k 271' 2:1 ml-
e—jkai _ e—jkz

hzz+a2

02'8 —jk0i_ Ze—jkz

+11 T252 +T
1(2 4)2 U244}?

+ (T124 + z2T3 + T5)h,-

 

 

(4.19)

73

The second term

Substituting Iterm2 in Eq. (4.17) into each double integral in Eq. (4.9) and perform-

 

ing the change the variable a = \/B2 — 22 — 5% produces the analytical expression

given by

Pterm2—2D(Iv 3” Z; k) =

4 mi —jk0- 1 2 2
Z (121' m e 2 g0h.i(4A1C'2 -— 414201 + 010'th- — 01020 )
i=1

1 2 / 2 2 1 2 2 U
—a(BlA2+§BlC20 )+ a +hi(B1A2+§-(a +hz')) U—ida, (4.20)

where the values of the (122' are gm 2 gm = —1 and gm 2 q23 = 1.

The third term

Substituting Iterm3 in Eq. (4.17) into each double integral in Eq. (4.9) and perform—

 

ing the change the variable a = \/[32 -- 22 — 5% produces the analytical expression

given by
Pterm3—2D($’ y, Z; k) I
:31? (go B2h113+éBlBgf€+ A1B2h1) (‘3—jkd2 " e—jkdl)
+327 (£013.27)? + $131th? + AlBth') (e_jkd4 ‘ ‘3—jkd3)

1:1;013203 + %BIBZU2 + A1320) 0

 

n. .A (
+22-lz3q2i/ ze—JL’Ul 0- do. (4.21)
7

nz- z

74

The final equation for the 2D quadratic apodization

The total pressure in Eq. (4.9) is the sum of the three terms given by Eq. (4.19), Eq.
(4.20) and Eq. (4.21), thus

P2D(xay12; k) = Pterm1—2D(x7 ya Z; k)

+Pterm2—2D($w y, z; k) + Pterm3—2D(xiy’ Z; k)' (4°22)

4.1.5 The 2D apodized FNM

To obtain the formula for the 2D apodized F N M, origin of the coordinate system needs
to be shifted from the center of the source to the lower left corner of the rectangular

source. The observation point (331,311, 21) is then defined as

xl=x+a
31:2.

The time-harmonic pressure generated by an apodized rectangular source computed

with the 2D apodized F N M is now given by

- 2 2 2 .
2a2b e_]k\/ali+zl+h12 e-Jk21

 

 

jwt _
. _ _pcr_)__0:
Papod($1,yl,21,k) — Z/[Tlihlz‘e 02 +h2 dfldV
2:10 0 12' 12'
_ / 2 2 .
_p—2__a10:th4jk a2z+z1+h2i_ €_‘7 1921
Z / TQihgze 2 do (424)
a .+h2 -
i=1nz 22 22
where the values of le' are T11 = T12 2 _0f 1;,1/ and T13 = T14 = _0f ’Z’V ,

the values of ah- are 0‘11 = (112 = p — 2:1 and 0:13 2 0114 = 1/ — yl, and the values

75

of h12’ are hll = y1,h12 = z/ — y1,h13 = 2:1, and h14 = p - $1. The values
of T22- are T21 = T22 = f(0,2b) and T23 = T24 = f(2a,0), the values of 022- are
0121 = 0122 = 0 - 2:1 and (123 = 0224 = a — yl, the values of h22’ are I221 = 3,11,
h22 = 2b — yl, I223 = x1, and hu 2 2a — 2:1, and the values of (ml-,ni) are
(m1,n1) = (7722,222) = (2a, 0) and (m3,n3) = (m4,n4) = (2b, 0) for 2' = 1,2,3,4.

4.1.6 The Rayleigh-Sommerfeld integral

The time-harmonic pressure generated by apodized rectangular source is also com-

puted with the Rayleigh—Sommerfeld integral [6] via

. 'wt 2a 2b _ 'kR
pr'voe] 6 J
PRayleighC‘lJ/lazvk) = T f(#,I/) R dpdu, (4.25)
0 0

 

 

where R = \/(a:1 — m2 +(y1 — V)2 + 2% is the distance between the observation

point (x1, y1, 21) and the source point (a, V, 0).

4.1.7 Error Calculations

Let P(:L', y, z; k) represent the computed pressure field, and let Pref(:r,y, z; k) repre-

sent the reference pressure field. The peak normalized error is given by

|P(.’L', 1 a 2:; k) _ Pref(xvyazi k)lmax(x_y,z)

 

777an = (4.26)

[13”?me y: 27; k)|max(;r,y,z)

76

Table 4.1. Simulation parameters that achieve peak normalized error values of 10%,
1%, and 0.1% for the 1D quadratic apodization function 21(2) = 122 — 2:2. Parameters
listed include the number of Gauss abscissas, the computation time, and the compu-
tation time relative to the polynomial apodized FNM for the Rayleigh—Sommerfeld
integral.

Peak Normalized 10% 1% 0.1%

Error

 

 

FNM Rayleigh FNM Rayleigh FNM Rayleigh
Gauss abscissas 8 32 X 32 10 74 x 74 11 122 x122
Time (seconds) 0.0761 1.4147 0.0924 7.5038 0.0991 20.3334

Computation time 1x 18.58x 1x 81.17x 1x 205.11x

relative to the

 

 

 

polynomial
apodized FN M

 

 

 

 

 

 

 

 

 

4.2 Results

4.2.1 1D quadratic apodization
Reference pressure field

Simulation are performed with a rectangular source, where each side is equal to 2a =4
wavelengths. The origin of the coordinate system is the center of the rectangular
source. The reference pressure field is computed for the apodization function u($) =
2:2 — a2 using the Rayleigh-Sommerfeld integral in Eq. (4.25) with 10,000 Gauss
abscissas in each direction. The reference pressure shown in Figure 4.1 is evaluated in
the y = 0 plane. The grid in the x direction extends from 0 to 4 wavelengths, and the
grid in the z direction is evaluated from 0.04 to 4 wavelengths. The sample spacing
in the x direction is 0.05 wavelengths and the sample spacing in the z direction is

0.04 wavelengths. The total number of grid points is 8100.

77

u
\\
‘.|‘“\|““\

”’2!

[II II 3

III '0 5‘“; “h‘

1,": ’II’I 1"": ”o o‘ ‘i\\‘\‘\ \‘\ “i“

, ”I I I I Q

' ’o’Ifl, 07:2, ,2, ‘0‘“
"’III $0....“

“‘|‘ :\
‘3‘ i“:
\\\:“n‘ 3“: ““5 \ ‘\“““\““

\\\.
‘\\
02.035 t‘ \‘ \\\\\\“\‘“
13:01.: not ““3“
o

‘ s
\\\\u\“ \\\ ‘ ‘\\\:“““ ““5. “no
\‘ \\\\\‘ ‘ r“ ‘
\‘\\‘\“ “N ‘\\ ‘\‘
\\\ “

normalized pressure

 

\\\\\‘\‘\
\\\“‘\‘ \\\“
\\\

2

 

x (units of k) 2 (units of A)

Figure 4.1. Normalized simulated time-harmonic reference pressure field in the y = 0
plane for a rectangular source with each side equal to 4 wavelengths. The pressure
field is computed with 10, 000 Gauss abscissas in each direction using the Rayleigh-
Sommerfeld integral with the 1D quadratic apodization 22(1) = x2 — a

Error comparisons and computation times

The numerical errors and computation times for the polynomial apodized FNM and
the Rayleigh—Sommerfeld integral are given in Figure 4.2. Figure 4.2 (a) shows that
for a given number of Gauss abscissas in one direction, the polynomial apodized
F NM always achieves a smaller error than the Rayleigh-Sommerfeld integral. The
polynomial apodized FNM achieves a maximum error of 10_14 with less than 80
Gauss abscissas while the Rayleigh-Sommerfeld integral only achieves a maximum
error of 10_5 when 200 X 200 Gauss abscissas are applied. Figure 4.2 (b) shows the

computation times for the polynomial apodized FNM and the Rayleigh-Sommerfeld

78

 

 

— 1D apodized FNM
'- ' - ' Rayleigh

 

 

 

normalized error

 

 

 

 

o 20 4'0 610 810100120140150180200
number of abscissas

a) Maximum errors.

 

G
O
.4

 

—1D apodized FNM .1
...... mmRayleigh I‘d

0|
O

 

 

 

.h
C
‘1
s

N

G

$
s

computation time (seconds)
A w
G O
3‘ :
a
5 5‘

 

 

 

O . L ; i
0 20 40 60 80 100 120 140 160 180 200
number of abscissas

b) Computation times

Figure 4.2. a) Maximum errors and b) computation times for the polynomial apodized
FN M and the Rayleigh-Sommerfeld integral.

79

 

-5

.........................................

 

 
   

 

_ —Polynomia| FNM
_ -----Rayleigh

normalized error

 

 

 

 

-15 Lg jig};

0.01 i 0.1 1 To 100
computation time (seconds)

"Ali i i;i""i

 

10

Figure 4.3. Time vs. error comparison between the polynomial apodized FNM
and the Rayleigh-Sommerfeld integral. For the same computation time, the polyno-
mial apodized FNM achieves smaller errors, and for the same error, the polynomial
apodized FNM requires less time.

integral as a function of N, the number of Gauss abscissas in each direction. in Figure
4.2 (b), the polynomial apodized FNM consistently requires less time for the same
value of N.

Figure 4.3 shows the computation times and the errors evaluated with the 1D
polynomial apodized FNM and the Rayleigh-Sommerfeld integral on the same plot.
In Figure 4.3, the polynomial apodized FNM achieves smaller errors for the same com-
putation time, and the apodized FN M uses less time to achieve the same error value.

Comparisons between the polynomial apodized FNM and Rayleigh-Sommerfeld inte—

80

gral are summarized in Table 4.1. For 10% peak normalized error, the polynomial
apodized FNM is 18.58 times faster than the Rayleigh-Sommerfeld integral, for 1%
peak normalized error, the polynomial apodized FNM is 81.17 times faster than the
Rayleigh-Sommerfeld integral, and for 0.1% peak normalized error, the polynomial
apodized FNM is 205.11 times faster than the Rayleigh-Sommerfeld integral. These
values were obtained for the grid, quadratic apodization function, and piston geome—

try used here.

4.2.2 2D Apodization function

The pressure field is also evaluated for a square piston with an apodization function
22(1):, y) = (:2:2 - a2)(y2 — oz). The reference pressure is computed using the Rayleigh-
Sommerfeld integral in Eq. (4.25) with 10, 000 Gauss abscissas in each direction. The
reference pressure shown in Figure 4.1 is computed in the y = 0 plane. The grid in
the a: direction extends from 0 to 4 wavelengths, and the grid in the z direction is
evaluated from 0.04 to 4 wavelengths. The sample spacing in the a: direction is 0.05
wavelengths, and the sample spacing in the z direction is 0.04 wavelengths. The total

number of grid points is 8100.

Error comparisons and computation times

The numerical errors and computation times for the 1D polynomial apodized F NM,
2D apodized F NM, and the Rayleigh-Sommerfeld integral are given in Figure 4.5.
The apodization function is a 2D function given by u(a:,y) = ($2 — a2)(y2 — a2).
Figure 4.5 (a) shows that, for a given number of Gauss abscissas in one direction,
the computation time of 1D the polynomial apodized FNM is always smaller than
the 2D apodized FNM and the Rayleigh-Sommerfeld integral. The 1D polynomial
apodized F N M converges very rapidly, and an error of 10‘14 is achieved with less

than 80 Gauss abscissas. For the same piston and grid geometry, the 2D apodized

81

Table 4.2. Simulation parameters that achieve peak normalized error values of 10%,
1%, and 0.1% for the 2D apodization function 22(23, y) = (x2 —a2)(y2 —a2), where a 2
2 wavelengths. Parameters listed include the number of Gauss abscissas, computation
time, and the ratio of the computation time relative to the polynomial apodized

FNM. (a) 10% peak normalized error, (b) 1% peak normalized error and (c) 0.1%
peak normalized error.

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

(a)
Peak normalized error 10%
1D FNM 2D FNM Rayleigh-Sommerfeld
Gauss abscissas 10 8x8 34x34
Time (seconds) 0.0816 0.1306 1.1828
Computation time relative to the 1 x 1.60x 14.49x
polynomial apodized FNM
0))
Peak normalized error 1%
1D F NM 2D F N M Rayleigh-Sommerfeld
Gauss abscissas 10 14 x14 78x78
Time (seconds) 0.0816 0.3708 6.1866
Computation time relative to the 1x 4.54x 75.78x
polynomial apodized FNM
(C)
Peak normalized error 0.1%
1D FNM 2D FNM Rayleigh-Sommerfeld
Gauss abscissas 12 28x28 126x126
Time (seconds) 0.0957 1.3986 16.3940
Computation time relative to the 1x 14.61x 171.25x
polynomial apodized FNM

 

 

82

 

 

 

u
‘m “‘ ‘
“u \
‘\\\\

\\\l

‘\
u
\\
\‘l‘\““‘\‘\
\\\\\“‘\\ ‘
\\\\“‘“\

\ u \

u“ \\\
u“ u u“ \
‘\\\:u““‘.\\u\\\\\\n\“‘\‘u

\ \\\“ u\ “u
“n“ “\\\u“ (m “\\\
\\\“ \\\ \\\ “u
u m

opoo
can-hmoo-s

, I
\
\|“‘\\‘ \
m“ “\\\\‘\‘\‘\\‘\\““
\\\l‘\‘\\\\\“““‘\‘
““\‘-‘\‘u\\\\\\\\\
u\\\|u\\\‘“““‘ \
\ ‘uu\uuu\\‘“‘“\\\“
“\\\““““\\\\\\\‘\‘\“\‘\\\\\
|‘ \\\\\\\\\\u\ ‘ m “
“\l ‘ “um“ ‘
l “um“ \\\uumumu
‘\\“‘“\\\\\\\\\\\\\umuum

\\ \

I“ m
\‘ u
\\““\\“‘|‘

normalized pressure

2

 

x (units of 7L) 5 0 2 (units of A)

Figure 4.4. Absolute value of the simulated time-harmonic reference pressure field
in the y = 0 plane for a square source with each side equal to 4 wavelengths. The
apodization function is a 2D function given by u(:r, y) = (x2 — a2)(y2 — (12), where
0 =2 wavelengths. The pressure field is computed with 10, 000 Gauss abscissas in
each direction using the Rayleigh-Sommerfeld integral.

FNM achieves an error of 10_6 with 100 X 100 Gauss abscissas, and the Rayleigh-
Sommerfeld integral acheives an error of 10‘5 when 200 X 200 Gauss abscissas are
applied.

Figure 4.2 shows the computation times and the errors evaluated with the 1D
polynomial apodized FNM, 2D apodized FNM and the Rayleigh-Sommerfeld integral
on the same plot. In Figure 4.2, the 1D polynomial apodized F NM converges to the
smallest errors in the least amount of time. Comparisons among the 1D polynomial
apodized FN M, the 2D apodized F NM and Rayleigh—Sommerfeld integral are sum—

marized in Table 4.2. For 10% peak normalized error, the 1D polynomial apodized

83

 

 

 

 

—1D Polynomial FNM
; 3 ; '-'-'Rayleigh
10° s ~-§- ------ i ------- - - -2n Apodized FNM ,

normalized error

 

 

 

 

 

 

 

 

 

 

0 20 40 60 80100120140160180200
number of abscissas
a) Maximum errors

80 l ' .' .r l .' F F F
1'; —1D Polynomial FNM i i E I
U , : : : ’
g '-'--Raylergh § § i f
860... -A--2‘DApodizedFNM "; ......
g ,é'
:2 40. ............... ,’_ ’1 ‘\’4
2 i E x i .2":
E s i .‘g' '
3.20. ”’ "\f“..'. 2
E a r," i '5
8 . . ,.;.;.~“~‘ ; . .

0 . ‘,,.' _ . . . .

 

 

 

 

0 20 40 60 80 100 120140 160 180 200
number of abscissas

b) Computation times

Figure 4.5. a) Maximum errors and b) computation times plotted as a function of
the number of Gauss abscissas for a 2D apodization function with the 1D polynomial
apodized F NM, the 2D apodized FN M and the Rayleigh-Sommerfeld integral.

84

 

 

 
  
 

 

 

normalized error

; ? —1D Polynomial FNM
'-'--Rayleigh
' - - -2D Apodized FNM

 

 

 

 

A A L

 

0.01 0.1 1 10 100 1000
computation time (seconds)

Figure 4.6. Time vs. error comparison between the 1D polynomial apodized FNM,
2D apodized F NM and the Rayleigh—Sommerfeld integral. The polynomial apodized
FNM converges to smallest errors with least time.

FNM is 1.6 and 14.49 times faster than the 2D apodized F NM and the Rayleigh-
Sommerfeld integral, respectively. The 1D polynomial apodized FNM is 4.54 and
75.78 times faster than the 2D apodized F NM and the Rayleigh-Sommerfeld integral
for a 1% peak normalized error, respectively. For 0.1% peak normalized error, the
polynomial apodized FNM is 14.61 and 171.25 times faster than the 2D apodized

FNM and the Rayleigh-Sommerfeld integral for the grid, 2D apodization function,

and piston geometry used in this paper.

85

4.3 Discussion

4.3.1 Advantages and disadvantages

The 1D polynomial apodized F NM successfully reduces the two-dimensional apodized
pressure calculation integral in Eq. (3.1) into a single integral in Eqs. (4.16) and
(4.22). As shown in Section 4.2, the performance of the 1D polynomial apodized FN M
outperforms the 2D apodized F NM and the Rayleigh-Sommefeld integral for peak
normalized errors of 10%, 1% and 0.1%. Furthermore, due to the rapid convergence
of this method, the 1D polynomial apodized FNM is an ideal method for reference
pressure field calculations.

The 1D polynomial apodized F NM equations in Eqs. (4.16) and (4.22) are applica-
ble to polynomials of order up to two. When the order of the polynomial apodization
function is increased, the polynomial apodized FNM formula needs to be rederived.
For 2D apodization functions, the current polynomial apodized FNM only handle
apodization functions that are separable in the spatial coordinates. For example,
the apodization function f (u, V) = V112 + V2 can not be evaluated with expressions
in Eqs. (4.16) and (4.22). For this apodization function, the 2D apodized FNM

expression is required.

4.3.2 Interpolation of the apodization function

The convergence of the 1D polynomial apodized FNM is very fast compared to the
2D apodized F N M and the Rayleigh-Sommerfeld integral. Though equations (4.16)
and (4.22) are specialized for the polynomial apodization functions, other apodiza-
tion functions can be approximated with polynomials. The most commonly used
apodization functions are sin(-), cos(-) and exp(-), which can all be approximated
using polynomials within a certain input range of parameters. For these apodiza-

tion functions, expressions for the 1D polynomial apodized FNM with higher order

86

polynomial apodization functions are needed.

4.4 Conclusion

A 1D polynomial apodized fast nearfield method for calculating nearfield pressures
generated by a polynomial apodized rectangular piston is derived and evaluated.
The derivation is based on the instantaneous impulse response for rectangular piston
given by Scarano et al [44]. The 1D polynomial apodized F NM expressions for two
apodization function, u(:r) = 11:2 —a2 and u(:r, y) = (:1:2 —az)(y2 —a2), are given. The
results of the 1D polynomial apodized F N M are compared with the 2D apodized F NM
and the Rayleigh-Sommerfeld integral. The results show that the 1D polynomial FNM
achieves the best performance among all of the three methods for the same number

of Gauss abscissas. Thus, the 1D polynomial apodized FNM is an ideal method for

calculating nearfield pressure generated by a polynomial apodized rectangular source.

87

CHAPTER 5

A Fast Nearfield Method for the
Numerical Evaluation of 3D

Potential Integrals

The evaluation of potential integrals over a volume domain are important calculations
in scattering problems. For example, Moraros and Popovic [23] optimize the volume
potential integrals involved in the moment-method analysis of 3D dielectric scatters.
The precorrected-FFT solution is proposed by Nie et. a1. [45] to solve the volume-
surface integral equation for scattering from conducting-dielectric objects. Potential
integrals are often singular and thus direct evaluation potential integrals may en-
counter numerical difficulties. For potential integrals with singularities, singularity
subtraction [23, 25, 26, 27, 28] or singularity cancellation methods [30, 31, 32, 33]
are often used. Both singularity subtraction methods and singularity cancellation
methods retain the same number of dimension over the integration as the original
potential integrals. Typically, those methods manipulate the integrands to eliminate
singularties and to achieve better performance. This approach is reasonable for gen-

eral potential integrals; however, when uniform potential integrals are considered,

88

more efficient methods can further improve the numerical performance by reducing
the number of dimension over the integration.

Pressure fields from uniformly excited planar sources can be computed using both
the Rayleigh—Sommerfeld integral [6] and the fast nearfield method (FNM) [8, 12].
Those two methods are mathematically equivalent but the numerical performance
differs. Of these, the FNM achieves a smaller error in less time because the FNM
eliminates the 1 / R singularity in the Rayleigh-Sommerfeld integral by subtracting
the singularity in the integrand. The integrands of the potential integrals and the
pressure integrals are the same. For potentials from uniform planar sources, the FNM
can be directly applied. However, the FNM expressions for uniformly excited volume
sources are also needed.

To improve the performance of potential integral calculations for volume sources, a
single FN M integral is derived from the FNM expression for uniformly excited polyg-
onal pistons. The derivation begins after the relationship between FNM expressions
and potential integrals is established. Then, FN M expressions for potential integrals
are obtained for a polygonal source. The prism is subdivided into five subdomains,
and the total potential is obtained by superposing by the potential produced by all of
the subdomains. After substituting the FN M expression for polygonal source into the
potential integral for the subvolume, performing a change of variable and exchanging
the order of integration, a 1D F N M volume integral is obtained for a volume source
that is uniformly excited. The resulting 1D FNM expression for the potential and the
singularity cancellation method are then evaluated on six single observation points
and on a large 3D grid. The results indicate that, when compared to calculations
performed with the singularity cancellation method, the FNM requires fewer sample

points to achieve a given accuracy.

89

5.1 FNM Calculations for Rectangular and Trian-

gular Sources

5.1.1 The potential integral

The time-harmonic pressure produced by a triangular transducer can be obtained

using the Rayleigh-Sommerfeld diffraction integral

, —jwt e—ijds
Pm, 2.2:): MW”: [[3 (5.1)

where s is area of the triangular source, to is the excitation frequency in radians per

 

 

second, p is the density of the medium, v0 is a constant normal velocity evaluated at
the surface of the triangular piston, k is the wavenumber, t is the elapsed time and
R is the distance from the source point to the observation point.

The potential integral considered is given by

<I>(a: ,=y,z) [IA/(31:2RW (5-2)

where V is the integration domain over which the integral is evaluated. The potential

 

in Eq. (5.2) is normalized with respect to —voejwt. The integrands in Eq. (5.1) and
Eq. (5.2) are the same, but the former is evaluated over a surface and the latter is
evaluated over a volume. Thus, the F N M expression that is analytically equivalent

to Eq. (5.1) can also be extended to solve for the potential generated by a 3D source.

5.1.2 FN M calculations for a planar source

Let AABC denote a triangular source, and a, b, and c are the length of the corre-
sponding sides that are located oppsite angles A, B, and 0. According to the fast
near-field method (FN M) [12], the potential integral evaluated over this triangular

source at the observation point above the point A is given by

90

¢@w,z)

 

 

d0,

=jOCe e—jkflds =—1—z figzgai£_e —jkt/;Z¥§7IT§_ e—jkz
47rR 47rk a2—bi__—02 02 + l2
(5.3)
where (a2 — b2 —— 02)/(2c) and (a2 — b2 + 62)/(2C) are the lower and upper limits of
the integration, respectively, 2 is the the orthogonal distance between the observation
point and the source plane, and 1 represents the orthogonal distance between the
vertex A and the base of the triangle (i.e., the height of the triangle).

The FNM expression for a planar source at an observation point (x,y,z) is ob-
tained by first subdividing the planar source into several triangular subelements that
share the vertex D that is defined by the projection of the observation point onto
the planar source. Then, the pressure generated by each subelement is computed
using Eq. (5.3), and finally the total pressure generated by the planar source is the

summation of the pressure generated by each subelement

. N
J E,x+fiy+G,
¢@JyJ) —; 2:

i=1 j/E22+FZ.2
2 2 2
a-—b.+c. _. /_—_2 2 2 .
475—113 3k 0 +2 +hi—e—sz

x 2 do , 5.4
[12412—02 02+h;2 ( )

Ci

 

where Ez-x + Fiy + G,- = 0 describes the line that is coincident with each side of
the planar source, ai, bi, and c,- are the parameters defined by Figure 5.1, and h,- =
Ell-3+1? y+G
i¢éth
edge of the trianguar source. The sign of h,- is defined such that the superposition

z is the distance f1 om the projection of the observation point D to the

of subelements correctly respresents the total potential generated by the flat planar
source. The parameter N describes the number of sides/vertexes in the flat planar

source where N = 3 for a triangular source and N = 4 a rectangular source.

91

 

 

 

 

 

 

 

b)

Figure 5.1. Parameters defined for F N M calculations with a triangular source AABC.
The lengths ai, bi, and c,- are defined for three different triangles that share the vertex
D in each subfigure. a). The projection of the observation point D is located inside
of the triangular source. b). The projection of the observation point D is located
outside of the triangular source.

92

 

 

 

 

a) Subtetrahedron. b) Subpyramid.

Figure 5.2. Subdividing a prism about an observation point to form subtetrahedrons
and subpyramids. (a). Subtetrahedron. (b). Subpyramid.

5.1.3 FNM calculations for a volume source

Three common examples of volume elements are tetrahedrons, bricks, and prisms.
In order to calculate the volume potential integral, the volume is subdivided with
respect to the observation point. For the volumes considered here, the domain can be
subdivided into subtetrahedrons and /or subpyramids according to whether the face
shape is a triangle or rectangle [31] as shown in Figure 5.2. Thus, efficient formulas

for the subtetrahedrons and subpyramids are obtained.

FNM for subvolumes

In order to utilize the F NM expression for a planar source, the subvolumes are sub-
divided into layers of triangles and quadrilaterals. The procedure for calculating
potential integrals for subvolumes is given by Khayat and Wilton [31] and is adapted

here for the fast nearfield method.

1. Subdivide the volume element about the observation point into subdomains

including subtetrahedrons and subpyramids.

2. For each subtetrahedron, choose the number of sample points along the height

93

4— Observation Point

        

--
—-
--——--
-’
-—

—‘—

   

 
 

J-

 
 

j--—

Figure 5.3. The geometric configuration indicating how the potential integral is eval—
uated for a subtetrahedron.

direction of the subtetrahedron. The cross—section at each sample point is a

planar element.

3. Find the vertices and the projection of the observation point for each planar

element.

4. Calculate each potential integral over the planar element using the fast nearfield

method (FNM).

5. Transform the observation point from the global coordinate system to the local

coordinate system.

94

For potential calculations within subdomains, a local coordinate system is used as
shown in Figure 5.3 where the projection of the observation point onto the base
plane is the origin of the local coordinate system and the base plane is located in the
z = 0 plane of the local coordinate system. Let (x, y, z) denote the coordinates of
the observation point in the local coordinate system. The potential integration over

a subdomain V is given by

e—ij hv e—ij I
= d -
<I>(,z:cy, )=///V 47rR dV f0 //3 41rR dsz, (55)

where hv is the distance from the observation point to the base of the subtetrahedron

 

 

(i.e., hv is the height of the subtetrahedron). In the local coordinate system, hv = z,
where z is the local coordinate of the observation point. The inner integration over the
source plane 3 represents the potential produced by a triangular source distribution,
which is the shape of the subtetrahedronal cross-section. Figure 5.3 illustrates the
geometrical configuration of the potential calculations over the subtetrahedron. The
base plane is denoted as sb, where b indicates the base plane, and the plane with a
triangular cross-section is Sc at z = z’, where c represents the cross plane. The result

of the inner integral is given by [12]

 

2 _ I2 2_ _
e JkR or? e Jk\/0 +(z z) +hCZ —e-3k(3_zl)
.//3 d8: 21h Ci/l; do,
47TH 47rk 0'2 + h;

(5.6)

 

 

where the subscript c indicates that the parameters are all calculated in the sc plane,
h,- is the distance between the projection of the observation point on the sc plane and
each side of the source in the 3c plane, and lil and l2? are the corresponding limits

of integration. Since these parameters are related to the value of z’, the relationship

/
EQ_ _ l_'czl _l_ci2 : Z " Z (57)

hm: 11721 [M2 2

95

is defined using the similar triangles in Figure 5.3. In Eq. 5.7, the subscript b indicates
the parameters that are all calculated in the 3b plane, and the subscript c indicates
the parameters that are calculated in the sc plane. Then, all terms containing the
subscript c are substituted for the terms containing the subscript b according to Eq.
(5.7), and then a change of variable according to the formula 0 = (z — z’)/z * Q is

performed, yielding

 

 

 

I
_-—__z 2 2 2 .
c _Jde [(222 e _Z_— \/Q +2 +hbi —e"3k(z—Z,)d
[[9 8 47k 2" “/1 Q'

(5.8)

Eq. (5.8) readily computes the potential generated by the cross-section Sc when the
position of the cross-section z’ is known. The parameters hi1 [2'11 and 1,-2 are computed
in the base plane and are the same for all cross-sections. Thus, the complexity
compared to Eq. (5.6) is greatly reduced. After substituting Eq. (5.8) into Eq. (5.5)

and exchanging the order of integration, Eq. (5.5) becomes

 

I
_, .z—z 2 2 2 ,

(5.9)

 

The inner integral (with respect to z’) is integrated analytically, so a single integral

for the potential is obtained, which is described by

lb22
2477/92 22% bl./lb1-1

<I>(:1:

 

 

 

_'k 2 22 h2. ’.
2(8 3 \/Q + + bl—1)—(e‘JkZ—l)\/Q2+z2+hgi
2 2 2
(Q +112 b.)\/Q +22 +11%

 

 

(162. (5.10)

96

Eq. (5.10) and Eq. (5.5) are analytically equivalent, but the numerical performance
of each expression differs. Eq. (5.5) is a triple integral that includes a 1 / R singularity.
In Eq. (5.5), the total number of sample points is the multiplication of the number
of sample points in each direction, which is much larger than the number of sample

points required to evaluate the single integral in Eq. (5.10).

Volume FNM for the prism

To calculate the potential generated by a prism, the prism is subdivided about the
observation point to form five subvolumes including three subpyramids and two sub—
tetrahedrons. The potential of the prism is the superposition of the contributions
from the five subdomains. Let Amrz: + Bmy + sz + Dm = 0 denote the general
form of the plane that passes through each face of the prism, where <I>m is the result of

the potential integral over each subdomain, then the overall potential over the prism

(pprz'sm is given by

(I) . _ Z5: (Amx+Bmy+sz+Dm)

 

5m. (5.11)

The term (Ama:+Bmy+sz+Dm)/|Am$+Bmy+sz+Dm| is the sign function,
which determines whether the contribution from a subdomain is positive or negative.
The values of Am, Bm, Cm and Dm are selected such that the contributions from

all the subdomains for an observation point inside the prism are all positive.

5.1.4 Error calculations

For potential field calculations, the numerical error 77(x, y, z) is defined as the nor-

malized difference between the reference field and the computed field according to

lq)(xi y: Z) _ (prefab 3’: Z“
ma$|¢reff$1y1zll

 

n(r,y12) = (5.12)

97

 

Table 5.1. The geometry of the prism.

 

Vertex 1 2 3 4 5 6
(x, y, z) [m] (0,1,0) (0,0,0) (1,0,0) (0,1,1) (0,0,1) (1,0,1)

 

 

 

 

 

 

 

 

 

 

 

where (pref-(cc, y, z) is the reference potential and <I>(:1:,y, z) is the computed pressure

field. The maximum error 77mm: is defined as

77771011? = 1193):, 7](.’L‘, y) 3)) (513)

3 1

and the number of significant digits n is given as

5.2 Results

5.2.1 Comparisons of potential evaluated at single points

The shape of the prism is shown in Figure 5.4, and the coordinates of each vertex are
given in Table 5.1. This is the same prism used in thesingularity cancellation method
[31]. The wavelength is 10m. The observation points are located at (x, y, z) = (1 / 3m,
1 / 3m, d) where d = 0.5, 1.0, and 1.25 [m]. The reference potential given in Table
5.2 is calculated using the singularity cancellation method with 30, 30, and 30 Gauss
abscissas in the :12, y, and 2: directions, respectively. The reference potential is accurate
to 15 digits. The input parameters and the values for the reference potential are given
in Table 5.2.

The number of significant digits in each result is plotted as a function of the
number of sample points for the F NM and the singularity cancellation method in
Figure 5.5, which shows that the uniformly excited 1D FN M expression for the volume

integral converges much faster than the singularity cancellation method. For the point

98

 

 

 

 

Figure 5.4. The prism geometry, where vertex B is coincident with the origin, and
vertices C, A, and E are located on the :13, y, and z axes, respectively.

(1 / 3m, 1/3m, 0.5m), the the uniformly excited 1D F N M expression for the volume
integral uses 9, 9, 6, and 9 times fewer sample points than the singularity cancellation
method for answers accurate to 2, 3, 4, and 5 significant digits, respectively. For
the point (1 /3m, 1/3m, 1.0m), the the uniformly excited 1D FNM expression for
the volume integral uses 6, 12, 27, and 15.75 times fewer sample points than the
singularity cancellation method for answers accurate to 2, 3, 4, and 5 significant
digits, respectively. For the point (1 /3m, 1/3m, 1.25m), the uniformly excited 1D

F NM expression for the volume integral uses 6, 5, 12, and 10.125 times fewer sample

99

Table 5.2. Reference potential fields are evaluated at three points (:17, y, z) =
(1 /3m, 1/3m, d), where d = 0.5m, 1.0m and 1.25m. The reference fields are com-
puted using the singularity cancellation method with 30, 30, and 30 abscissas in the
:12, y, and 2 directions, respectively.

 

Abscissas Observation point Rea1{<I>} Imaginary{<l>}

 

x y z x y z
30 30 30 1/3 1/3 0.5 01101990007812 -0.0246818749055

30 30 30 1 / 3 1/ 3 1.0 0.07701456144055 —0.02427846658870
30 30 30 1/ 3 1/ 3 1.5 0.04833922443213 -0.02377978649190

 

 

 

 

 

 

 

 

 

 

 

 

 

 

points than the singularity cancellation method for answers accurate to 2, 3, 4, and
5 significant digits, respectively.

Next, the potentials are compared for values computed with the FNM and the
singularity cancellation method where the observation points have small 2 values.
The observation points are located at (:13, y, z) = (1/3m, 1/3m, d), where d =
0.1m, 0.01m, and 0.0001m. The reference potentials at the three observation points
are given in Table 5.2.

The number of significant digits are plotted as a function of the number of sample
points for the FN M and the singularity cancellation method in Figure 5.6. The FN M
demonstrates similar performance to that demonstrated in Figure 5.5 for larger 2
values in terms of convergence for all three observation points with small 2 values. In
contrast, the singularity cancellation method converges much more slowly for smaller
2 values such as 2 :2 0.01m and z = 0.0001m. A comparison between Figure 5.6
and Figure 5.5 indicates that, as the z coordinate approaches zero, the uniformly
excited 1D F NM expression for the volume integral and the singularity cancellation
method converge more slowly. Table 5.5 summarizes the ratio of the number of
sample points required to achieve a given accuracy with the uniformly excited 1D
F NM expression for the volume integral and the singularity cancellation method for

three points (1:, y, z) = (1/3m, 1/3m, d), where d = 0.1m, 0.01m, and 0.0001m.

100

 

_LA—L.‘
ON-hm

 

 

M.- 2 (1/31/3-0.5)(m]-
.. i§. . (1I31I31nm1 ..
1 - (mm-251m

number of significant digits

 

 

 

 

 

 

 

acre-ham
C
II

0 100 1,000 10,000 100,000
number of sample points

a) FNM

 

,6_,,. 2:22:53? f
. (1l31/30.5)[m] 5 ;;332235
14" I ‘ (1’31/31llm] . ;;';I;'.-'-,i-- 2 2: EEE
12» - (1/31/31.25)[m] - 5 2222525
10-, 1355:3; 3,2 2,232.25: " "

 

    

 

 

 

  

number of significant digits
on

 

 

 

O 100 1,000 10,000 100,000
number of sample points
b) Singularity cancellation method

Figure 5.5. The number of significant digits in the computed potential for a prism
plotted as a function of the number of sample points for three observation points
(:13, y, z) : (1/3111, 1/3m, d), where d = 0.5m, 1.0m, and 1.25m.

101

Table 5.3. Simulation parameters that achieve between 2 to 5 significant digits with
the F NM and the singularity cancellation method for the potential evaluated on the
observation points (:13, y, z) 2 (1/3, 1/3, d), where d 2 0.5, 1.0, and 1.25 [m].
Parameters listed include the number of sample points for each observation point
and the ratios of the number of sample points required to achieve a specific accuracy
relative to the number required with the FNM. (a) Significant digits 2 and 3 and (b)
Significant digits 4 and 5.

(a)

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Observations Significant digits 2 3
F N M Sigularity F NM Sigularity
cancellation cancellation
d 20.5 Sample points 36 324 54 324
Sample points 1x 9 x 1x 6 x
relative to F NM
(1 21.0 Sample points 36 216 36 432
Sample points 1x 6x 1 x 12 x
relative to FNM
d 21.25 Sample points 54 324 72 360
Sample points 1 X 6 x 1 x 5 x
relative to FNM
0))
Observations Significant digits 4 5
F NM Sigularity F NM Sigularity
cancellation cancellation
d 20.5 Sample points 90 54 108 972
Sample points 1x 6x 1x 9x
relative to FNM
d 21.0 Sample points 36 972 72 1134
Sample points 1x 27x 1x 15.75x
relative to FNM
d 21.25 Sample points 108 1296 144 1458
Sample points 1x 12 x 1x 10.125 x
relative to FNM

 

 

 

 

 

 

 

Table 5.4. Reference potentials evaluated at three points (:r, y, z) 2
(1/3m, 1/3m, d), where d 2 0.1m, 0.01m, and 0.0001m. The reference potentials
are computed using the singularity cancellation method with 30, 30, and 30 abscissas
in the 2:, y and 2: directions, respectively.

 

Abscissas Observation point Re{<I>} Imag{<I>}

 

x y z x y z
30 30 30 1 / 3 1 / 3 0.1 -0.09045651713588 0.02442323471174
30 30 3O 1 / 3 1 / 3 0.01 -0.07858100147339 0.02429436586981
30 30 30 1/ 3 1 / 3 0.0001 -0.07703049150655 0.02427862714031

 

 

 

 

 

 

 

 

 

 

 

 

 

 

For the point (1 /3m, 1/3m, 0.1m), the the uniformly excited 1D F NM expression
for the volume integral uses 6, 7.7, 6.3, and 7.38 times fewer sample points than the
singularity cancellation method to achieve 2, 3, 4, and 5 significant digits, respectively.
For the point (1 / 3m, 1/3m,0.01m), the uniformly excited 1D FNM expression for
the volume integral uses 8, 10.8, 10.125, and 17.14 times fewer sample points than the
singularity cancellation method to achieve 2, 3, 4, and 5 significant digits, respectively.
For the point (1 /3m, 1/3m, 0.0001m), the uniformly excited 1D F NM expression
for the volume integral uses 12, 27, 25, and 10 times fewer sample points than the

singularity cancellation method to achieve 2, 3, 4, and 5 significant digits, respectively.

5.2.2 The potential evaluated on a 3D grid

In this section, the potential is evaluated on a 3D grid. The a: and y directions extend
from -0.5m to 1.5m with an interval of 0.02m, and the z direction extends from 0
to 1.5m with an interval of 0.03111. The reference is computed using the singularity
cancellation method evaluated with 1000 Gauss abscissas in each dimension. Figure
5.7 shows the number of sample points and the number of significant digits for the
uniformly excited 1D FNM expression for the volume integral and the singularity
cancellation method evaluated on a 3D grid. Table 5.6 summarizes the ratio of the

number of sample points required to achieve a given accuracy with the uniformly

103

 

 

   

 

.53“ . (1I31l30.01)[m]

3:25; ; g; - (1/31/3o.ooo1)[m]g'§‘

number of significant digits

 

 

 

 

 

 

8

6
4-.. , I (1/31/3 0.1)[m]

2

o

1

0 100 1,000 10,000 100,000
number of sample points

a) FNM.

 

15-.. ,, ; 5,335.33. ,_
. (1/31/30.1)[m]

A (1/31/30.01)[m]
- (1/31/3 0.0001)[m]

 

 

 

 

10*

number of significant digits

 

 

 

..

0 I, W553; 5333352:
10 100 1,000 10,000 100,000
number of sample points

 

b) Singularity cancellation method.

Figure 5.6. The number of significant digits in the calculated potential achieved for
the prism in Figure 5.4 plotted as a function of the number of sample points for
three observation points (3;, y, z) 2 (1/3m, 1/3m, d), where d 2 0.1m, 0.01m, and
0.0001m.

104

 

Table 5.5. Simulation parameters that achieve between 2 to 5 significant digits with
the uniformly excited 1D FNM expression for the volume integral and the singularity
cancellation method for the potential evaluated at the observation points (:13, y, z) 2
(1/3m, 1/3m, d), where d 2 0.1m, 0.01m, and 0.0001m. Parameters listed include
the number of sample points required at each observation point and the number of
sample points required to achieve a specified error relative to the F NM. (a) Significant
digits 2 and 3. (b) Significant digits 4 and 5.

(a)

 

 

 

 

 

T-L

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Observations Significant digits 2 3
F NM Sigularity FNM Sigularity
cancellation cancellation
d 20.1 Sample points 72 432 126 972
Sample points 1x 6 x 1x 7.7x
relative to FN M
d 20.01 Sample points 54 432 90 972
Sample points 1x 8 x 1 x 10.8 x
relative to FNM
d 2 0.0001 Sample points 36 432 36 972
Sample points 1x 12 x 1x 27x
relative to FNM
0))
Observations Significant digits 4 5
F NM Sigularity FNM Sigularity
cancellation cancellation
d 20.1 Sample points 180 1134 234 1728
Sample points 1x 6.3x 1x 7.38x
relative to FNM
d 20.01 Sample points 144 1458 126 2160
Sample points 1 x 10.125 x 1x 17.14 x
relative to FNM
d 2 0.0001 Sample points 72 1080 108 1080
Sample points 1x 15x 1x 10x
relative to FNM
105

_\
O5

 

 

0 FNM

  

i
_L
A

 

 

ig ts

 

...........................................................

_xa
ON

number of significant d

 

 

 

no»
i

o 31100 3,000 136,000

number of sample points

Figure 5.7. The number of significant digits achieved in calculations of the potential
of the prism shown in Figure 5.4 plotted as a function of the number of sample points
evaluated over a 3D grid. The results show that the FNM is accurate to a large
number of significant digits than the singularity cancellation method for the same
number of sample points.

excited 1D FNM expression for the volume integral and the singularity cancellation
method on the 3D computational grid. For this computational grid, the uniformly
excited 1D FNM expression for the volume integral uses 3, 5, 7.2, 9.33 times fewer
sample points less than the singularity cancellation method for results accurate to 2,

3, 4, and 5 significant digits, respectively.

106

 

Table 5.6. Simulation parameters that achieve between 2 and 5 significant digits with
the uniformly excited 1D F NM expression for the volume integral and the singularity
cancellation method for the potential evaluated on a 3D grid. Parameters listed
include the number of sample points for each observation point and the number of
sample points required to achieve a specified error relative to the uniformly excited 1D
FNM expression for the volume integral. (a) Significant digits 2 and 3. (b) Significant
digits 4 and 5.

(a)

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Significant digits 2 3
FNM Singularity F N M Singularity
cancellation cancellation
Sample points 72 216 108 540
Sample points relative to the F NM 1x 3x 1 x 5x
(b)
Significant digits 4 5
F N M Singularity FN M Singularity
cancellation cancellation
Sample points 180 1296 324 3024
Sample points relative to the F N M 1x 7 .2 x 1 x 9.33x

 

 

 

5.3 Discussion

5.3. 1 Other geometries

The potential integral in Eq. (5.1) is a triple integral, and the uniformly excited
1D FN M expression for the volume integral reduces this integral to a single integral
when the excitation is uniform. The single integral needs far fewer sample points
to converge than the triple integral that is evaluated for the singularity cancellation
method. Thus, the uniformly excited 1D F N M expression for the volume integral is an
ideal method for calculating the reference potential. The fast approach demonstrated
here can be easily extended to any volume source with polygonal faces. The potential
integrals over those domains are the sum of the contributions from each side of the

polygonal sources or each face of the volume sources. As seen in Eq. (5.4), the number

107

of terms required for a triangular source is 3 and for a rectangular source is 4.

5.3.2 Sample point calculations

Both the singularity cancellation method and the uniformly excited 1D FNM expres-
sion for the volume integral subdivide the prism into five subdomains including two
subtetrahedrons and three subpyramids. Both methods evaluate three integrals for
each subtetrahedron and four integrals for each subpyramids. In sum, both methods
evaluate 18 integrals at each point for the prism shown in Figure 5.4. Each integral
in the singularity cancellation method is calculated with the same number of sample
points. If 113;, rig, and 11.2 denote the number of Gauss abscissas used for each integral
in the singularity cancellation method, then the number of sample points used for
each integral is 113 2 Tll'nynz. Thus, the total number of sample points for the sin—
gularity cancellation method is 18113. For the uniformly excited 1D FNM expression
for the volume integral, each integral is a single integral, and the lower and upper
limits are given by lbz'l and lbz'2 in Eq. (5.10). In each calculation, the integral is
further subdivided into the summation of two integrals if lbz’l and lbz'2 have opposite
signs, where the limits of the subintegrals are from 0 to [(712 and from lbz'l to 0. When
lbz’l and lb172 have the same sign, no subdivisions are needed. Thus, if the number of
abscissas for each integral in the uniformly excited 1D FNM expression for the volume
integral is n, then the total number of sample points for the uniformly excited 1D
F N M expression for the volume integral calculations performed at a point due to a
prism volume source is between 1811 and 3611.. Thus, the FN M calculations is C(11),

and the singularity cancellation method is 0(713) for calculations at a single point.

5.3.3 Future work

The FNM algorithm is capable of calculating the uniformly excited potential integral.

However, the FNM formulas have not yet been derived for non-uniform potential

108

sources. This difficulty can in part addressed by applying the polynomial apodization
from Chapter 4 to the potential integrals. Thus, certain apodized potential integrals
can be solved. Meanwhile, a large variety of functions can be approximated very well

using polynomials.

5.4 Conclusion

A fast 1D integral for calculating uniformly excited 3D potential integrals is derived
and evaluated. The 1D integral is based on the FN M expressions that were derived
for acoustic radiation problems [8, 12]. The FNM subdivides the prism volume source
into five subdomains, and on each subdomain, the FNM eliminates the singularities
and reduces each triple integral into a single integral. The FNM expressions can be
extended to different volume sources which consist of polygons on each face. The
results are compared at single points and on a large three-dimensional grid. For the
six observation points evaluated here, the FNM reduces the number of sample points
by a factor of 5 to 27 relative to the singularity cancellation method, where the number
of significant digits ranges in each result from 2 to 5. The FN M reduces the number of
sample points by a factor of 3 to 9.33 relative to the singularity cancellation method
when evaluated on a 3D grid and the result is accurate between 2 and 5 significant
digits. Thus, the F N M is an ideal method for calculating uniformly excited potential

integrals.

109

CHAPTER 6

A Fast Nearfield Method for
Numerical Evaluation of Surface
Integrals with Polynomial

Apodization

Numerical calculation of surface integrals involved in integral equations are very im-
portant in scattering problems. Polynomial functions which are popularly used to
approximate the electric and magnetic currents. For example, the most common
used Rao-Wilton-Glisson basis function is also a linear basis function [24]. Potential
integrals over a surface domain are often singular and direct evaluation of poten-
tial integrals may encounter numerical difficulties. Singularity subtraction methods
[23, 25, 26, 27, 28] or singularity cancellation methods [30, 31, 32, 33] are often used
to acheive better performance. Typically, those methods don’t diminish the number
of integraion of the potential integrals but manipulate the integrands to eliminate sin-
gularities instead. This approach is reasonable for general surface integrals, however,

more efficient expressions can be acheived when polynomial apodized surface integrals

110

 

are considered. By ultilizing the fast nearfield method (FNIVI) [8, 12] for planar source
that is uniformly exicited, which is proved to achieve much better accuacy than the
Rayleigh-Sommerfeld integral, the 1D F NM expressions for the surface integrals with
polynomial apodization is obtained.

To improve the numerical performance of surface integral calculations with the
polynomial apodization, a single integral FNM expression is derived based on the
FNM expressions for a planar source that is uniformly excited. The FNM expres-
sions for planar sources with linear, quadratic, and cubic apodization functions are is
obtained by utilizing the two derivatives. After the 1D FNM integrals are obtained,
the 1D FNM integrals and the Rayleigh-Sommerfeld integral are both evaluated on a
2D grid. The results indicate that the 1D FNM expressions requires fewer abscissas

than the Rayleigh-Sommerfeld integral to achieves a given accuracy.

6.1 Method

The potential generated by a planar source that is uniformly excited is given by

e—ij
<I)(,:1:y,z )/2/e 47rR ds (6.1)

where R 2 \/ (:17, — l‘)2 + (y, - y)2 + (z — z’)2 is the distance between the observa-

 

 

tion point (1:,y,z) and the point (:1:’,y’,z’) on the piston face and s is the area of
the planar source. The potential in Eq. (6.1) is normalized with respect to —v08jwt,
where 110 is the constant normal particle velocity evaluated on the source of the planar
source, a) is the excitation frequency in radians per second. Two derivatives can be
evaluated with respect to R,

0R _ (.1." -:1:) 81?. (y, — y)

___— d——2 , .2
81’ R an 63/ R (6)

 

Thus the following two equations are easily obtained,

111

 

 

 

(.17, __ 3:)6—ij _ 188—.ij d (y! _ y)e—ij _ j 68—.ij

R ‘1. 6:1." an R' 7% 611’ (63)

The potential in Eq. (6.1) is rewritten as

ff: (ad—M (1%) e—ijd ’d" <64)
(Ham 2): 47er 32]: L1(y’) 471R :1: y '

. . . _ B yLC
where b 18 the her ht of the trian ular or the rectan ular source, L I 2
g g g 21 1 4.2—2

 

 

and L1(y )=_BJ%’__1 for a triangle which is formed by the three lines y— —— 0,
Alx+Bly+Cl 20 and A22+Bgy+02— — 0, and L2(y ’)- — a and L1(y ’) 2 0 for

a rectangle that has width a and height b.

6.1.1 Uniformly excited source

The F NM for the uniformly excited planar source is equivalent to the Rayleigh-

Sommerfeld integral given by the following expression,
—ijds
F N M
[[96 47er

.2 .2— r—
+ b2 .- 2 2 2 .
_Z_2_1_Z_e—_]k‘ a +2 H1 —e_9k3
:41er gag/2b? 62 022.12 do , (6.5)
__Z__2J.__Z_ Z

Ci

 

 

where 5,1,- 2 (Bic): + Fig + G,)/‘ / E? + F 2'2 combines the Sign and the height terms
within a single expression, N 2 3 for the triangular source and N 2 4 for the
rectangular source. Calculations with Eq. 6.5 compute the values of C2“, E2“: F11 and
C,- only once for each edge of the planar source, whereas the values of a,- and b,- are
calculated once for each (:15, y) pair. Unlike the expressions for the impulse response
that change depending on the spatial coordinate, Eq. 6.5 is a general formula that

computes the nearfield pressure generated by any triangular source. The distance R

112

 

 

is given by R 2 \/(:r' — 11:)2 + (y’ - y)2 + 22, where (1.", y’, 0) and (:13, y, 2) represent
the source point and the observation point, respectively. In the following derivation of
surface integrals with polynomial apodization, the following FNM expression is also

needed.

 

 

FNA/[2— —
f/se‘1 =247rk 2§S,-l,

x
a2—b2C: c2

2‘32
-'k (32 2 [2 .
e '7 +0 +1(jk‘/z2+02+122+1)—e-~7kz(jkz+1)

02+lz‘2

 

dd, (6.6)

 

where N 2 3 for the triangular source and N 2 4 for the rectangular source.

6.1.2 Linear apodization

Let f (:13) represent an apodization function defined in the :1: direction and f (12’ ) 2 11:, .

I

First, f (1:, ) 2 :r’ is written as :1: 2 (1r, — 2:)+:1:. The potential for the linear apodized

triangular source is the sum of the following two integrals,

64ij I I
—ij —ij
I e C I I
2 — 2d ' . .7
//S(:c 2:) 47rR dr:+//S:1: 471R drrdy (6)

The first of the two double integrals on the right hand side of Eq. (6.7) is reduced to

 

 

 

a single integral by

L2(y')

. b .
——]kR —ij
I 8 _ /_ 6 1,1
[[805 :r) 47rR ds—/ (:1: :1:) 471R dxdy
0 L1(y’)

 

 

113

 

is given by R = \/(£L, — :1:)2 + (y’ -— y)2 + 2:2, where (:c’,y’, 0) and (:1:, y, 2) represent
the source point and the observation point, respectively. In the following derivation of

surface integrals with polynomial apodization, the following FNM expression is also

needed.

02+ 624,2

FNMZ2/fseijCL9=4 771111.21 Silt/2“,)C ficzx

Z

_‘k / 2 2 12 .
e ‘7 2 +0 +1 (jk(/22+02+122+1) —e_9kz(jkz+1)

2 2
0 +12

 

 

dd, (6.6)

where N 2 3 for the triangular source and N 2 4 for the rectangular source.

6. 1.2 Linear apodization

Let f (:10) represent an apodization function defined in the :1: direction and f (:1:, ) 2 :1:, .

I

First, f (x, ) 2 :1:, is written as :r 2 (:1:, - x) + :1:. The potential for the linear apodized

triangular source is the sum of the following two integrals,

I [J _// I —
3:,y,2)— a:
s
—ij —ij
I 8 6 I I
2 r— :1- . .
[[901 :1:) 47rR d11y+//Sa: 471R dasdy (67)

The first of the two double integrals on the right hand side of Eq. (6.7) is reduced to

 

.7
R dx’dy’

 

 

a single integral by

L I)
e—ij b 2(9) e—th 1 1
/ (mi—:1:) 471R d8=/ (:1: —:1:)e4 Rdasdy
s
0

L1(y’)

 

 

113

 

b L2('.III) . +1,
' fl“? - . I: , I
= _L/ gig—T—dedyI2 L/e—Jlexl L2(y,)dyI
0 L1(y’) 0

 

 

 

+b
__J'_/( e-J'k\/(L2(y’)-:r)2+(y’—y)2+22_ —jk\/(L1(y’)-rc)2+(y’—y)2+z2)d.
k e y.
0
(6.8)
Denote Ml(u )2£f0u+b e—JkR|:,: _iQE: ,gdyI, thus
e—JkR
f/Scc’ — :1:) ds = Ml(1). (6.9)

The second of the two double integrals on the right hand side of Eq. (6.7) is also

reduced to a single integral easily by

 

e—ij
d 2 FNM .
[[813 47rR S :1: (610)

Eqs. (6.9) and (6.10) are added together and, the pressure for the linear apodized

triangular source is given by

(1)1(23I,yI,zI) -_- Ml(1) + xFNM. (6.11)

6.1.3 Quadratic apodization

Let f (:1:) represent a quadratic apodization function defined in the 113 direction and
f (:rI ) 2 27,2. First, f (xI ) 2 (112’)2 is written as (33’)2 2 :cI(:rI — 2:) +£L'1'I . The pressure

for the quadratic apodized triangular source is the sum of the following two integrals,
@( I26 I
2( :rI ,I,y zI) ff (I) —————dx dyI
R I I
(12: dy . (6.12)
l

—ij —-J
_ ./ I_ e _I ,I I6
2//Sz(:1: :r)47TRd:Ldy+//S$334R

 

 

114

The first double integral on the right hand side of Eq. (6.12) is reduced to a single

integral by

 

b L2(y’) .1 b L2(y’) .
/ xI(a:I — x)e—]IRdedyI = _j_ / xI—ae—IdexIdyI
47rR 47rk 823’
0 L1(yI) 0 L1(yI)
_ +6 L( 6 L2(?/)
=_L ’e-ij‘CIr‘ZIy —ij I I
47rk/x :1:/2(yL1dyI_47rk/ 8 dx dy (6.13)
0 0 L1(y’)

Denote A12(u)2 ‘LEfO [iffy/I; ue’ijdedyI, thus
I
b L2(y ) _ij
/ / xI(:I:I — ”634m dedyI = Ml(a:I) — 1142(1) = M1(:rI) — FNM2. (6.14)
0 L1(y’)

The second of the two double integrals on the right hand side of Eq. (6.12) is reduced

 

to a single integral by applying the the result from Eq. (6.11)

i
0 L1
Eqs. (6.14) and (6.15) are added together, and the pressure for the quadratic apodized

L2(y’) k
I, xxI— RRdx’dy’= 23Ml(1)+a:2FNM. (6.15)
(y’)

triangular source is given by
(D2(:1:I,yI, zI) = Ml(:1:I) — FNM2 + 22611 (1) + xQFNM. (6.16)

6.1.4 Cubic apodization

Let f (:1:) represent an cubic apodization function defined in the a; direction and

f(:z:I) 2 (:L'I)3. First, f(:L'I) 2 (2’)3 is written as (:rI)3 2 xI2(:rI — :1:) + a:(a:I)2.

115

 

The pressure for the cubic apodized triangular source is the sum of the following two

 

 

integrals,
— des
3I(:r ,yI, zI) 2//S(2:I) ,)3e J
47TR
b L2(yI) b L2('yI)
— / (xI)2(:1:I — :1:)e—ijded I + :1:] (xI)2£fi-R:d$Id I (6 17)
_ 47rR III 47rR y I I
0 L1(y’) O L1(y’)

The first double integral on the right hand side of Eq. (6.17) is denoted as

 

2L —ijl‘I =L2(y I)
47rk/ (x I28) ll‘I 33I=L1(y I)

= Ml(:1:I2) — 21120!) (6.18)

where M2(:I: 2 47th f0 le L(2 (I)? 2’ e_ijd:z:IdyI and is rewritten as

b L2(yI) . b L2(yI)
AI2(I) 2 —j—k/ (.rI — 2:) e )ijdedyI-F +fi/ we __ijde dyI
0 0

L1(yI) L1(yI)

(1111(3) — 2.1/11(1)) + xii/12(1). (6.19)

PPM-lb

The second double integral on the right hand side of Eq. (6.17) is also reduced to a

single integral by

116

 

b L2(yI)
313/ xI2——RRId:1:Idy 223(Ml(:rI)—M2(1)+:1:M1(1)+x2FNM) (6.20)
0

L1('y’)
Combining Eqs. (6.18), (6.20) and (6.19), the pressure for the cubic apodized trian-

gular source is represented by

(p3(xI2 y” zI) = M1 (332)

—%(le(R) — 3:1:M2(1) + IMl (1r) +(1L‘2 - —)M1(1 )+ x3FNM). (6.21)

6.1.5 Error calculations

The numerical error 77(36, y, z; k) is defined as the normalized difference between the
reference potential <I>ref(:1:, y, z; k) and the computed potential <I>(:z:, y, z; 16) according

to

I‘I>(I,y,z; k) -‘I’Tef(=r,y,2; k)|

 

"Mani/,2; k) = (6.22)
ma$|¢ref($1ya 2; k)|
The maximum error Tlmax is defined as
717an — 31:11ng 77(33 ya Z)- (623)

6.2 Results

The results are evaluated for a triangular source located in the z 2 0 plane, where
the positions of the three peaks are (0,4A,0), (—2A,0,0) and (2A,0,0). The refer-

ence potential is computed for the linear, quadratic, and cubic apodization functions

117

 

Table 6.1. Simulation parameters that achieve 10%, 1% and 0.1% peak normalized
error in the potential obtained with the fast nearfield method and the Rayleigh-
Sommerfeld integral computed using linear apodization for a triangular source. Pa-
rameters listed include the number of Gauss abscissas, the computation time, and the
ratios of the computation times relative to the times obtained with the fast nearfield
method.

 

Linear Apodization

 

 

 

 

 

 

10% 1% 0.1%
F NM Rayleigh FNM Rayleigh FNM Rayleigh
Abscissas 12 42 x 42 14 99 x 99 15 190 x180
Time 0.14s 1.39s 0.16s 7.718 0.175 25.475
Computation 1 x 9.97 x 1 x 46.85 x 1 x 145.82 x
Time relative
to the FNM

 

 

 

 

 

 

 

 

 

 

 

f (1‘) described in the Section 6.1, namely f (:1:) 2 :1:, f(a:) 2 x2, and f(:1:) 2 11:3,
respectively. The potential is evaluated in the a: 2 1.0) plane, where the grid in
the a: 2 1.0A plane extends from -—2)\ to 6A in the y direction, and the grid spacing
in the y-direction is 0.08/\. The grid extends from 0.02A to 1.0A in the z direction
with a spacing of 0.02/\. The potential in the z direction is shifted by 0.02) relative
to the piston source in the z 2 0 plane for both the Rayleigh-Sommerfeld integral
and the fast nearfield method in order to avoid the most severe singularities in this
location for the Rayleigh-Sommerfeld integral. The reference potential is obtained

using the Rayleigh-Sommerfeld integral evaluated with 100,000 Gauss abscissas in

each direction, and the results are computed on a 50x101 spatial grid.

6.2.1 Linear apodization

The linear apodization function evaluated here is f (1:) 2 ac. The results for the fast
nearfield method and the Rayleigh-Sommerfeld integral are plotted in Figures 6.1

and 6.2. The numerical errors and the computation times are plotted in Figure 6.1 a)

118

Table 6.2. Simulation parameters that achieve 10%, 1%, and 0.1% peak normal-
ized error in the potential obtained with the fast nearfield nearfield method and the
Rayleigh-Sommerfeld integral computed using quadratic apodization for a triangular
source. Parameters listed include the number of Gauss abscissas, the computation
time, and the ratios of the computation times relative to the times obtained with the
fast nearfield method.

 

Quadratic Apodization

 

 

 

 

 

10% 1% 0.1%
FN M Rayleigh FNM Rayleigh F N M Rayleigh
Abscissas 12 33 X 33 14 94 x 94 15 190 x 180
Time 0.173 0.883 0.223 7.403 0.213 25.963

 

Computation 1x 5.25x 1 x 34.01 x 1x 123.18x
time relative
to the FNM

 

 

 

 

 

 

 

 

 

 

 

and Figure 6.1 b), respectively as a function of the number Gauss abscissas used in
one integration direction for the fast nearfield method and the Rayleigh-Sommerfeld
method. The time vs. error comparison is also plotted in Figure 6.2. Figures 6.1 and
6.2 show that the fast method converges much faster than the Rayleigh-Sommerfeld
integral for the linear apodization function evaluated here. For 10% peak normal-
ized error, the fast nearfield method is 9.97 times faster than Rayleigh-Sommerfeld
integral. For 1% peak normalized error, the fast nearfield method is 46.85 times
faster than Rayleigh-Sommerfeld integral. For 0.1% peak normalized error, the fast

nearfield method is 145.82 times faster than Rayleigh-Sommerfeld integral.

6.2.2 Quadratic apodization

The quadratic apodization function evaluated here is f (:1:) 2 $2. The results for the
fast nearfield method and the Rayleigh-Sommerfeld integral are plotted in Figures
6.3 and 6.4. The numerical errors and the computation times are plotted in Figure

6.3 a) and Figure 6.3 b), respectively, as a function of the number Gauss abscissas

119

h

 

 

 

 

_ §'-'-'Rayleigh

 

 

 

normalized error

 

 

 

 

o 2.0 46 610 810100120140160180200
number of abscissas

a) Maximum errors for linear apodization.

w
01

 

 

 

 

 

 

 

 

 

’u'? _FNM I I I I I I
cg 301... ._I,_,Rlayleligh ....... ....... ....... ....... ........ ”9!.
8 ' i ‘ ' i : i i
£25 . ........... _ ’3‘;
d) : : : : :l '
20 ............... 2 ...................... . ....... _ ....... ‘ ....... , ............. -

g éo""3
S 15 .............. .................... ....... My; .................... -
H ' .‘ i\l
1:310 ............... \'\’ ........................... .
Q Q 3 ; I"
g 5 ................ ....... «‘4‘ ................................... 1
° ; 2 .m?’ ‘

0.4—‘2 '. . 1 . ; g g

0 20 40 60 80 100 120 140 160 180 200

number of abscissas

b) Computation times for linear apodization.

Figure 6.1. a) Maximum errors and b) computation times obtained with the fast
nearfield method and the Rayleigh—Sommerfeld integral for a triangular source with
linear apodization.

120

 

 

,.l . . . Yt‘r‘rr]

 

normalized error

 

—FNM _
'-'-'Ray|eigh

 

 

-15 :sg. . .1
10 0.1 1 10 100
computation time (seconds)

 

 

 

 

 

Figure 6.2. Numerical errors plotted as a function of the computation time for a
triangular source with linear apodization. The results show that the fast nearfield
method achieves much better convergence performance than the Rayleigh—Sommerfeld
integral.

used in one integration direction for the fast method and the Rayleigh-Sommerfeld
integral. The time vs. error comparison is also plotted in Figure 6.4. Figures 6.3
and 6.4 show that the fast nearfield method converges much faster than the Rayleigh-
Sommerfeld integral for the quadratic apodization function evaluated here. For 10%
peak normalized error, the fast nearfied method is 5.25 times faster than Rayleigh-
Sommerfeld integral. For 1% peak normalized error, the fast nearfield method is 34.01
times faster than Rayleigh-Sommerfeld integral. For 0.1% peak normalized error, the

fast nearfield method is 123.18 times faster than Rayleigh-Sommerfeld integral.

121

 

 

normalized error

 

-15 i

l

I

 

 

i—FNM

 

 

'----Rayleigh

1o 1"... ...... ...... ...... ...... ...... ....... ...... .

 

A

Jwfi

. 1:»!qu

 

 

 

10 r

0204

to
0|

0 610 810100120140160130200

number of abscissas

a) Maximum errors for quadratic apodization.

 

I

I

 

O

 

—FNM ;
" '-'-'Rayleigh

 

 

computation time (seconds)
N N w
o a:
.,3

 

02040608

1

_......................................................... ......

: : 3
............................. {........}....".-
E i!
. f
j I.- ...... -
: : l"
E if 3
I O
2,. ;
.2"

.3
0|
1
5
i

_t
O
I
t

i

0 0|
3;
$

0 100 120 140160 180 200

number of abscissas

b) Computation times for the quadratic apodization.

 

Figure 6.3. a) Maximum errors and b) computation times obtained with the fast
nearfield method and the Rayleigh—Sommerfeld integral for a triangular source with

quadratic apodization.

122

 

normalized error

 

—FNM _
'---'Rayleigh

    

 

 

 

 

 

-15 “.2.

0.1 1 16 foo
computation time (seconds)

 

10

Figure 6.4. Numerical errors plotted as a function of the computation time for a
triangular source with quadratic apodization. The results show that the fast nearfield
method achieves much better convergence than the Rayleigh-Sommerfeld integral.

6.2.3 Cubic apodization

The cubic apodization function evaluated here is f (1r) 2 1:3. The results for the
fast nearfield method and the Rayleigh—Sommerfeld integral are plotted in Figures
6.5 and 6.6. The numerical errors and the computation times are plotted in F ig-
ure 6.5 a) and Figure 6.5 b), respectively, as a function of the number Gauss
abscissas used in one integration direction for the fast nearfield method and the
Rayleigh—Sommerfeld integral. The time vs. error comparison is also plotted in Fig—
ure 6.6. Figures 6.5 and 6.6 ShOW that the fast nearfield method converges much
faster than the Rayleigh-Sommerfeld integral for cubic apodization function evalu-

ataed here. For 10% peak normalized error, the fast nearfield method is 5.22 timas

123

Table 6.3. Simulation parameters that achieve 10%, 1% and 0.1% peak normalized
error in the potential obtained with the fast nearfield method and the Rayleigh—
Sommerfeld integral computed using cubic apodization for a triangular source. Pa-
rameters listed include the number of Gauss abscissas, the computation time, and the
ratios of the computation times relative to the times obtained with the fast nearfield
method.

 

Cubic Apodization

 

 

 

 

 

 

10% 1% 0.1%
F NM Rayleigh FN M Rayleigh F NM Rayleigh
Abscissas 12 33x33 13 ' 89x 89 15 190x 180
Time 0.403 2.073 0.433 14.943 0.503 55.253
Computation 1X 5.22x 1x 34.76x 1x 111.28x
time relative
to the FNM

 

 

 

 

 

 

 

 

 

 

 

faster than Rayleigh-Sommerfeld integral. For 1% peak normalized error, the fast
nearfield method is 34.76 times faster than Rayleigh-Sommerfeld integral. For 0.1%
peak normalized error, the fast nearfield method is 111.28 times faster than Rayleigh-

Sommerfeld integral.

6.3 Discussions

6.3.1 Other Polygonal Sources

The fast nearfield method for poynomial apodized triangular source and rectangular
source are derived in this paper. The derivation process can be generalized to other
polygonal sources. The F NM for the uniformly excited polygonal source is readily
computed using Eq. (6.5) and Eq. (6.6), and the number of 1D integrals evaluated is
equal to the number of sides of the polygonal source. The computation of Eq. (??)

also needs to be adapted according to the specific polygonal source.

124

 

 

 

1o .' 7 7 l l 7 f ' '

I I i j j :_FN
5 2 ? ""‘Rayleiyh
10 “i ..... ...... ...... ...... ..... ...... ....... ...... _

 

 

 

normalized error

 

-15 : a i I ; 1 '22
1o

 

 

 

o 20 4o 60 80100120140160130200
number of abscissas

a) Maximum errors for cubic apodization.

 

on
O

Y I T V T I I T T

I I '. I I i
_ FNM : E : l ’\

 

 

 

 

O
O

computation time (seconds)
N a
O O
i, i
5 ‘ 5
i is

 

 

if 1 L i L i 1 i
0 20 40 60 80 100120140160180 200
number of abscissas

b) Computation times for cubic apodization.

Figure 6.5. a) Maximum errors and b) computation times obtained with the fast
nearfield method and the Rayleigh-Sommerfeld integral for a triangular source with
cubic apodization.

125

 

normalized error

 

: 2 223325: 2 sséE '---'Ray|eigh
1 I 1::Zliii Z 2:..ZCIE : Z'.:::::
0 0.1 1 10 100

computation time (seconds)

 

 

 

 

 

 

 

Figure 6.6. Numerical errors plotted as a function of the computation time for a trian-
gular source with cubic apodization. The results show that the fast method achieves
much better convergence performance than the Rayleigh-Sommerfeld integral.

6.3.2 Higher order polynomials

In this paper, the FN M expressions for surface potential integrals with linear,
quadratic, and cubic polynomial apodization are derived. The F N M expressions of the
surface potential integrals with higher order polynomial apodization can be derived
using the same strategy. Many apodization functions, like 3in(), cos(), exp(), etc.,
can be approximated using higher order polynomials, thus potential integrals with
such apodizations can be approximated with higher order polynomials combined with
appropriate 1D FNM expressions. Since the FNM expressions converge very rapidly,
the accuracy of numerical evaluation of the potential integrals with general apodiza-

tion is mainly dependent on the accuracy of the approximation to the apodization

126

function.

6.4 Conclusion

Single integral expressions for calculating the potential generated by polynomial
apodized planar sources are derived and evaluated. Starting from the FNM for
the uniformly excited planar source, three polynomial apodizations including lin-
ear, quadratic, and cubic apodizations are reduced to a 1D integral. Higher order
polynomials can be obtained using the same strategy. These fast expressions remove
the singularities that appear in the Rayleigh-Sommerfeld integral and thus achieve
very fast convergence with a very small number of Gauss abscissas. Simulation results
are compared between the fast nearfield method and the Rayleigh-Sommerfeld inte
gral. For a given peak normalized error, the fast nearfield method always converges
much faster than the Rayleigh-Sommerfeld integral when polynomial apodization is

evaluated with the appropriate 1D FN M expression.

127

 

CHAPTER 7

A Fast Nearfield Method for
Numerical Evaluation of Volume
Integrals with Polynomial

Apodization

Numerical calculation of volume potential integrals involved in integral equations are
of great importance in scattering problems and the polynomial function is a very
popular function to approximate the electric and magnetic currents. For examples,
three-dimensional polynomials are adopted to approximate the entire-domain normal-
ized current density by Moraros and Popovic [23] when applied to optimize the volume
potential integrals involved in the moment-method analysis of 3D dielectric scatters.
Potential integrals over a volume domain are often singular and direct evaluation of
potential integrals may encounter numerical difficulties. Often, Singularity subtrac-
tion methods [23, 25, 26, 27, 28] or singularity cancellation methods [30, 31, 32, 33]
are used to improve the accuracy of the potential integrals. Those methods retain the

same number of dimensions over which the integration is performed. This approach

128

 

l‘T
l

is reasonable for general potential integrals, however, more efficient expressions can
be acheived when polynomial apodized volume potential integrals are considered.
Potentials generated by uniformly excited planar sources are computed with both
the Rayleigh-Sommerfeld integral [6] and the fast nearfield method (FNM) [8, 12].
The existing FNM achieves much higher accuracies by eliminating the 1 / R singularity
from the Rayleigh-Sommerfeld integral and by decreasing a 2D integral to a 1D inte-
gral. Numerical evaluation of potential integrals with polynomial apodization is an
important problem for scattering calculations, which motivates the derivatiaon of a
similar fast nearfield method to improve the numerical performance of those integrals.
A single F NM integral is derived based on the FNM exprassions for uniformly
excited and polynomial apodized polygonal sources in order to improve the perfor-
mance of potential integral calculations for volume sources with polynomial apodiza-
tion. The F NM expressions for uniformly excited sources are first reviewed, and then
the FNM expressions for the polynomial apodized planar source are obtained. The
volume source is then subdivided into five subdomains, and local coordinate system
is defined for calculations of the potential for a subdomain. The potential gener-
ated by the entire volume source is the superposition of the potentials evaluated over
each subdomain, where each contribution is represented by a single integral. Thus,
1D FNM expressions for a polynomial apodized volume source is obtained. The 1D
FN M expressions and the singularity cancellation method are evaluated on a 3D grid.
The results indicate that the FNM requires less time to achieve a given accuracy than

the singularity cancellation method for a polynomial apodized volume source.

7 .1 Method

The volume integrals evaluated in this article are given by

129

 

 

 

 

a) Subpyramid. b) Subtetrahedron.

Figure 7.1. Subdomains defined for a prism where the shared vertex is defined at the
observation point.

(7.1)

 

f/f'

where V is the volume integration domain, ./\(rI)47r represents the apodization over
the integration domain, R is the distance between the observation point (x, y, z) and
the source point (2’, y’, 25’), and k is the wavenumber. The potential in Eq. (7.1) is
normalized with respect to —vOejIUt, where '00 is the constant normal particle velocity
evaluated on the volume source, 62 is the excitation frequency in radians per second.
Here, the special case A(rI) 2 xI is considered.

The volume is first subdivided about the observation point into subdomains. Fig-
ure (7.1) illustrates a subpyramid and a subtetrahedron for a prism.

The volume integral shown in Eq. (7.1) is first rewritten as

z>=///.,w' v-—-.// . 2” 2...,va mi?

(7.2)

 

 

 

Then, two derivatives are evaluated with respect to R,

130

 

0R_(:rI-—a:) 8R_(yI—y).

 

 

 

 

 

_ _ __ _ 7.3
0:1:I R and By, R ( I
The following equations are then obtained,
(:1:I _ x)e_ij _ 188'.ij and (yI _ y)e-ij _ 1598—ij (7 4)
R _ k 8:13, R _ k 331’ E
The potential integral in Eq. (6.1) for a triangular base is rewritten as
3111— —C /B L (39/) —ij
<I>.(:1: y,z 2)‘//A(I) d3 2/ 1 1/ 2 A(TI)?———d.rIdyI
e47rR 0 L1(yI) 47rR
o L3(:1:I) e _ij —C2/A2 L4(:1:I) e—ij
= A(rI)———— dedyI + f / A(rI) dedyI,
[_01 /A1/() 47f R 0 0 47rR
(7.5)

Where L1(yI) = ('3in - C1)/A1,L2(1UI)= (-B2yI - C2l/A2, L3($I) = (—A1$I-
Cl)/Bl and L4(:rI) 2 (~A2xI — C2) / B2 for a triangular base which is formed by
the three lines 3; 2 0, A11: + Bly + Cl 2 0 and A22: + B2y + C2 2 0. The potential

integral in Eq. (6.1) for a rectangular base is rewritten as

—ij -ij
<I>,(:1: y,z) =//3A(I) )8 (13 2L: [_bbAO: )e RdedyI, (7.6)

where the center of the rectangular source is the origin of the coordinate system, and

 

 

the width and height of a rectangular source is 2a, and 2b, respectively.

7 .1.1 Uniformly excited planar source

The F N M for a planar source that is uniformly excited is equivalent to the Rayleigh-

Sommerfeld integral given by the following expression,

e—ij
VM 2
Fl [[9 471R d3

 

131

 

 

Figure 7.2. Geometric parameters defined for the potential calculations in a tetrahe-
dral subdomain.

a2 62_ 2
+ b- e_- /2 2 2 .
=4Tk 2i 2_b2c: C2 02+? , .

Ci

 

where 521,- 2 (E11: + Fz-y + at” E? + F 22 combines the sign and the height terms
within a single expression. The number of 1D integrals is given by N 2 3 for a
triangular source and N 2 4 for a rectangular source. Calculations with Eq. 7.7
compute the values of Ci, E11 F2" and C,- only once for each edge of the planr source,

whereas the values of az- and b,- are calculated once for each (2:, y) pair.

132

 

 

The fast nearfield method expressions for the uniformly excited potential integrals

derived previously for a prism, denote by (Pu, are given by

 

 

<I> _/// e—ijdV- % (Amx+Bmy+sz+Dm)q) (78)
U V 47TR _1 lAm$+Bmy+CmZ+Dm| ”In, I

where M 2 5 is the number of faces that the volume has, Amx+Bmy+sz+Dm 2 0
denotes the general form of the plane that passes through each face of the prism, and

q’um is the potential generated by each subdomain. The expressions for (bum are

___1 N
<I> 2 —— h -
um 47rk2 2.22:] m

 

__'k/ 2 2 ’22. _
162246 I Q +2 + m_1)_(e—sz_1)\[Q2+z2+h§,

[1611 (Q2 + ha.) Q2 + 2.2 + I232“

 

 

dQ, (7.9)

for each subvolume, where hv 2 z is the distance from the observation point to the

base of the subdomain, and the parameters hM1lb11’Ibi2 are shown in Figure (7.2).

7 .1.2 Linear apodization for a planar source
Linear apodization in the x direction

The potential integral for the apodization function f (:L'I ,yI ,zI ) 2 mI applied to a

planar source is given by

—ij —ij —ij
I8 I I I e I I e I I
d = — , . ——
//Sx 47rR dz y //S(:I: :1:) 47rR dzrdy+//Sx 47rR drcdy

. _ij
= LIT/E / 95537—116’111,’ +1rFNMu
S

 

 

m'form' (7'10)

For a triangular source, substituting Eq. (7.5) into Eq. (7.10) yields

133

 

 

 

le—ij I
//8x 471' R ——d3: dy 2:1:FNM uniform+

 

j f-CI/Bl —jk\/(L2(y’>—x>2+<y'—y)2+(z—z')2
47rk 0 e

 

—e_jk\/((((yL1,)(l—x)2+y “W Z)2 dy’ (7.11)

For a rectangular source, substituting Eq. (7.6) into Eq. (7.10) yields

I8"ij I j
82: 47f R —————d:cdyl =xFNMuniform+m

 

 

x [:(e-J'k\/(a-:v)2+(y’-:tI)2+(z-Z’)2 __ e-jk\/(a+w)2+(y’-y)2+(z-Z’)2)dyr.
(7.12)

Linear apodization in the y direction

The potential integral for the apodization function f(:z:',y’,z’) 2 y, applied to a

planar source is given by

—ij -—ij —ij
I6 I I I 6 I I e I I
d d' 2 . — d d d
[[83, 4,, x y f/St mm, mdy+//Sy47rR x y

Be —ij I
24779 j/fs—br —-——d:L‘ dy +yFNMum'form' (7.13)

 

 

 

For a triangular source, substitutinyg Eq. (7.5) into Eq. (7.13) yields

 

. 0
’8'.)de ’dz’ FNM j
3y 47TR I y :y uniform'i'm
"Cl/Al

134

 

 

 

_,-.. I_ 2 I_ 2 ,_I2 _,- I_.2 .2 .,_.,/2
e liq/<2: 22) +<L3(a:) y) +<~ z) _e Jk\/($ 2) +y +<~ ~> d3,

 

+_J;_ -02/A2 e—jk\/<x'—x>2+<L4(x'>—y)2+<z—z’>2
47rk 0

 

_- I_ 2 2 ,_I2
-e “((93 5”) +3” +(“ z) (133’. (7.14)

For a rectangular source, substituting Eq. (7.6) into Eq. (7.13) yields

 

’eTijd ’d ’— FNM j
3y 47rR :r y —y uniform'i'm

 

 

[0 e-jk\/(x’—:r)2+(b—y)2+(z—z’)2__e—jkfis’—x)2+(b+y)2+(z—z’)2 (135'.

(7.15)

Linear apodization in the z direction

The potential integral for the apodization function f (3:, ,y’ ,z’) 2 z’ applied to a

planar source is given by

 

I€_ij I I I ,
Sz 47rR d1: dy 2:. FNjwuniform' (7.16)

7 .1.3 Linear apodization for a subdomain

The apodization function is defined as f (23’ ,y’ ,z’ ) for a subdomain of the volume
source in the local coordinate system shown in Figure 7.2. The origin of the local
coordinate system is the projection of the observation point onto the base plane of

the subdomain. The base plane of the subdomain is either a triangular source which

135

 

is formed by three lines 3) 2 0, A111: + Bly + Cl 2 0 and A23: + 82y + Cg 2 0 or a

rectangular source with the width and height 2a, and 2b, respectively.

Linear apodization in the x direction

Consider the apodization function f (:r’ ,y’ ,z’ ) 2 l" for a subdomain of the volume
source in the local coordinate system shown in Figure 7.2. The potential generated

by the subdomain is the linear superposition of the following two integrals,

’63—ij I e—ij e—ij
ff/Vx 47TH dV2/ffl/(3: —a:) 47rR dV+///Vz 417R (1V. (7.17)

The first volume integral on the right side of Eq. (7.17) is computed as

’ e_ijdV h” ’ {Maw 718
///V($—x) 477R —/(; //3(2:—x) 47rR sz. (.)

For a triangular base plane, substituting Eq.(7.11) into Eq. (7.18), then exchanging

 

 

 

 

 

the order of integration and analytically evaluating the innermost integral gives

///V<w'—x>ii’;Rdv

’1 ' _ 7’ —C B -,z—z’ ~ z—z’
0 0

 

 

m: z

j _Cl/BI e“jk(R1+R2)

=_ jk(R1+R2) 2- 27
4m» 0 k2RgRg (8 (R2 Rll”)

 

e-jk(R1+R2)
_ .2 2 2
A R1122

 

(eikngz — 79ij1sz + jR1R2k3(R28ij2 — Rleijl))dy’

(7.19)

 

 

where R2 2 \/(L2(y’) — :1:)2 + (y’ — y)2 + 2:2 and
R1 = \/(L1<y'> — a)? + (y' — y)? + 32

136

 

For a rectangular base plane, substituting Eq.(7.15) into Eq. (7.18), then exchanging

the order of integration and analytically evaluating the innermost integral gives

///V<x'-x>::’:fdv

h - _ I b .‘ _ ’ - z—z’
"f ”(7171/ b<e“9kzzsz2—e-9k z Rb1>dy’>dz’
0 71' Z _

 

 

(ejk(R1+R2)(Rg _ 3%)3)

 

_ j /b e-jk(R1+R2)
‘_ 2 2 2
47k —b k R1122

e-jk(Rl+R2)
2 2 2
k R1R2

 

(eijQRgz -— eijl R¥z +jR1R2kz(R2eij2 — Rleij1))dy’
(7.20)
where R2 2 \/(a —— :1:)2 + (y’ - y)2 + 22 and R1 2 \/(a +1?)2 + (y’ - y)2 + 22.

 

 

The second volume integral on the right side of Eq. (7.17) is computed as

e—ij
ff/Vx R dV2xFNM. (7.21)

7 .1.4 Linear apodization in the y direction

 

The apodization function is defined as f (x’, y’, 2’) 2 y, for a subdomain of the volume
source in the local coordinate system shown in Figure 7.2. The potential generated

by the subdomain is the linear superposition of the following two integrals,

,e—ij I e—ij e—ij
///Vy 47rR dV2///V(y —y) 42R dV+///Vy 47rR dV. (7.22)

For a triangular base plane, the first triple integral on the right hand side of Eq.

 

 

 

(7.22) is given by,

137

 

 

 

ff/v 9’ —>y >—ijdV= j f—CQWe—JMRfiRfl<eJk<R1+R2)(R%—2 R¥>z>
477R 47rk 0 k2R%R%

e-jk(Ri+R2)
2 2 2
k R1122

 

(eij? R32 — eij1R%z + jR1R2k2(R2€ij2 — Rleij1))d:r’

 

j 0 e‘3k(R1+R3) jk(R1+R3) 2_ 2 '
2 2 6 (R3 31);.)
431/211 #3133

-jk(Rl+R3) - ,- ,- -
e k2R2R2 (63191321232: — eJkR1R2z + jR1R3kz(R3eJkR3 — 1216316121»de L;
1 3 S
(7.23)

where R2 2 \/(:1:’ - x)2 + (L4(:r’) — y)2 + Z2 , R1 2 \/(—:r’ - 2:)2 + y2 + Z2 and
R3 2 \/(:r’ - :1:)2 + (L3(:r’) - y)2 + 22.

For a rectangular base plane, the first triple integral on the right hand side of Eq.

 

 

 

 

 

(7.22) is given by

 

 

— kR - -'k R +R
[#wa _ye J dzv J [a e H 1 2ng ij(R1+R2)(2 422))
47rR 47rk a k2R%R%

 

(eijQRgz —~ eij1R%z + jR1R2k2(R2eij2 — Rleij1))d:r’,
(7.24)
where R2 2 \/(iL‘,—SE)2 +(b——y)2 +2:2 and R12 \/(:r’-—:r)2 +(b+y)2+22.

 

 

The second triple integral on the right hand side of Eq. (7.22) is equivalent to,

e—ij
[ffvy R dV2yFN.M. (7.25)

 

138

7 .1.5 Linear apodization in the z direction

The apodization function is defined as f (113’, y’, z’) 2 z, for a subdomain of the volume

source in the local coordinate system shown in Figure 7.2. The potential is given by

_ij Z M
ff/vzleer dV:/z’//SC———4WR dsdz’. (7.26)
0

Substituting the F N M expression for the surface source into Eq. (7.26) and exchang-

 

ing the order of integration yields

M ”3.2%

 

 

 

l ' ZZ Z, 2 2 2 .
Z , j N m2]! 6 _jk_ VQ +z +hb2— e_Jk(Z—Z,)d d I
=/zm. /“ Q2+h2. QZ
0 2:11 , b2
bzl
, _J
j N ’bz2 z , —2k-Z—22\/Q2+22+h2-_e-jk(z_z') ,
247%.: f2)!“ Q2+h2 dZdQ
”117771 0

 

z - —'k‘/ 2+z2+h2.
hb N (”2 22(1—jk\/Q2+z2+hl2n.—e 3 Q bl)

2 2 2 2
bz'izll k (Q +Z +hb7)

 

 

 

_ "k: _ —jkz
_(1 J k2 e )) dQ. (7.27)

Then, the 1D integral for the linear apodization in the z direction is obtained.

7 .1.6 The prism geometry

The prism source evaluated here is shown in Figure 7.3, and the coordinates of each

vertex of the prism are given in Table 7.1.

139

 

Table 7.1. Vertex locations defined for the prism source..

 

Vertex A B C D E F
(x, y, z) [m] (0,1,0) (0, 0,0) (1,0,0) (0,1,1) (0,0,1) (1,0,1)

 

 

 

 

 

 

 

 

 

 

 

7 .1.7 Global and local systems

The global coordinate system is denoted by (:c,y,z), and the rotated coordinate
system is denoted by (X, Y, Z). The rotation is characterized by Euler angles using
the Z-X-Z convention, where

a is the angle between the x—axis and the intersection of the 3:3; and the X Y coordinate
planes.

B is the angle between the z-axis and the Z —axis.

 

7 is the angle between the intersection of the 1133/ and the X Y coordinate planes and
the X —axis.
The relationship between the global coordinates (:1:, y, z) and the rotated coordi-

nates (X, Y, Z) is given by

 

 

 

 

 

 

, -
X
Z
r- - r- 1
cosacos'y — cosfisinasirw sinacosv + cosflcosasin'y sinfisin'y x
= —cosasin7 —— cosfisinacosv cosficosacos'y — sinasiny sinficos'y y
sinflsina —sin,3cosa cosfi z

140

 

 

 

 

Figure 7.3. The prism geometry, where vertex B is coincident with the origin, and
vertices C, A, and E are located on the :r, y, and z axes, respectively.

The plane defined by the vertices ADE B

The Euler angles for the plane defined by the vertices ADEB are (a, fin) =
(7r / 2, 7r / 2, 0). After rotating the global coordinate system :ryz according to the Euler
angles ((1, fl, '7), the origin of the rotated system X 1Y1 Z 1 is translated to the center
of the rectangle ADEB which is denoted by ($1,y1,z1). The values of (331,311, 21)

are given by 6 %, 0). Substituting the Euler angles into Eq. (7.28) yields

141

 

 

 

 

 

 

 

 

 

‘ Y]
E

u D

l

I

: >-
F : Ol (x1.y1,zl) X]

i 21

I

,L ................... -
x B A y
C

Figure 7.4. The local coordinate system for the plane ADE B.

 

 

 

 

 

 

 

 

 

 

- 1 - - _ . _ - - _
X1 0 1 0 :1: x1 y —— %
Y1 = 0 0 1 y ‘ y1 = Z—% (7'29)
21 j 1 O 0 z 21 :1:
Rewriting Eq. (7.29) yields
' ' - T
£1? 21
y = XI + % . (7.30)
2: Y1 + 1
_ J 2 _

 

 

 

 

142

 

Y2

 

 

02 (x2.y2,22)

 

 

 

,’ ZZ

 

 

Figure 7.5. The local coordinate system for the plane ADF C .

The plane defined by the vertices ADF C

The Euler angles for the plane defined by the vertices ADF C is (013,7) =
(37r / 4, 7r / 2, 0). After rotating the global coordinate system xyz according to the Euler
angles (a, ,8, 7), the origin of the rotated system X2Y222 is translated to the center
of the rectangle ADF C which is denoted by ($2,312, 22). The values of (332,342, 22)
are given by (0, %, é). Substituting the Euler angles into Eq. (7.28) yields

143

 

 

0

fin

 

L

_—\/§/2 73/2 ol

Rewrite Eq. (7.31) yields

0 1
fi/zo
1'

y =
Z

 

 

 

 

_

 

 

 

:L‘ 1‘2
y - 312
Z Z2

Y2+%

The plane defined by the vertices BEF C

 

pgfiagfi+i
§X2+3§Z2+%

 

 

%% ~

 

(7.31)

(7.32)

The Euler angles for the plane defined by the vertices BEFC is (a, fi, ’7) = (7r, 1r / 2, 0).

After rotating the global coordinate system :cyz according to the Euler angles (0, fl, '7),

the origin of the rotated system X 3Y3Z3 is translated to the center of the rectangle

BEF C which is denoted by ($3,y3,z3). The values of (x3,y3,z3) are given by

(—%, 1%, 0). Substituting the Euler angles into Eq. (7.28) yields

 

_

 

d

 

Rewriting Eq. (7.33) yields

N

 

O

1

 

1 P

 

 

 

1' 2133

y 93

 

 

Z Z3

144

 

 

p-

 

—:r+%

2—3

T

 

(7.33)

(7.34)

 

 

 

 

 

E
n D
l
Y3 .
F X3 :
: Z3
()3(x3,y3,z3) )_ ___________________ >
B A y

 

 

Figure 7.6. The local coordinate system for the plane BEF C.

7.1.8 Evaluating potential integrals generated by each sub-

domain
To calculate the potential integrals generated by the volume source shown in Figure
7.3, the volume is subdivided into five subdomains. The potential generated by the

prism , denoted by (Pa, is obtained by superposing the contributions generated by

each subdomain according to

145

 

(Amx + Bmy + sz + Dm)a
lAmSL‘ + 87111] + CmZ + Dml

(7 35)

 

 

ME

ij
Rdx’dy’dz— —

all/f

where M = 5 is the number of faces that the volume has, Amx+Bmy+sz+Dm = 0

m:l

denotes the general form of the plane that passes through each face of the prism, and

(Dam is the potential on each subdomain. The expressions of (Dam is

cram: ff] f(x’, Y’ Jig—Z) ————dX ’,dY’dz’ (7.36)

for each subvolume, where f (X I , Y’, Z I ) is the apodization function in the local coor-
dinate system, which can be obtained by substituting the expressions for the :1:], y’, z,

terms defined by Eqs. (7.30), (7.32), and (7.34).

7.1.9 Error calculations

For potential calculations, the numerical error 77(23, y, z; k) is defined as the nor-
malized difference between the reference potential (Dre f(:c, y, z) and the computed

potential <I>(:1:, y, :5) according to

l¢<x7yv Z) _ (1)7-ef(1',y, Z)|

 

77(27, y, z; k) = (7-37)
max|¢ref(x,y,z; k)!
The maximum error Umax is defined as
7,777,017 : 11,1:ng 77(27, y: Z): (7'38)
and the number of significant digits n is given by
n = [0910(77ma1r) (7-39)

146

Table 7.2. Simulation parameters that achieve between 2 and 5 significant digits using
the FNM and the singularity cancellation method to evaluate a potential integral with
linear apodization evaluated on a 3D grid. Parameters listed include the number of
significant digits achieved and the ratio of the computation time relative to the time
required for the F NM. (a) Significant digits 2 and 3. (b) Significant digits 4 and 5.

(a)

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Significant digits 2 3
FNM Singularity FNM Singularity
cancellation cancellation
Time (seconds) 0.095 0.83 0.11 1.70
Computaton time relative to FNM 1x 8.74>< 1x 15.45x
(1))
Significant digits 4 5
FNM Singularity FNM Singularity
cancellation cancellation
Time (seconds) 0.16 4.96 0.25 8.69
Computation time relative to FNM 1x 31.00x 1x 34.76x

 

 

 

 

 

 

 

7 .2 Results

The potential integral with linear apodization is evaluated on a 3D grid. The x
and y directions extend from -0.5m to 1.5m with an interval of 0.02m, the z direction
extends from 0 to 1.5m with an interval of 0.03m. The reference potential is computed
using the singularity cancellation method evaluated with 1000 Gauss abscissas in each
direction. Figure 7.7 contains a plot of the number of significant digits achieved with
the F N M and the singularity cancellation method evaluated on a 3D grid as a function
of computation times. Table 7 .2 summarizes the ratio of the computation times
obtained with the F N M and the singularity cancellation method when each achieves
between 2 and 5 significant digits on a 3D computational grid. For this computational
grid, the FN M is 8.74, 15.45, 31.00, 34.76 times faster than the singularity cancellation

method when each achieves 2, 3, 4, and 5 significant digits, respectively.

147

 

 

 

 

 

 

  

    

 

 

 

3 I Singularity Cancellation

'5, 0 FNM g

E 2 g :I

812_. ................. .........

E E E:

c 1 :‘

_g) 9_ ............................... g .........

‘0 : 5 3:

“5 . 3 3 33

'- 6_. ..........

.8 . 5 3%

E 3 ......... ._.-“'

3 .0 .. ...

= ”ééééi i "1"" 3]..- ...
o I 1:22:13; 'I-Zliilii I IIIIiIIi I 22123::
0.01 0.1 1 10 100

computation time (second)

Figure 7.7. The number of significant digits achieved in calculations 0d the potential
integrals linear apodization over the prism plotted as a function of computation time
for evaluations on a 3D grid. The results show that the FNM achieves a larger
number of significant digits than the singularity cancellation method for the same
computation time.

7 .3 Discussion

7 .3.1 Other volume sources

The potential integral with linear apodization in Eq. (7.1) is a triple integral and
the 1D FNM integral reduces this triple integral to the evaluation of a small number
of single integrals. The 1D F NM integral converges much faster than the singularity
cancellation method, and the 1D FNM integral needs much less computation time to
achieve a given accuracy than the singularity cancellation method. The fast nearfield

method calculations developed here are easily extended to any volume source with

148

polygonal faces. First, the volume source is subdivided into subdomains, and the
number of subdomains is equal to the number of polygonal faces. Then, the potential
generated in each subdomain is computed in the local coordinate system, where the
Euler angle is first obtained for each of the base planes. The total potential generated
by the volume source is obtained by superposing the potentials produced by each of

the subdomains.

7 .3.2 Higher order polynomials

The 1D F NM expressions for volume potentials with linear apodization are derived
here. The FNM expressions of the volume potential integrals with higher order poly-
nomial apodization can be derived using the same strategy. The Euler angle is first
obtained for each of the base planes, and the apodization function in the global co—
ordinate system is translated into the the function in the local coordinate system.
Then, after obtaining the fast nearfield method for each apodization in the local co-
ordinate system and superposing the contributions from each apodization function in

local coorinate system yields the total potential.

7 .4 Conclusion

1D FNM expressions for calculating volume potentials with linear apodization are
derived and evaluated. The 1D FNM integral for a potential integral is first given,
and FNM expressions for a planar source with linear apodization is then derived. The
volume source is then subdivided into several subdomains including subpyramids
and subtetrahedrons. The potential generated by each subdomain is computed in
the local coordinate system and the total potential over the volume source is the
superposition of the potential generated by each subdomain. After the 1D FNM

expressions are obtained, the results obtained with the FNM are compared with the

149

 

 

 

singularity cancellation method on a large 3D observation grid. The results indicate
that, for the 3D grid used here, the FNM reduces the computation time by a factor
of 8.74 to 34.76 relative to the singularity cancellation method when each achieves
between 2 and 5 significant digits on a 3D grid. Thus, the F NM is an ideal method

for linear apodized potentials.

150

 

CHAPTER 8

Conclusion

This thesis derives and evaluates fast expressions for calculating pressures gener-
ated by planar pistons and potentials generated by volume sources. The 1D fast
nearfield method (FNM) expressions are compared to several existing methods in-
cluding the implulse response method, the Rayleigh-Sommerfield integral, the Field
II program and the singularity cancellation method. Chapter 2 introduces 1D fast
nearfield method expressions for the time-harmonic and transient pressures generated
by triangular pistons. These fast nearfield method expressions remove singularities
from a 1D integral and therefore converge very quickly. The transient calculations
are further accelerated by time-space decomposition. The fast nearfield expressions
for a triangular source are compared with the impulse response method, the Field 11
program, and the smoothed impulse response. The comparison results indicate that
the fast nearfield method achieves smaller errors than the other three methods in less
time.

Analytical 2D integral expressions for fast calculations of time—harmonic and tran-
sient nearfield pressures generated by apodized rectangular pistons are given in Chap-
ter 3. A simplified time space decomposition method is also introduced to further

reduce the computation time for transient pressure fields. The results are compared

151

 

with the Rayleigh-Sommerfeld integral and the Field II program to show that the 2D
FNM integral converges much faster than the Rayleigh-Sommerfeld integral and the
Field 11 program.

As a special case of Chapter 3, A fast nearfield method for calculating pressure
generated by a polynomial apodized rectangular pistion is obtained based on the in-
stantaneous impulse response in Chapter 4. Two kinds of apodization functions are
considered in the derivation process. A trigonometric transform of the integrand is
performed and the order of integration is exchanged to obtain the 1D fast nearfield
method expressions. The results show that the 1D polynomial apodized FNM con-
verge faster than the Rayleigh-Sommerfeld integral and the 2D apodized FNM.

Chaper 5 introduces a 1D fast nearfield method for the calculations of uniformly
excited volume potential integrals. The results are compared with the singularity
cancellation method at six different observation points and a volume grid, and the
results show that the fast nearfield method needs less sample points than the singu-
larity cancellation method to achieve a given number of significant digits. A 1D fast
nearfield method to calculate potentals generated by surface integrals and volume
integrals with polynomial apodization are introduced in Chapters 6 and 7. These
expressions remove the singularities from the Rayleigh—Sommerfeld method integral.
The results compared with the singularity cancellation method indicate that the 1D
fast nearfield method achieves much faster convergence with much less computation

time.

152

BIBLIOGRAPHY

[1] F. Oberhettinger, ”On transient solutions of the ’baffled piston’ problem,” Jour-
nal of Research of the National Bureau of Standards—B. Mathematics and Math-
ematical Physics, 65B(1):1-6, 1961.

[2] P. R. Stepanishen, ” Transient radiation from pistons in an infinite planar baffle,”
J. Acoust. Soc. Am., 49(5):1629—1638, 1971.

[3] J. C. Lockwood and J. G.Willette, ”High-speed method for computing the exact
solution for the pressure variations in the nearfield of a baffled piston,” J. Acoust.
Soc. Am., 53(3):735-741, 1973.

[4] J. A. Jensen, “Ultrasound fields from triangular apertures,” J. Acoust. Soc. Am.,
100(4), pp. 2049-2056, 1996.

[5] J. Zemanek, “Beam behavior within the nearfield of a vibrating piston,” J.
Acoust. Soc. Am., 49(1):181-191, 1971.

[6] K. B. Ocheltree and L. A. Frizzell, “Sound field calculation for rectangular
sources,” IEEE Trans. Ultrason. Ferroelect. Freq. Contr., 36(2):242—248, 1989.

[7] R. J. McGough, T. V. Samulski, and J. F. Kelly, ”An efficient grid sectoring
method for calculations of the nearfield pressure generated by a circular piston,”
J. Acoust. Soc. Am. Vol 115 (5), Pt. 1, May 2004, pp. 1942-1954.

[8] R. J. McGough, “Rapid calculations of time-harmonic nearfield pressures pro—
duced by rectangular pistons,” J. Acoust. Soc. Am., 115(5), Pt. 1, 2004.

[9] J. A. Jensen, “A new calculation procedure for spatial impulse responses in
ultrasound,” J. Acoust. Soc. Am. 105(6), 3266-3274 (1999).

[10] J. A. Jensen, “Field: A program for simulating ultrasound systems, ” Med.
Biol. Eng. Comp., 10th Nordic-Baltic Conference on Biomedical Imaging, Vol.
4, Supplement 1, Part 12351-353, 1996b.

[11] J. A. Jensen and N. B. Svendsen, “Calculation of pressure fields from arbitrarily
shaped, apodized, and excited ultrasound transducers,” IEEE Trans. Ultrason.,
Ferroelec, Freq. Contr., 39:262-267, 1992.

153

[12] D. Chen, J. F. Kelly, and R. J. McGough, “A fast nearfield method for calcula-
tions of time-harmonic and transient pressures produced by triangular pistons,”
J. Acoust. Soc. Am. 120 (5), 2450—2459 (2006).

[13] X. Zeng and R. J. McGough, “Evaluation of the angular spectrum approach for
simulations of nearfield pressures,” J. Acoust. Soc. Am. 123(1), 68-76 (2008).

[14] T. D. Mast, “Fresnel approximations for acoustic fields of rectangularly symmet-
ric sources,” J. Acoust. Soc. Am. 121 (6), 3311-3322 (2007).

[15] G. R. Harris, “Review of transient field theory for a baffled planar piston,” J.
Acoust. Soc. Am. 70(1), 10—20 (1981).

[16] M. Greenspan, “Piston radiator: Some extensions of the theory,” J. Acoust. Soc.
Am. 65(3), 608-621 (1979).

[17] C. Lafon, F. Prat, J. Y. Chapelon, F. Gorry, J. Margonari, Y. Theillere, and D.
Cathignol, “Cylindrical thermal coagulation necrosis using an interstitial appli-
cator with a plane ultrasonic transducer: in vitro and in vivo experiments versus
computer simulations,” Int. J. Hyperthermia 16, 508-522 (2000).

[18] D. R. Daum, and K. Hynynen, “A 256-element ultrasonic phased array system
for the treatment of large volumes of deep seated tissues,” IEEE Trans. Ultrason.,
Ferroelect. Freq. Contr. 46(5), 1254-1268 (1999).

[19] Y. Li and J. A. Zagzebski, “A frequency domain model for generating B-mode
images with array transducers,” IEEE Trans. Ultrason., Ferroelect. Freq. Contr.,
46(3), 690-699 (1999).

[20] G. Cincotti, G. Cardone, P. Gori, and M. Pappalardo, “Efficient transmit beam-
forming in Pulse-Echo ultrasonic imaging,” IEEE Trans. Ultrason., Ferroelect.
Freq. Contr. 46(6), 1450-1458 (1999).

[21] Y. Lin, J. M. Dodson, J. D. Hamilton, J-U. A. Kluiwstra, C. Cain, and K. Grosh,
“Theory and experiment for the design of piezoeletric element for phased arrays,”
Proceedings of Ultrasonics Symposium, 2, 1697-1700 (1997).

[22] J. L. Prego Borges, F. Montero de Espinosa, J. Salazar, J. Garcia-Alvarez, J.
A. Chavez, A. Turo, and M. J. Garcia-Hernandez, “Diffraction aperture non-
ideal behaviour of air coupled transducers array elements designed for NDT,”
Ultrasonics, 44, Supplement 1, e667—e672 (2006).

[23] B. M. Notaros and B. D. Popovic, “Optimized entire-domain momentmethod
analysis of 3D dielectric scatterers,” Int. J. Numer. Model, vol. 10, pp. 177—192,
1997.

154

 

[24] S. M. Rao, D. R. Wilton, and A. W. Glisson, “Electromagnetic scattering by
surfaces of arbitrary shape,” IEEE Trans. Antennas Propagat., vol. AP-30, pp.
409-418, May 1982.

[25] C. M. Butler, “Evaluation of potential integral at singularity of exact kernel in
thin-wire calculations,” IEEE Trans. Antennas Propag., vol. AP-23, pp. 293-295,
Mar. 1975.

[26] D. H. Werner, J. A. Huffman, and P. L. Werner, “Techniques for evaluating the
uniform current vector potential at the isolated singularity of the cylindrical wire
kernel,” IEEE Trans. Antennas Propag., vol. 42, pp. 1549—1553, Nov. 1994.

[27] L. Rossi and P. J. Cullen, “On the fully numerical evaluation of the linear-
shape function times the 3-D Green’s function on a plane triangle,” IEEE Trans.
Microwave Theory Tech., vol. 47, pp. 398-402, Apr. 1999.

[28] T. F. Eibert and V. Hansen, “On the calculation of potential integrals for linear
source distributions on triangular domains,” IEEE Trans. Antennas Propag., vol.
43, pp. 1499-1502, Dec. 1995.

[29] M. A. Khayat and D. R. Wilton, “Numerical Evaluation of singular and near-
singular potential integrals”, IEEE Trans. Antennas Propag., vol. 53, no. 10,
2005.

[30] M. G. Duffy, “Quadrature over a pyramid or cube of integrands with a singularity
at a vertex,” SIAM J. Numer. Anal., vol. 19, no. 6, pp. 1260—1262, 1982.

[31] M. A. Khayat and D. R. Wilton, “Numerical Evaluation of singular and near-
singular potential integrals”, IEEE Transactions on antennas and propagations,
vol. 53, no. 10, 2005.

[32] A. H. Stroud, Approximate Calculation of Multiple Integrals. Englewood Cliffs,
NJ: Prentice-Hall, 1971, pp. 31—32.

[33] R. D. Graglia, “Static and dynamic potential integrals for linearly varying source
distributions in two- and three-dimensional problems,” IEEE Trans. Antennas
Propag., vol. AP-35, pp. 662—669, Jun. 1987.

[34] M. Frigo and S. G. Johnson, “FFTW: An adaptive software architecture for the
FFT,” Proceedings of the ICASSP, vol 3, no. 1, pp.1381-1384, 1998.

[35] J. D’hooge, J. Nuyts, B. Bijnens, B. De Man, P. Suetens, J. Thoen, M.-C.
Herregods, and F. Van de Werf, “The calculation of the transient near and far

field of a baffled piston using lwo sampling frequencies,” J. Acoust. Soc. Am.
102(1), 78-86 (1997).

155

 

[36] J. F. Kelly and R. J. McGough, “A time-space decomposition method for calcu-
lating the near field pressure generated by a pulsed circular piston,” IEEE Trans.
Ultrason. Ferroelect. Freq. Contr., 53(6), pp. 1150-1159, 2006.

[37] J. F. Kelly and R. J. McGough, “An annular superposition integral for axisym-
metric radiators,” J. Acoust. Soc. Am. 121(2), 759-765 (2007).

[38] T. L. Szabo, “Generalized fourier transform diffraction theory for parabolically
anisotropic media”, J. Acoust. Soc. Am., vol. 63(1), pp. 28-34, 1978.

[39] X. Zeng, Private communication, 2008.

[40] J. J. Wen and M. A. Breazeale, “A diffraction beam field expressed as the su-
perposition of Gaussian beams,” J. Acoust. Soc. Am. 83(5), 1752-1756 (1988).

[41] D. Ding, Y. Zhang, and J. Liu, “Some extensions of the Gauss beam expansion:
Radiation fields of the rectangular and the elliptical transducer,” J. Acoust. Soc.
Am. 113(6), 3043-3048 (2003).

[42] B. D. Cook and W. J. Arnoult III, “Gaussian-Laguerre/ Hermite formula for the
nearfield of an ultrasonic transducer, ” J. Acoust. Soc. Am. 59(1), 9—11 (1976).

[43] D. Chen and R. J. McGough, “A 2D fast nearfield method for calculating
nearfield pressures generated by apodized rectangular pistons”, accepted by J.
Acoust. Soc. Am., 124(3), pp. 1526—1537,2008.

[44] G. Scarano, N. Denisenko, M. Matteucci, and M. Pappalardo, “A new approach
to the derivation of the impulse response of a rectangular piston”, J. Acoust.
Soc. Am. 78 (3), 1985, pp. 1109—1113.

[45] Xiao—Chun Nie, Le-Wei Li, Ning Yuan, Tat Soon Yeo, and Yeow-Beng Gan,
“Procorrected—FFT solution of the volume integral equation for 3—D inhomege-
neous dielectric objects,” IEEE Trans. Antennas Propagat., vol. 53, pp. 313—320,
January 2005.

[46] D. Chen and R. J. McGough, “A Fast Nearfield Method for Calculating Pres-
sures Generated by Rectangular Pistons with Polynomial Apodization”, IEEE
Ultrasonics Symposium, pp. 134-137, 2007.

156