An ellipsoidal acoustic infinite element

Computer methods in applied mechanics and engineerlng

ELSEWIER

Comput.

Methods Appl. Mech. Engrg. 164 (1998) 49-76

An ellipsoidal acoustic infinite element David S. Burnett*,

Richard L. Holford

Bell Laboratories, Lucent Technologies, Whippany, NJ 07Y81, USA

Abstract In previous papers the authors presented new 3-D time-harmonic prolate and oblate spheroidal acoustic infinite elements, based on new prolate and oblate spheroidal multipole expansions, for modeling acoustic fields in unbounded domains. Here, we present the development of an ellipsoidal infinite element, which is the logical generalization of those elements. This development is also based on a new multipole expansion as well as a new system of ellipsoidal coordinates. The element stiffness, radiation-damping and mass matrices are developed and presented in sufficient detail to enable their software implementation. Both the new coordinate system and the new element include the previous spheroidal coordinate systems and spheroidal elements as limiting cases. Therefore, all the previously reported performance data for the spheroidal elements apply to this ellipsoidal element when used in spheroidal form. Since the three axes of an ellipsoid can be chosen independently, an ellipsoid can circumscribe any structural shape at least as closely, and generally more closely, than a prolate or oblate spheroid. The resulting reduction in size of the finite computational domain will result in even greater computational speeds than those already reported for the spheroidal elements. The element may be used to model problems in free-space (471 steradians), half-space (2rr), quarter-space (7~) or eighth-space @r/2). Since this infinite element provides maximum computational efficiency for structures of all shapes, 0 1998 Elsevier Science S.A. All software element libraries would need only this one element for all problems in unbounded domains. rights reserved.

1. Introduction In the field of computational structural acoustics, the problem of efficiently modeling the acoustic field in unbounded domains has motivated an extraordinary amount of theoretical research and algorithm development since the 1960s. Many techniques have evolved (see overviews in [1,2]). The boundary element method (BEM) was the method of choice for most of this period, but currently there is a resurgence of interest in the infinite element method (IEM). To clarify the subject matter, an infinite element is a finite element [3] that covers a semi-infinite sector of space. In general finite element semantics, a finite element is a region of space in which an unknown field variable is represented as a trial function that is defined over the entire element. (The trial function comprises an expansion in terms of known basis functions, each multiplied by a coefficient, often called a degree of freedom, whose value is to be determined in the ensuing numerical computation.) Thus, an acoustic infinite element is a finite element in which the field variable, typically acoustic pressure, is represented by a trial function over an entire semi-infinite sector of space. We now take a brief look at the history of acoustic infinite elements. A comprehensive chronology of infinite elements, acoustic and otherwise, from 1973 to 1992, is provided by Bettess [4]; here we mention only the milestones in that period and developments since then. Following the first attempt in 1973 [5], Bettess [6] introduced a trial function with exponential decay in the radial direction. This did not conform to the correct decay law but nevertheless was able to provide some degree of approximation. During the 1980s there were many papers on the so-called ‘mapped’ infinite elements [7] and ‘wave-envelope’ infinite elements [S]. The mapped elements used the correct radial decay (i.e. eP”‘/ r in 3D) but their conditions for convergence were not * Corresponding

author.

0045.7825/98/$19.00 0 1998 Elsevier Science S.A. All rights reserved PII: SOO4S-7825(98)00046-2

50

D.S. Burnett, R.L. Holford I Comput. Methods Appl. Mech. Engrg. 164 (1998) 49-76

understood. The wave-envelope elements used complex conjugate test functions (a Petrov-Galerkin approach), along with mapping, to simplify the numerical integration by eliminating the oscillatory behavior in the radial direction. This work is still continuing [9]. In 1994 Burnett [l] used a radically different approach based on a new multipole expansion for the exact field in terms of prolate spheroidal coordinates, and Holford [lo] provided a proof of convergence for this expansion, which defined the proper conditions for use of the infinite element. This approach used a Bubnov-Gale&in formulation, i.e. non-conjugated test functions. These prolate spheroidal infinite elements yielded an extraordinary performance improvement over the BEM: over 400 times faster to achieve the same solution to the same accuracy. The elements retain their accuracy and numerical stability to dimensionless frequencies as high as ku = 100 [I 11. This approach has been extended to oblate spheroidal infinite elements [ll]. The high accuracy of this formulation has enabled the infinite elements to be placed very close to a structure, viz. typically within half a wavelength, to achieve pointwise errors of about l%, or smaller if desired. There have been several studies of the convergence and performance of the Burnett infinite element, both in non-conjugated and conjugated forms [ 12- 171, application of wave-envelope elements with spheroidal meshes [ 181, and modification of the Burnett element for improved performance with iterative solvers [ 191. All the above infinite elements are time-harmonic, i.e. frequency-domain. Interestingly, infinite elements have been recently extended to transient analyses, i.e. time-domain, by inverse Fourier transforming the Burnett element [20] and the wave-envelope element [21]. At the present time we see two prominent trends in infinite element formulations: (1) the use of multipole expansions in curvilinear coordinate systems; and (2) the use of both non-conjugated and conjugated formulations of multipole-based elements, with continuing investigations of their relative merits. The motivation for the current work was to create a multipole-based infinite element that would be more general than the authors’ previous prolate and oblate spheroidal infinite elements, yet would include them as limiting cases. Such as element would, ipso facto, perform identically to the previous elements when used in those same configurations, but in all other configurations would generally be faster, sometimes dramatically faster. In addition, just one infinite element would be sufficiently flexible to be applicable to virtually all structural shapes, thereby simplifying software implementation. It will be recalled from [l ,lO or 1 l] that multipole-based infinite elements must lie outside a minimal closed surface that just circumscribes the structure. Fig. 1 shows the closed coordinate surfaces belonging to four of the eleven so-called confocal quadric systems [22]. (The other seven systems have one of these closed surfaces or no closed surfaces at all.) The ellipsoidal surface is the most general of all, including the other three as limiting cases. Its three orthogonal principal cross sections may all be ellipses of different eccentricities, including circles as degenerate forms (which are the spheroids and sphere). Since the lengths of the three axes of an ellipsoid can be chosen independently, an ellipsoid can circumscribe any structural shape at least as closely, and generally more closely, than a prolate or oblate spheroid, resulting in a reduction in the size of the finite computational domain and hence a faster computation time.

(d)

Fig. 1. The closed coordinate ellipsoids.

surfaces

of the four coordinate

systems

for: (a) spheres;

(b) prolate

spheroids;

(c) ablate

spheroids;

(d)

D.S. Burnett, R.L. Holjord

I Comput. Methods Appl. Mech. Engrg.

164 (1998) 49-76

51

Below, Section 2 derives a system of ellipsoidal coordinates (r, 0, 4) relative to an ellipsoid, S, surrounding the structure of interest, that will serve as the base for constructing the infinite elements. In this system the surfaces r = constant are confocal ellipsoids (with r serving as the ‘radial’ variable), and the surfaces 0 = constant and 4 = constant are confocal hyperboloids (with 6, 4 serving as ‘angular’ variables). These coordinates have the property that they reduce to the spheroidal coordinates employed in [ l,lO, 1 I] when the ellipsoid has two of its axes equal. In Section 3, we show how the radiation pressure field, p, exterior to S can be represented by a multipole expansion in ellipsoidal coordinates in the form

(1.1) where S is the ellipsoid r = r,. Section 4 shows how a truncated form of this expansion is incorporated in an infinite element. The underlying global and differential geometry and governing physics equations are presented first, follow by a detailed development of all mathematical expressions, culminating in the final expressions for the element matrices. Section 5 enumerates the many advantages that this ellipsoidal element affords. No numerical results are presented in this paper. However, all performance data for the spheroidal elements in [ l,lO,l l] apply to this ellipsoidal element since the latter includes the former as limiting cases. This paper is a logical sequel to [ 111. There is a strong parallelism between the theoretical development of the ellipsoidal and spheroidal infinite elements (although algebraic details are virtually all different). This parallelism is emphasized by deliberately maintaining the same format, wherever possible, so that the two can be treated as companion papers.

2. Ellipsoidal

coordinates

The basis for our development of an ellipsoidal infinite element is a multipole expansion for in a region exterior to an ellipsoid surrounding the structure of interest [23]. Our purpose in derive a suitable system of ellipsoidal coordinates for effecting this. We begin by orienting Cartesian coordinates (x, y, z) with the x-axis aligned with the longest structure and the z-axis aligned with the shortest. Then, choosing the origin of coordinates at circumscribing ellipsoid, the latter will be given by the equation

the pressure field this section is to dimension of the the center of the

(2.1) where 2a, 26 and 2c are the major, intermediate,

and minor axes of the ellipsoid,

a>b>c.

respectively,

such that (2.2)

(For the time being, we demand strict inequality between a, b and c, as shown in Eq. (2.2). Subsequently, these inequalities will be relaxed.) Fig. 2 shows the positive octant of the ellipsoid defined by Eqs. (2.1) and (2.2), together with the foci defined by its principal cross-sections, as explained below. In the plane y = 0, the cross-section of the ellipsoid given by Eq. (2.1) is an ellipse with foci at x = +,L where f.2=&&

(2.3)

In the plane z = 0, the cross-section

is an ellipse with foci at x = kg,

where

g2 = a2 - b2. And in the plane x = 0, the cross-section h2 = b2 -

c=.

(2.4) is an ellipse with foci at y = -+h, where (2.5)

D.S. Burnett, R.L. Holford I Comput. Methods Appl. Mech. Engrg. 164 (1998) 49-76

52

Y b

Fig. 2. Positive octant of ellipsoid,

showing

semi-axis

lengths, positive foci, and their interrelationships.

From Eq. (2.2), we note that

(2.6)

f’g, and from Eqs. (2.3), (2.4) and (2.5), we have

h2=f2-g2,

(2.7)

so that the positions of only two sets of foci are independent. Now, for any t, we see that the surface X2

-=z2

Y2

1

(2.8)

a2+t+b2+t+C2+t is a quadric conformal the form

with the given ellipsoid.

For given values of x, y and z, this equation may be rewritten in

E(t) = 0 with E(t) = (a2 + t)(b2 + t)(c2 + t) - x2(b2 + t)(c2 + t) - y2(c2 + t)(a2 + t) - ~‘(a’ + t)(b2 + t) qt

- h)(t - p)(t - v) ,

(2.9)

where we have denoted the roots of the equation E(t) = 0, considered as a cubic in t, by t = A, ,u, v (in descending order). A consideration of the sign of E(t) at t = -a’, -b2, -c2 and ~0, shows that these three roots are distinct, and are such that -_a’
putting

x2 =

t = -a’,

-b2,

(2.10) and -c2

(a2+ A)(a2+ p)(a2+

v)

(b2 - a2)(c2 - a’) y* =

’

(b2 + A)(b2 + p)(b2 + v)

(2.11)

’

(c2 - b2)(a2 - b2) z2 =

in turn in Eq. (2.9), yields

(c* + A)(c2 + p)(c2 + V) (a2 - c2)(b2 - c2)

’

In Eqs. (2.1 l), (A, /_L,V) are ellipsoidal

coordinates

for the point (x, y, z). Inspection

of Eq. (2.8) shows that the

D.S. Burnett, R.L. Holford I Comput. Methods Appl. Mech. Engrg. 164 (1998) 49-76

53

quadrics given by t = A (A > -c2) are ellipsoids, with A < 0 corresponding to the interior of Eq. (2.1) and A > 0 corresponding to the exterior; the quadratics given by t = ,u (-b* < p < -c2) are hyperboloids of one sheet; and the quadrics given by t = v (-a2 < v < -b2) are hyperboloids of two sheets. These three sets of quadrics are mutually orthogonal. Although the coordinates (A, p, V) are useful in certain applications, Eqs. (2.11) show that their relationship to (x, y, z) is not one-to-one. To make it so, and simultaneously to define a radial ellipsoidal coordinate for which the multipole expansion to be derived in Section 3 has the desired convergence properties, we introduce new ellipsoidal coordinates (r, 8, 4) according to A = r2 - L12) (2.12)

p = -b2 cos20 - c2 sin28, v = -a2 sin24 - b2 cos’q5.

Substituting Eqs. (2.12) in Eqs. (2.11) yields X = r( 1 - 4 cos28)“2 cos 4 ) (2.13)

y = (r2 - g2)“2 sin 8 sin 4 , z = (r2 -f2)“’

cos e( 1 -p cos2+)1’2,

r = con&.

(b)

Fig. 3. (a) Ellipsoidal

coordinates;

(b) exploded

view of surfaces

r = constant,

0 = constant

and 4 = constant

D.S. Burnett, R.L. Holford

54

I Comput. Methods Appl. Mech. Engrg.

164 (1998)

49-76

where f and g are given by Eqs. (2.3) and (2.4), and p=&f2,q=

(2.14)

l-p.

Here, the ranges for the variables

r, 0, and 4 are given by r ~5

and either

OSeS2n,

OGf$S?r,

(2.15a)

OGl9Cn,

OCC#IG2lr.

(2.15b)

or

The coordinate surfaces given by r = const. (ellipsoids), 0 = const. (one-sheeted hyperboloids), and $I = const. (two-sheeted hyperboloids), are illustrated in Fig. 3. The ellipsoid corresponding to Eq. (2.1) is given by r = a. We note that the radial variable r is the same as the variable p used by Hobson [24], but the introduction of the angular variables 6 and C#Iappears to be new. 2.1. Spheroidal

coordinates

The strict inequalities given in Eq. (2.2) were necessary to render the coordinate transformation (2.11) in terms of A, p, V, well defined. With the introduction of the coordinate transformation

given by Eqs. given by Eqs.

b-4

Fig. 4. Sequence of ellipsoidal coordinate (e) 5:4.9:1; (f) 5:5:1 (oblate spheroid).

surfaces for values of a:b:c given by (a) 5: 1:l (prolate spheroid);

(b) 5:l .l:l;

(c) 5:2:1; (d) 5:4: 1;

D.S. Burnett, R.L. Holford

I Comput. Methods Appl. Mech. Engrg.

(2.13),

in terms of r, 0, 4, these inequalities can now be relaxed. In particular, (i) b = c, corresponding to a prolate spheroid; and (ii) a = b, corresponding

2.1.1.

164 (1998) 49-76

55

we can consider the special cases: to an oblate spheroid.

Prolate spheroid

When b = c, we have f = g, implying x=rcos+,

y = (r2 -f

p = 1 and q = 0. Eqs. (2.13) reduce to

2)‘12 sin 8 sin 4 ,

z = (r2 -f 2)“2 cos 8 sin 4 ,

for which the choice given by Eq. (2.15a) is appropriate, sections (y/z = const.). 2.1.2. Oblate spheroid When a = b, we have g = 0, implying x = r sin 0 cos f#~,

Note that the surfaces 0 = const. degenerate

(2.16) into plane

p = 0 and q = 1. Eqs. (2.13) reduce to

y = r sin 8 sin 4 ,

z = (r2 -f

2)“2 cos 6 ,

(2.17)

for which the choice given by Eq. (2.15b) is appropriate. Note that the surfaces $J = const. degenerate into plane sections (x/y = const.). We note that Eqs. (2.17) are identical with those given in Eqs. (2.18) of [ 111. Also, Eqs. (2.16) are equivalent to those given in Eqs. (2.13) of [ 1 l] under the transformations x +jz and 19@ 4. Fig. 4 illustrates a progression of ellipsoidal coordinates, starting with a:b:c = 5:l:l (prolate spheroid), through intermediate values 5: 1.1 :l, 5:2: 1, 5:4:1, 5:4.9: 1, and ending with a:b:c = 5:5: 1 (oblate spheroid). Figs. 3 and 4 were drawn using MATHEMATICAD software [25].

3.

Ellipsoidal multipole expansion

We shall now derive the multipole expansion in ellipsoidal coordinates that forms the basis for the infinite element representation to be presented in Section 4. This derivation is modeled after that for spheroidal coordinates given in Holford [lo] and summarized in Burnett and Holford [ill. Let S be the ellipsoid specified by Eq. (2.1), forming an artificial surface surrounding a structure that is vibrating with angular frequency w. Let p(r) elwt be the acoustic field generated in the exterior of S, where p(r) satisfies the Helmholtz equation V2p + k2p = 0,

and the Sommerfeld

radiation

dP

condition

as R+@J,

x+ikp=o(l/R)

with k(=w/c) the acoustic wavenumber, and R = ]rl. We shall call any such function p, which is assumed to be at least twice differentiable on, and exterior to, S, a ‘radiation function’ in the exterior of S. With respect to the ellipsoidal coordinates defined by Eqs. (2.13), for which S is the surface r = a, we can now state the following theorem. EXPANSION

THEOREM.

p(r) can be expanded

Let p(r) be a radiation function

in ellipsoidal

coordinates

in the region exterior to the ellipsoid S(r = a). Then,

(r, 0, 4) in the form

mG,@, 4; k)

p(r) = emikr C n=l

rn

where the series converges

(3-l)

’ for r > a, and converges

r > a + E(E > 0). The series may be differentiated absolutely and uniformly

in r > a i-

E.

absolutely

any number

and untfonnly

in r, 0, and $J in any region

of times, and the resulting series all converge

D.S. Burnett, R.L. Holford I Comput. Methods Appl. Mech. Engrg. 164 (1998) 49-76

56

The proof takes as its starting point the Helmholtz integral representation, which expresses exterior of S in terms of its values, and those of its normal derivative, on S:

the field in the

(3.2) where n’ is the outward normal to S, p = Ir - r’), and du’ is the element of area at position r’ on S. We now wish to introduce the ellipsoidal coordinates given by Eqs. (2.13), which we shall rewrite in the form: z = (2 -f

y = (r2 - g*)“*B,

x=rA,

*y2c,

(3.3)

with A=(~-~c~~*~)“~cos~, B = sin 8 sin 4 , C=cos0(1

(3.4)

-pcos2C$)1’2,

and where, from Eqs. (2.3), (2.4), and (2.14), f*=&c*, with g S$

g2 = a2 -b*,

q=l-p,

p=g21f2,

(3.5)

Note from Eqs. (3.4) that

A*+B*+C*=l.

(3.6)

For a point r’ = (x’, y’, z’) = (a, 0’, 4’) on S, we have also x’ = aA’ ,

y’ = bB’ ,

Z’ = CC’ )

(3.7)

where A’, B’, C’ are given by Eqs. (3.4) with 19,C#Jreplaced by 8’, 4’. In ellipsoidal distance function p = ]r - r’l appearing in Eq. (3.2) is given by p = [r* - 2aAA’r - 2bBB’(r* - g*)“* - 2cCC’(r* -f*)“’

coordinates,

+ d]“*,

therefore, the

(3.8)

where d = a2 - g2(B2 + B’*) -f2(C2

+ C’*),

(3.9)

and where we have made use of Eqs. (3.5) and (3.6). To proceed with the derivation of Eq. (3.1), consider which we shall denote by J(r). Then, we have

the second term on the right-hand

side of Eq. (3.2),

pik(p-r) J(r)

eikr =

II

H(r’) e S

P

dcr’ ,

(3.10)

where H(r’) = (HIT-’ ap(r’) / an is a continuous function. Put now C = 1 /r and consider f to be a complex variable. (This is the key step in the proof, as in [lo] and [ 111.) The kernel of the J integral is s(5) = p-l

- 41

exp]-ik(p

= lPeLi2 exp

-; [

(p”* _ 1)

1

(3.11)

with P = 1 - 2aAA’l-

2bBB’[( 1 - g212)“*

- 2cCC’5( 1 -f212)1’2

+ dl* .

(3.12)

are chosen having the value + 1 at l= 0; and in In Eq. (3.12), those branches of (1 - g*l*)“* and (1 -f’~*)“* Eq. (3.11), that branch of P1’2 is also chosen having the value + 1 at 5 = 0. Now from Eqs. (3.5), we have g “f < a, and it follows that the branch points of P at C = t g - ’ and 4’ = tf -’ lie outside the circle (51 = a-’ in the complex l-plane. Let us assume that any zeros of P also lie outside, or possibly on, the circle ([( = a-’ (see below). Then the branch points of PI” generated by these zeros also lie in the region ([IS a-‘. From this it

D.S. Burnett, R.L. Holford I Comput. Methods Appl. Mech. Engrg. 164 (1998) 49-76

51

follows by inspection of Eq. (3.11) that the function g( 5) defined therein is analytic inside the circle 1[I = a Consequently, it may be expanded in the form of a convergent power series: in )l(
g(5) = n#, s,Y”

,

’.

(3.13)

where the constant term (g,) vanishes since g(0) = 0. The remainder of the proof follows exactly along the lines given in Holford [lo]. Namely, if the above expansion for g(f) is substituted in Eq. (3. lo), the series can be integrated term by term, showing that Jelkr = c a,(& 4)/r” *=I Moreover, this power of convergence, where term on the right-hand same manner. Hence,

in r>a

.

series converges absolutely and uniformly in any circle Ir-‘l s (a + E)-’ inside its circle it can be differentiated with respect to r, 0 and 4 any number of times. Finally, the first side of Eq. (3.2) is similar in structure to J, and may therefore be treated in exactly the on combining the two terms, the expansion given by Eq. (3.1) results.

3.1. The zeros of P What remains to complete the proof of the expansion theorem is to show that the zeros of P lie outside (or on) the circle I[( = a-‘, as we have assumed. This is equivalent to the assertion that the zeros of p*, as a function of r, all lie inside (or on) the circle jr1 = a in the complex r-plane. Putting z = r/a in Eq. (3.8), where z is a complex variable (not to be confused with the Cartesian coordinate z), this is equivalent to the assertion that the roots of the equation (3.14)

z2 - 2AA’z-2bBB’&-&2:CC’JzT--72+6=0 lie in the unit disc lz[ G 1, where we have introduced b=b/a,

Z=cla,

gLg?/&

l-62,

S=d/a2,

the dimensionless

parameters (3.15)

and 72 =f 2/a* = 1 -C*.

(3.16)

In terms of these, S may be written 6 = 1 - a*(B* + Bf2) - 7*(C2 + C’*) .

(3.17)

q=l-p.

(3.18)

Also, p = (T2/T2 )

Eqs. (3.15) through (3.18), together with Eqs. (3.4) and their primed counterparts, show that Eq. (3.14) depends on six independent parameters: a, r, 6, 4, 8’, and 4’. The first two satisfy 0 < u < r < 1, and specify the shape of the ellipsoid S. The remainder range over the values shown in Eq. (2.15a or b), and specify the angular positions of the variable points r and r’. To proceed with the analysis of Eq. (3.14), we shall rewrite it in simpler notation in the form

z2- az-pp-T-y\+8=0,

(3.19)

where (Y= 2AA’,

p=2bBB’,

y = 2CCC’ .

(3.20)

As it stands, Eq. (3.19) is irrational in z, which complicates the analysis of its roots. In the special cases of a prolate spheroid ((T = 7) or an oblate spheroid ((T = 0), Eq. (3.19) can be rationalized with the simple substitution z = G-cash x, yielding a quartic for the variable A = eX. In Holford [lo], it was shown analytically

D.S. Burnett, R.L. Holford I Comput. Methods Appl. Mech. Engrg. 164 (1998) 49-76

58

that the roots of this quartic do indeed satisfy ]z] G 1. In the general case 0 < cr < T < 1, no such simplification possible, and we must proceed differently. Let us introduce the functions 11

is

=z2 -cuz+s,

.=pdZT,

(3.21)

W=#7. Then, Eq. (3.19) is equivalent

to the equation (3.22)

u-v-w=o. This equation

can be rationalized

by replacing

it with the equation

(U - Y - w)(u - V + w)(u + V - w)(u + I/ + w) = 0. This equation

is equivalent

(cd2 + v*-

(3.23)

to

w*>* - 4U2V2 = 0 )

which, in turn, may be written as (U+v-w)2-4Uv=o,

(3.24)

where

u = (z2-

az +

s>2,

v= p’(z2- G-‘>)

(3.25)

w= y2(z2- T2). Eq. (3.24) is seen to be an equation of the eighth degree (i.e. an octic) in z. Before continuing, it should be pointed out that Eq. (3.23) does not introduce any new roots over and above those given by Eq. (3.22). To see this, let us fix the values of the variables 8’ and +‘, and let the ranges of 0 and 4 be as given by Eq. (2.15a). Now, let z = z, be a root of u - v - w = 0 [i.e. Eq. (3.22)] corresponding to specific values of 13and #J lying in the interval (0, n/2). Then, z = z, will be a root of u - v + w = 0 [i.e. the second factor in Eq. (3.23)] if 19 is replaced by rr - 19,since this changes the sign of C in Eqs. (3.4), but not those of A or B-i.e. it changes the sign of y ‘in Eqs. (3.21), but not those of LYor p (or 6). Similarly, z = z0 will be a root of u + v - w = 0 [the third factor in Eq. (3.23)] if 6J is replaced by 2rr - 13;and z = z0 will be a root of u + v + w = 0 [the fourth factor in Eq. (3.23)] if 0 is replaced by IT + 8. Finally, if the roots of Eq. (3.23) lie in the disc ]z] 6 1, so also do the roots of the equation obtained by replacing z by -z therein. From Eqs. (3.24) and (3.25), this is seen to be equivalent to replacing (Y by - LYin Eq. (3.23), or, equivalently, by replacing 4 by IT - 4. Moreover, the same results obtain if we fix 13 and 4 and consider the same set of replacements for 0’ and 4’. The conclusion to be drawn from these observations is that we may restrict our analysis of the roots of Eq. (3.23), or equivalently Eq. (3.24), to the range 0 G 0, &8’, Writing

(3.26)

4’ G n/2.

out Eqs. (3.24) and (3.25), we now formulate

the following

conjecture.

CONJECTURE. For any values of u and r such that 0 ~cr~r* + P2(z2 - (T2) - y2(z2 - T2)]’ - 4P2(z2 - (r2)(z2 - Lyz + q2 = 0 all lie in the unit disc (z( < 1, where (Y, fi, y, 6, u, and r are given by Eqs. (3.20). For the special cases g = r and u = 0, mentioned

above, this conjecture

4,0’,and

+‘, the

(3.27)

(3.4), (3.15) through (3.18), and

was proven

in [lo]. For general

D.S. Burnett,R.L. Holford I Comput. Methods Appl. Mech. Engrg. 164 11998) 49-76

s9

values of u and 7, we have no formal analytical proof, but we have gathered strong numerical evidence in support of it. That is, for selected values of u and r, we have evaluated the roots of Eq. (3.27) for discrete values of 0, 4, 0’ and 4’, using MATHEMATICAB software [25], and in all cases the conjecture has been found to be true. The cases examined ranged over values of a:b:c varying from 3:2:1 to 1OO:lO:l. The results of a representative calculation are shown in Table 1. Here, the generating ellipsoid S has its axes in the proportion a:b:c = 10:3:1, implying IY*= 0.91 and r2 = 0.99. The first four columns of the table give the Table1 Minimum and maximum valueof themoduliof theroots ofJ?q. (3.27) fora:b:c= 10:3:1 H

8’

6

d’

0 0

0

0

0

0

30

0

0

0

60

0

0

0

90

0

0

30

0

0

0

0

0

0

0

0

0

0

0

30

0

30

0

30

0

30

0

30

0

30

30 30 30 60 60 90 0 0 0 0 30 30 30 60 60 90 0 0 0 0 30 30 30 60 60 90 0 0 0 0 30 30 30 60 60 90 0 0 0 0 30 30 30 60 60 90

0

0

30

0

30

0

30

0

30

0

60

0

60

0

60

0

60

0

60

0

60

0

60

0

60

0

60

0

60

0

90

0

90

0

90

0

90

0

90

0

90

0

90

0

90

0

90

0

90

30 30 30 30 30 30 30 30 30 30

30 30 30 30 30 30 30 30 30 30

60 90 60 90 90 0 30 60 90 30 60 90 60 90 90 0 30 60 90 30 60 90 60 90 90 0 30 60 90 30 60 90 60 90 90 0 30 60 90 30 60 90 60 90 90

Min 0.840 0.642 0.232 0.300 0.390 0.135 0.564 0.547 0.879 1.000 0.861 0.672 0.267 0.265 0.413 0.112 0.542 0.529 0.865 0.988 0.910 0.746 0.362 0.173 0.463 0.059 0.494 0.490 0.836 0.962 0.959 0.832 0.487 0.100 0.491 0.012 0.466 0.469 0.820 0.949 0.880 0.703 0.308 0.228 0.436 0.038 0.508 0.500 0.844 0.962

Max

0

8’

4

4’

1.000

30

60

0

0

0.955

30

60

0

30

0.713

30

60

0

60

0.300

30

60

0

90

1.000 0.888 0.564 1.000 0.879

30

60

30

30

60

30

60

30

60

30

60

I .ooo

30

60

0.999 0.943 0.685 0.265 0.990 0.873 0.542 0.988 0.865 0.988 0.990 0.902 0.605 0.173 0.964 0.838 0.494 0.962 0.836 0.962 0.959 0.832 0.487 0.100 0.947 0.818 0.466 0.949 0.820 0.949 1.000 0.929 0.654 0.228 1.000 0.876 0.548 1.000 0.870 1.000

30

90

30

90

30

90

30

90

30

90

30

90

30 30 30 60 60 90 0 0 0 0 30 30 30 60 60 90 0 0 0 0 30 30 30 60 60 90 0 0 0 0 30 30 30 60 60 90 0 0 0 0 30 30 30 60 60 90

30

90

30

90

30

90

30

90

60

60

60

60

60

60

60

60

60

60

60

60

60

60

60

60

60

60

60

60

60

90

60

90

60

90

60

90

60

90

60

90

60

90

60

90

60

90

60

90

90

90

90

90

90

90

90

90

90

90

90

90

90

90

90

90

90

90

90

90

60 90 60 90 90 0 30 60 90 30 60 90 60 90 90 0 30 60 90 30 60 90 60 90 90 0 30 60 90 30 60 90 60 90 90 0 30

60 90 30 60 90 60 90 90

Min

Max

0.924 0.781 0.417 0.163 0.482 0.080 0.480 0.482 0.829 0.949 0.970 0.844 0.508 0.173 0.510 0.153 0.494 0.495 0.836 0.962 0.960 0.838 0.502 0.228 0.524 0.201 0.508 0.510 0.844 0.962 0.990 0.868 0.545 0.265 0.546 0.261 0.542 0.542 0.865 0.988 1.000 0.879 0.564 0.300 0.564 0.300 0.564 0.564 0.879

0.996 0.887 0.572 0.163 0.989 0.859 0.527 0.988 0.856 0.988 0.970 0.844 0.508 0.173 0.962 0.835 0.494 0.962 0.836 0.962 1.000 0.875 0.553 0.228 1.000 0.871 0.548

I .ooo

1.ooo

I .ooo 0.870 1.000 0.990 0.868 0.545 0.265 0.988 0.865 0.542 0.988 0.865 0.988

1.ooo 0.879 0.564 0.300

1.000 0.879 0.564

1.ooo 0.879

D.S. Bumerr, R.L. Holford I Comput. Methods Appl. Mech. Engrg. 164 (1998) 49-76

60

values of the angular variables in degrees, and the last two columns give the minimum and maximum values of the moduli of the roots of Eq. (3.27), accurate to three decimal places. In all cases, we see that the maximum modulus does not exceed unity, as conjectured. Note that we have restricted the ranges of the angular variables as in Eq. (3.26), together with 0 6 8’ and 4 s 4’. The reason for this latter restriction is that the coefficients (Y, /?, y and 6 are each symmetrical in 0, 8’- and 4, $‘. Two observations concerning the results shown in Table 1 are in order. The first has to do with conditions under which the maximum modulus achieves the value unity. In this regard, note from Eq. (3.19) that z = 1 is a root thereof if l-a-p&G+6=0, or, from Eqs. (3.17) and (3.20), 1 - 2AA’ - 2b*BB’ - 2c*CC’ + 1 - (1 -b’)(B* Making

use of Eq. (3.6) and its primed counterpart,

+ B’*) - (1 -c2)(C2 this condition

+ C’2) = 0.

may be written

(A - A’)* +b*(B - B’)2 +c’(C - C’)2 = 0, which implies A = A’, B = B’ and C = C’. These, in turn, imply 0 = 8’ and $ = 4’. This condition is seen to obtain for all entries with Max = 1.OOOin Table 1. It is readily verified that z = 1 is, in fact, a double root of Eq. (3.19) when 0 = 8’, 4 = 4’. The second observation relates to entries in the table for which Min = Max; i.e. the roots all have the same modulus. This is seen to occur in the following cases: (i) 4’ = 90” and either 8 = 0 or 4 = 0 (implying (Y= p = 0); (ii) 4’ = 90” and 8 = 90” (implying (Y= y = 0); and (iii) 0 ’ = 90” and 4 = 0 (implying p = y = 0). In each case, two of the coefficients (Y, p and y vanish, and the octic equation then degenerates. To elaborate, consider case (i) with 4’ = 90” and /3 = 0, for which (Y= 0, p = 0, and y = 2Z( 1 - p cos2@“2

cos 8’ , (3.28)

6 = 1 - a2 sin*0’ - r2( 1 -p

cos’+ f cos28’) .

Eq. (3.27) reduces to [(z2 + a)2 - y2(z2 - r2)]’ = 0. Putting

1+4 = z2, the roots of this equation

(3.29) are obtained

+‘-~~+rn=O

by solving

the quadratic

equation (3.30)

with f=

m = y27*

y2-28,

+ a*.

(3.31)

It is readily shown that the discriminant A = e2/4 - m of Eq. (3.30) is negative. Explicitly, from Eqs. (3.28) and (3.31), A = - y*[( 1 - c~* sin*4) sin’8’ + p cos’r$ cos*@‘]. Hence, Eq. (3.30) has complex conjugate roots given by $zi+i(m-$)‘I*, for which we have Ifi)’ = m. That is, the roots of Eq. (3.29) all have the same modulus, given by ]z] = m”4. Moreover, we can show that m S 1. For, putting 5 = sin’8 in Eqs. (3.28) and (3.31), we have m = (s - t)‘&*

-

2(ks + kt - 2st)& + k2,

where s = 1 - g2, t = 1 - 72 and k = 1 - cr2 COS’~ are such that 0 G t < s 6 k s 1. It is readily shown that the graph of m vs. 5 has its minimum in 5 2 1. It therefore follows that the maximum value of m in 0 G 5 G 1 occurs at 5 = 0, where m = k2 G 1. A similar analysis holds in case (ii). Case (iii) is an even more degenerate case. With p = y = 0, Eq. (3.27) reduces to

D.S. Burnett, R.L. Holford

I Comput. Methods Appl. Mech. Engrg.

164 (1998) 49-76

61

(22-(YZ+S)4=0.

(3.32)

Again, it can be shown that the discriminant (a’/4 - 8) of the quadratic factor is negative, implying that the roots of Eq. (3.32) all have modulus ]z] = 8”2. Moreover, we have explicitly S = 1 - (T’ sin2+’ - 7’9 cos28 in this case (i.e. for 8’ = 90”, 4 = 0), from which it is evident that 6 G 1 for all 8 and 4’. (The condition p = y = 0 also arises in case 0 = 0 and 8’ = 90”. In this case, however, it is found that the discriminant may be positive or negative, depending on the values of 4 and +‘.) We have dwelt on these observations for the reason that they are general results that obtain for special cases of the angular variables. For general values of these variables, although methods are available for evaluating the number of roots of a polynomial equation lying inside the unit disc [26,27], the algebra involved is formidable. Until a formal proof of our conjecture is forthcoming, therefore, we are forced to rely on the numerical evidence cited in support of it. Before leaving this discussion, we point out that the expressions for A, B and C given in Eqs. (3.4) can be simplified by introducing new angular variables, 0 and @, according to A = sin 0 cos @ ,

(3.33)

B = sin 0 sin (D ,

c = cos 0, where 0 G 0 < rr, 0 G @ G 21r. That is, there is a l-l transformation between (8, 4) and (0, @) defined by Eqs. (3.4) and (3.33). These angular variables are similar to those recently employed by Kallivokas et al. [28], for example, in their development of an absorbing boundary condition for an ellipsoidal truncation boundary. Note that, except in the limiting case of an oblate spheroid, when 0 = B and @ = 4, the coordinates (I, 0, @) do not form a mutually orthogonal system of ellipsoidal coordinates as (r, 13,4) do.

4. Ellipsoidal infinite element 4.1. Global geometry:

shape and orientation

of elements

The infinite element is shown in Fig. 5. One face of the element, the ‘base’, must attach to, and conform to the shape of, an ellipsoidal surface of ellipsoidal radius rl (previously denoted by a in Sections 2 and 3) surrounding the structure, called the ‘infinite element ellipsoid’. The base (Fig. 6) may be a curvilinear quadrilateral, as shown here, or triangle; it need not conform to the 0, 4 coordinate lines. This permits one to construct a completely general 2-D mesh on the ellipsoid comprising curvilinear quadrilaterals and triangles of any shape and orientation. Any type of FE representation may be used over the base, e.g. Lagrange, serendipity or hierarchic polynomials of any degree. (The quadratic Lagrange nodal pattern shown in Figs. 5 and 6 and later in Fig. 7, is for illustrative purposes.) The infinite element ellipsoid uniquely determines three pairs of foci, which in turn define a system of

INFINITE

3ASEOFELEMErn

Fig. 5. The ellipsoidal

infinite element.

Fig. 6. The base of the infinite element on the infinite element spheroid.

Dashed lines are 0, 4 coordinate

lines.

D.S. Burnett, R.L. Holford I Comput. Methods Appl. Mech. Engrg. 164 (1998) 49-76

62

Fig. 7. Illustration of local node numbering convention for the particular case of a quadratic directions and a quadrupole (m = 3) in the radial direction (N = m X n = 27).

Lagrange

quadrilateral

(n = 9) in the angular

confocal ellipsoids and hyperboloids (cf. Fig. 3) that are used in the construction of the rest of the element. The side faces of the element are the loci of confocal hyperbolas emanating from the sides of the base. A multipole of order m[p = (a, /r + * -. + a,/r”) e -ikr] requires m layers of nodes that are on confocal ellipsoids of ellipsoidal radii r, , r2, . . . , r,,,. The nodal pattern is identical on each ellipsoid. The value m = 2 corresponds to a dipole, m = 3 to a quadrupole (as shown in Fig. 5), m = 4 to an octupole, etc. Finally, there is an outer ellipsoidal face, $ at ellipsoidal radius ?, that recedes to infinity in the theoretical development. Conventional 5,~ coordinates (illustrated in Fig. 6 for a quadrilateral) are defined over the element cross section, i.e., on ellipsoids, by the mapping

where n is the number of nodes in the quadrilateral or triangle, 8, and +,, are the ellipsoidal angular coordinates of the vth node, and ,Y”,< 5, Q) are interpolation polynomials. (Alternatively, blending functions could be used if the elements are hierarchic.) Since the base, intermediate nodal layers, and outer face conform to confocal ellipsoids and the side faces are the loci of confocal hyperbolas, the element is a right cylinder in r, 0, &space (or r, 5, g-space). If the geometry for the finite-size elements is generated using a conventional polynomial mapping (rather than blending functions), there will be a discontinuity (in 5,~ coordinates and therefore in the dependent variable, pressure) between the infinite elements and adjacent finite-size elements that are on the inside of the infinite element ellipsoid. The finite-size elements will not conform exactly to the shape of the ellipsoid, whereas the infinite elements do conform exactly. This does not cause a problem because the error due to the geometry discontinuity converges to zero, as the mesh is h or p-refined, at the same rate as does the (well-known) error at the boundary of a curved domain when using conventional polynomial-mapped elements. Thus, there is asymptotic continuity in pressure across the infinite element ellipsoid. 4.2. Physics 4.2.1. Governing physics equations To recapitulate, the time-harmonic (e’“‘) behavior of the acoustic field is governed by the 3-D Helmholtz equation, V;?, + k2p = 0

(4.2)

where k is the wavenumber (‘w/c), c is sound speed (=mp), B is bulk modulus, p is density, and p is the complex-valued amplitude of scattered and/ or radiated pressure:

D.S. Burnett, R.L. Holford I Comput. Methods Appl. Mech. Engrg. 164 (1998) 49-76

c. scat

P

P=

P

‘p total_p’“”

rad=

scat+rad

if scattering, if radiation,

P’“‘“’ =

total

P

(4.3)

-P i”c if both .

To ensure uniqueness of the solution, ‘boundary’ at infinity: $-+ikp=o(l/r),

63

the pressure must satisfy the Sommerfeld

radiation

condition

r--+w.

at the outer

(4.4)

The r in Eq. (4.4) is a spherical r. However, ellipsoidal radii approach spherical radii as r -+ m, so that Fq. (4.4) can be used with ellipsoidal coordinates in the limit as r-+w. 4.2.2. Finite element representation of pressure The scattered and/or radiated pressure in the infinite

element

is represented

as follows:

m GJ5,rl;k) ~(5,

rl, r) = eeikr

C &L=I

rF

(4.5)

where (4.6) This is the multipole expansion of Section 3 [see Eq. (3.1) et seq.] modified in two ways: (i) the infinite series is truncated at m terms (cf. [ 1,l l] for practical values of m); and (ii) each of the m unknown C” functions G,(O, 4; k) of Section 3 is approximated by GL( 6,~; k), which are expansions in conventional polynomial shape functions, #t( 5, v), with 5, v being local angular variables related to 13,4 by the coordinate mapping in Eqs. (4.1). In the following development we begin with a more general expression, and at the end of the development it will be shown how the resulting expression is in the form of Eqs. (4.5) and (4.6). General expression; DOF numbering We begin with the expression

where $(k

rl, r) = KC& 59$;(r)

v=l,2

,...,

n;p=l,2

,...,

m;

nXm=N

(4.8)

Here, $“,( 5,~) are ‘angular’ shape functions [the same as in Eq. (4.6)] that interpolate p over ellipsoidal surfaces confocal to the infinite element surface, and (//L(r) are ‘radial’ shape functions that interpolate p along hyperbolic rays. (Hierarchic functions could also be used; see comment below.) Node numbering begins with the nodes on the base of the element and proceeds radially outward, a layer of nodes at a time. This is summarized in Table 2 and illustrated in Fig. 7. Angular shape functions The functions I/J”,<&,v) are conventional 2-D interpolation-polynomial serendipity, which provide Co continuity between the sides of adjacent Radial shape functions The function @L(r) uses a truncated mth order multipole expansion:

shape functions, infinite elements.

form of the radial part of the multipole

expansion

e.g. Lagrange

or

of Section 2, viz. an

D.S. Burnett, R.L. Holford I Comput. Methods Appl. Mech. Engrg. 164 (1998) 49-76

64 Table 2 Note numbering

pattern relating nodal (degree of freedom)

On ellipsoidal surface at r, (base of infinite element)

surface at r2

On ellipsoidal

surface at rm

On ellipsoidal

s;(r)

=

e-ik(r-rII)

i

index j to radial index p and angular

j

CL

v

1 2

1 1

1 2

n

;

n

nfl ni2

2 2

1 2

in

;

n

(m-l)n+l (m-l)n+2

m m

2

N

m

n

hWp’

/L’=I (/try’

p=l,2

,...,

index v, N = m X n

1

m

(ma2).

(4.9)

The phase factor eikrp does not need to be included; if omitted, it will simply appear in each h,,, in Eq. (4.12) below. The factors k“’ m ’ the denominators are also not necessary; they are included only to make the hWLP8 nondimensional. Eq. (4.9) represents m mth degree polynomials in the reciprocal variable l/r in the finite interval 1/r E (0, 1 /r,], with the constant terms omitted because t&(r) + 0 as 1 /r +O. Hence, there are m coefficients h ppf in each polynomial, which are determined by two considerations: (i) the requirement, for convergence, of asymptotic Co continuity (cf. end of Section 4.1) between the base of the infinite element and adjacent finite-size elements on the inside of the infinite element ellipsoid; and (ii) the desirability (not necessity) of having the degrees of freedom represent nodal quantities, for ease of postprocessing interpretation. Both features are achieved by transforming the polynomials to Lagrange interpolation polynomials in 1 /r with m interpolation points at 1 /r,, . . . . , 1 lr,, using the interpolation property (4.10)

$L(r,,) = aWLc(, . Applying Eq. (4.10) to Eq. (4.9) yields m sets of m linear algebraic equations: hS=Z

(4.11)

where S, = (krj)-’ and Z is the identity matrix. Therefore, h=S-‘.

(4.12)

This procedure defines m layers of nodes (with n nodes on each layer) lying on ellipsoids of ellipsoidal radii rr, r2,. . . rm (cf. Fig. 7). The resulting Lagrange polynomials in 1 /r are, effect, ‘interpolation multipoles’ in r over a semi-infinite interval. To illustrate, for a dipole element (m = 2) t&(r) = e- ik(r-rJ(!$+-f%)

p = 1,2

(4.13)

and h=A

-krf r2 - rl

kri

k2rir2 -k2r,ri

1’

(4.14)

D.S. Burnett, R.L. Holford I Comput. Methods Appl. Mech. Engrg. 164 (1998) 49-76

65

COMMENTS

(1) It is easily shown that the final expression for p( &T, r) is in the form of Eq. (4.5). Substituting into Eqs. (4.7) and (4.8) yields a linear transformation between Q,(k) and P,

Eq. (4.9)

(4.15) where (4.16) and P, and P “+, are identical lists of nodal pressures, the former using a single index and the latter using two indices, as explained in Table 2. (2) The above development employs interpolatory shape functions. For a hierarchic formulation, the angular directions would employ the standard 2-D hierarchic shape functions for quadrilaterals and triangles [29]. The radial direction would use the mapping c = l -2r, /r, which linearly maps the interval 1 /r E (0, 1 / Y,] to the interval ,!JE [- 1, I), and would then employ the standard 1-D hierarchic shape functions in ,g+’ [29], excluding the linear function that is unity at infinity (i.e. at i’= 1). 4.3. Differential 4.3.1.

geometry

of ellipsoidal

coordinates

Metric tensor

To simplify

the analysis

throughout

the remainder

of Section 4, we rewrite Eqs. (2.13):

x = rQ cos 4 y = rG sin 0 sin 4

(4.17)

z = rF cos 0 P where

P=(l

-pcos2+)“2

Q = (1 - 4 ~0~~8)“~

(4.18)

L!

p’J-

f’

q=l-p

f .Y= - = the maximum r

eccentricity

of all three ellipsoidal

The quantities F, G, P, Q, p, q, E are all dimensionless The metric tensor, g,,, is g a@=a a ‘Up

.

with values in the closed unit interval

a, P = rr @,+

where a, are the ellipsoidal ax

cross sections.

ay .

covariant

[O,l].

(4.19) base vectors, given by

az (4.20)

and i, j, k are the unit Cartesian

base vectors. Using Eqs. (4.18), we have

D.S. Burnett, R.L. Holford I Comput. Methods Appl. Mech. Engrg. 164 (1998) 49-76

66

a,=Qcos$i+G-’

sin8sin4j+F-‘Pcos$k

a, = rqQ_’ sinOcosOcos~i+rGcos$sin~j-rFPsinOk

(4.21)

a,=-rQsin4i+rGsin0cos4j+rFpP-‘cosOsin4cos4k It follows immediately g@

=

0

from Eqs. (4.21) that the off-diagonal

terms of g,,

are zero: (4.22)

ff#P

i.e. the ellipsoidal

coordinates

are orthogonal.

The scale factors, h,, are determined

by the diagonal

terms:

0 = r, l3,q75.

k, = Ia,1 =vL

These are a bit tedious to evaluate.

(4.23)

To illustrate

the procedure,

hf = a, *a, = Q2 COS’C#J + G-’ sin20 sin’+

consider

the case (Y= r.

f Fb2P2 cos213.

(4.24)

Multiply both sides by F2G2 to create a polynomial in even negative powers of r. Then trigonometric expressions in the coefficients of each power. The final results are as follows:

simplify

the

h = (1 - E’Q=)“*[~ - ~‘(1 -P*)]“’ r FG h = r(1 - E*Q*)“~(P~ + Q* - l)“* (4.25)

Q

0

h = r[l - ~‘(1 - P2)]“2(P2 + Q’ - 1)“’ 4 P 4.3.2. Jacobians Volume Jacobian J, = a, *a, X a+ = h,e, . hoe0 X h,e+

(4.26)

where e,, e,, e4 are unit vectors in the r, f3, 4 directions, respectively. These ellipsoidal right-handed orthogonal system, so e, . e, X e, = + 1. Hence, using Eqs. (4.25),

J, = h,h,h,

= r

2 (1 - &‘Q’)[l - ~‘(1 - P2)](P2 + Q’ - 1) FGPQ

We observe that Jv reduces to the two spheroidal l

coordinates

are a

(4.27)

forms [ 1,111:

Prolate spheroidal f=g,p=l,q=O,P=sin4,Q=l,whichyields J, = (r2 - f2 cos2t$) sin 4

(4.28)

where I#Jis the polar angle from the x-axis, which is the major axis. l

Oblate spheroid g=O,p=O,q=l,P=l,Q=sin& J = r(r2 -f V

(r2

which yields

* sin’0) sin 0 _f2)“2

where 13 is the polar angle from the z-axis, which is the minor axis.

(4.29)

D.S. Burnett, R.L. Hoiford / Comput. Methods Appl. Mech. Engrg.

164 (1998) 49-76

67

Surface Jacobian J, = la, x a41 = h,h, = r2

(1 - E2Q2)“*[l

- F2[1 - F2)]“*(P2

+ Q2 - 1) (4.30)

PQ

J, will be needed later on the outer face of the element as the face recedes to infinity (while f remains constant). Using lim,,, E = 0 yields (4.31) which, as with J,, reduces to the two spheroidal

forms [I, 111.

4.4. St@ness, mass, and radiation-damping matrices 4.4.1. Formal expressions for integrals The element matrices are derived by starting with a finite-size element, i.e. with the outer face on an ellipsoid of ellipsoidal radius i (cf. Fig. 5), and then taking the limit as i + 03. Thus, applying the Bubnov-Gale&in weighted residual method [30] to Eq. (4.2) over a single element yields lim i+=

(SV*p -t- 02pp)$

dV= 0

(4.32)

i=l,2,...,N

using k2 = “‘p/B. The first integral is converted into a surface integral and another volume integral using the Eq. (4.7) into the volume identity @,V’p = V. (eiVp) - V+. -VP and the divergence theorem. Substituting integrals yields the following element matrix equation: (4.33)

(K - w2kf)P = F where the stiffness Kij =lim i-t= Fi =lim i+r

matrix, K, mass matrix, M, and pressure BV$i - Ve dV

gradient

vector, F, are, respectively,

Mij =lim F-P

aP B+i ,dS. ur1

(4.34)

The surface integra!( F F is over the entire boundary of the element, Se). It is split into two integrals: “(Cl over the outer face, S e , and the other over the remaining faces, Str) - S , F, =lim i--+X

aP B&zdS.

one

(4.35)

Consider the first integral in Eq. (4.35). It can be shown that as i + 00, ellipsoids approach spheres, ellipsoidal radii approach spherical radii, dp/dn 4 dp/&-, and from Eq. (4.31), dS +r2((P2 + Q2 - l)/PQ) dt? d#, where r is spherical or ellipsoidal (they are identical in the limit), B and 4 are ellipsoidal angles, and P and Q are functions of ti and # but not of r [cf. Eqs. (4.18)]. To evaluate i@/dr as i -+a, substitute Eqs. (4.7), (4.8) and (4.9) into Eq. (4.4), which yields a stronger form of the Sommerfeld condition, (4.36) where the upper-case 0 means ‘at least as fast as.’ Substitute dp/ &- in Eq. (4.36) for aplan in the first Integral. The 0( 1 /r2) term makes no contribution to the integral because &. is 0( 1 /r) and, as noted above, dsmr’ d8 d+ as i + m. In the remaining term ikp, substitute Eq. (4.7) for p. In the second integral in Eq. (4.35), substitute the balance of linear momentum for dp/&z, i.e. @/an = w2pu,, where u, is the amplitude of the normal component of particle displacement. Making these substitutions, Eq. (4.35) becomes F=-iwCP+D,

(4.37)

D.S. Burnett, R.L. Holford I Comput. Methods Appl. Mech. Engrg. 164 (1998) 49-76

68

where C, = pc lim

i-m II s”, *i# dS 7

Di = (wpc)* lim i-F

(4.38)

ku,, dS .

The C matrix is the radiation-damping matrix, representing radiation loss to infinity. The D vector permits specification of u, on the side or bottom faces of the element. These faces are either on the boundary of, or interior to, the domain. Face on boundary of domain: A nonzero value corresponds to a moving boundary with a known normal displacement (or velocity) amplitude. Along the side faces, a nonzero u, must decay fast enough in the radial direction for the integral in the second of Eqs. (4.38) to be finite. However, the usual condition in practice is a zero value, corresponding to a rigid boundary (e.g. a plane of symmetry or an approximation to the ground or a large baffle). Face in interior of domain: Every interior face is common to two adjacent infinite or finite-size acoustic elements. Assembly of the D vectors from two adjacent elements yields, on the interface, a sum of two surface integrals with opposite sign because the outward normals point in opposite directions [31]. The pair of integrals therefore represents the difference of u, across the interface. A nonzero difference corresponds to a jump discontinuity in the normal displacement, i.e. a separation of the fluid. A zero difference corresponds to preservation of material continuity, the usual condition in practical applications. Substituting Eq. (4.37) into Eq. (4.33) yields for the element matrix equation: (K + ioC - w2A4)P = D

(4.39)

where the stiffness, mass and radiation-damping matrices are given by Eqs. (4.34) and (4.38). The remainder of Section 4.4 will transform the integrals in these three matrices to expressions that can be numerically evaluated. (The Di integrals will not be considered further since, as just argued, they are zero in virtually all practical applications.) 4.4.2. Transformation of integrals to ellipsoidal coordinates The stiffness, mass and radiation-damping matrices in Eqs. (4.34) and (4.38) may be written

Mij = lim

i-G

~ti$J, dr de W

(4.40)

where Jv and J, are given by Eqs. (4.27) and (4.30). The gradient operator in general curvilinear coordinates is (4.41) It is needed only in K,, where it occurs as VI,$.Vt,$J,. Combining the first of Eqs. (4.27) with Eq. (4.41) yields (4.42) Using Eqs. (4.25), we have

D.S. Burnett, R.L. Holford I Comput. Methods Appl. Mech. Engrg. 164 (1998) 49-76

c2 Q(1

Vt+$+‘t&l, = r2FG

P 1 P

+

- P')

E’

a*.

-L_-L

[ FGTj-FGPQ 1 a4

a+. a+

69

1--L--L a+. a+.

a0 a0

(4.43)

where terms have been separated so that each monomial coefficient of a pair of derivatives is the product of a radial function and an angular function. Eqs. (4.43), (4.27) and (4.31) can now be substituted into K,, M, and Cij, respectively, in Eqs. (4.40). In addition, Eq. (4.8) is substituted for I,+;and $, introducing primes to distinguish the indices v and ,u that are associated with the index i:

Performing these substitutions, we then make two observations: (i) all the integrands are sums of terms, each of which is a product of an angular function and a radial function [the angular shape functions are implicitly functions of 0, 4 because of the coordinate mapping in Eqs. (4.1)]; and (ii) since the base of the element conforms to an ellipsoid and the side faces are the loci of confocal hyperbolas, the element is a right cylinder in r, 0, &space (or r, 5, v-space). Therefore, the 3-D volume integrals for K, and Mij separate into products of 2-D angular integrals, A:!,, and 1-D radial integrals, Rx!,; and the 2-D surface integral for C, separates into a product of a 2-D angular integral and a radial function. The resulting expressions are as follows:

(4.45)

where _ A (1) v’v _ A (2) “‘V= A (3) Y’Y

(4.46) A

(4) V’”

_ -

= _ A (5) “‘Y -

(4.47)

70

D.S. Burnett, R.L. Holford I Comput. Methods Appl. Mech. Engrg. 164 (1998) 49-76

(4.48) and a’“’ is the ‘ellipsoidal cross section’ of the element, i.e. any confocal ellipsoidal surface inside the element that is bounded by the side faces. There is only one ellipsoidal cross section in r, 0, +-space because the 0, 4 coordinates of the boundary of the cross section are independent of r. (Thus, S^@)is equivalent to rr@) in the surface integration for C,.) The symbols Rz.‘, and Rz?fi include an overbar because a term will later be removed from each integral. The remaining terms will be labeled R’,‘!, and RzJp, so that all four radial integrals will have uniform symbols (no overbars) in the final expressions. The tensor products A,, yRfi,/r are N X N matrices that are constructed by multiplying each term in the m X m R matrix by the entire n X n A matrix. 4.4.3. Angular integrals The angular integrals in Eqs. (4.46) may be evaluated numerically by transforming them to the local 5, r]-coordinates defined in Eqs. (4.1). Thus, d0 d$ = J dc dv, where the Jacobian is

a8 a4 --

J=

a.5 a-6

--a0 a7

(4.49)

a4 87)

and the transformation is for quadrilateral cross sections, (4.50) for triangular cross sections . All five integrals can be numerically integrated using Gaussian rules since the integrands are smooth and bounded. Thus, if the three ellipsoidal axes are all unequal, then 0 < P < 1 and 0 < Q < 1. If two or three axes are equal, viz. spheroids and spheres, the integrands are still smooth and bounded [ 111. The derivatives of fi,“, are evaluated in the conventional manner,

,

(4.51)

4.4.4. Radial integrals l

R’,L.:,

From Eq. (4.9),

W’

L=ke dr

m+l -ik(r-rp) Iz

a /“a

(4.52)

a=~ W”

where aFe = -ih,,

-((Y - l)h,,,_,

p = 1,2,. . . , m a=1,2,...,m+l

(4.53)

and hpo = h, m+, = 0. Also (4.54)

D.S. Burnett,

R.L. Holford

I Comput.

Methods

Appl.

Mech.

Engrg.

164 (1998)

49-76

71

where b, = 2

(4.55)

aWfya,,p-Y

y=l

and where, from Eq. (4.53), a,, = 0 for cy > m + 1. Substituting Eq. (4.54) into the first of Eqs. (4.47) yields FGe -i2krkdy + $ bp+* p=1

(4.56)

where LP’fL =e

ikCrpj+rg)

(4.57)

lk.

A factor k is introduced into the integrals (and cancelled by the 1 /k in L,,,) to make them dimensionless. The is separated from the others because the upper limit yields an undefined term (which will later be cancelled by similar terms from Mij and C,j). Thus, b, integral

e-‘2k’k dr . ’ (FG - 1) e-i2k’k dr + lim lim FGe -i2krk dr = lim i-x 1;r, r’+x I i-, i-_)=lir,

(4.58)

on the right-hand side of Eq. (4.58) is well defined because - 1 < FG - 1 < 0 for as r + ~0, so the limiting process only requires changing the upper limit from i to ~0. The second integral may be written, formally, The first integral

rI s r < m, and FG - 1 = 0( 1 /r*)

lim e i-+= 1;II

-12krk

dr

=lim

t

e->2ki

_

i

e-i2kr,

.

(4.59)

i-tm

The other integrals in Eq. (4.56), following the summation, are well defined, so again the limiting requires changing the upper limit from i to ~0. Eq. (4.56) becomes, finally,

process only

(4.60) where

-i2’ +

2

ba+21a

p=0

1 ,

(4.62)

5= kr, , _i2k

and a60 is the Kronecker

(4.61)

‘kdr,

p>O

(4.63)

Delta (= 1 if /3 = 0; =0 if p > 0).

. RF?, From Eq. (4.9) -i2kr

k

(4.64)

where P-1 cp=

c h,?%.P-, y=l

(4.65)

and h,, = 0 for cy > m. Substituting Eq. (4.64) into the second of Eqs. (4.47) and letting the upper limit be ~0, since all integrals are well defined, yields

72

D.S. Burnett, R.L. Holford I Comput. Methods Appl. Mech. Engrg. 164 (1998) 49-76

(4.66) where p 2 2.

Jp =

(4.67)

Substituting Eq. (4.64) into the third of Eqs. (4.47) and letting the upper limit be ~0, since all integrals well defined, yields R ;!, = L,&’

2

c~J~+~.

are

(4.68)

p=2

Substituting Eq. (4.64) into the fourth of Eqs. (4.47) and proceeding in the same manner separating the c2 integral from the others and subtracting and adding unity to the integrand) -c41

R P’fi

=

(4) R P”F

=L

L

as for RE?@ (viz. yields

_ 1 _c,i lim e-i*ki + RF)p ILIP l2 2 i-+m

where

1

(4.70)

/3 2 0.

(4.71)

Zm-2 -j2’

@‘IL

+

c~+~J~

c

p=o

Jp =

The second integral is repeated from Eq. (4.67) but the range of /3 now includes p = 0 and 1. For p = 0, the integral is well defined because the integrand + emi”“/ r2 as T+W. For /3 = 1, the integral is also well defined because the integrand + e-i2kr/ r as r+w, which is the integrand of the well defined exponential integral. 4.4.5. Radial function on outer face The radiation-damping matrix, Cij, in Eqs. (4.45) includes

the term lim;,,

[$L,, I&]$~.

Using Eq. (4.64),

(4.72) Taking the limit as i + CCeliminates lim [~&,*;]~f’

= L,.,

2

all the cp terms except c2, which contains

f+i e-i2ki

an undefined

oscillatory

term: (4.73)

i+m

4.4.6. Vanishing of undejned oscillatory terms The above expressions for all the radial integrals and the radial function on the outer face can now be substituted into the stiffness, mass and damping matrices in Eqs. (4.45). Adding those three matrices together as in the element matrix equation, Eq. (4.39), yields K + iwC - w2M = K” - w2Ma + undefined

oscillatory

terms

(4.74)

where K” and M” comprise all the well-defined integrals and are summarized in the next section. The undefined oscillatory terms comprise all terms containing the expression lim,,, e-i2ki. There are three such occurrences: Eqs. (4.60), (4.69) and (4.73), which contribute to K, M and C, respectively. Thus,

D.S. Burnett, R.L. Holford / Comput. Methods Appl. Mech. Engrg. 164 (1998) 49- 76

undefined

oscillatory

terms = BAI!‘&

I~ %

lim e-i2ki ?+a

+ iopcA$‘&,,

?

- ~~pr~A~.‘~~,,

lim emizki i+= 1 c,i - lim edizki l’ 2 i+m

=(b2 + c2)B i Ay,‘$fi,, using w = ck and c2 = B/p. c2 = h,,,h,, = -b2. Hence, undefined

oscillatory

From Eqs. (4.53)

73

!im e-i2ki ,+=

and (4X),

(4.75)

b, = aF.,aP, = -h,,,h,,,

and from Eq. (4.62),

terms = 0 .

(4.76)

4.4.7. Final expressions for element matrices From Eqs. (4.74) and (4.76), the final form of the element matrix equation is (K” - o”M==)P = D,

for the ellipsoidal

infinite element

(4.77)

where

and the indices i and j are related to Y’, /_L’and V, p, respectively, as in Eqs. (4.8), (4.44) and Table 2. The vector D, cf. Eqs. (4.38), is zero in virtually all practical applications. The angular integrals, AT!,, i = 1,2, 3,4,5, are given in Eqs. (4.46). The integration over @,C$is transformed to 5, 7 using Eqs. (4.49)-(4.51). The radial integrals are as follows: + 2 bP+JP p=o

(4.79)

The quantities E,, J, LPrr, b,, cp, a_, and h,., are given by Eqs. (4.48), (4.62), (4.57), (4.55), (4.65), (4.53) and (4.12), respectively. The Zp and JP integrals are given by Eqs. (4.63) and (4.71), respectively. They are not independent, but are related by the expression Zp = JP - (1 +p)&;JP+2

+p[4~;JP+4,

p >O.

(4.80)

The range of /3 for Ip in Eqs. (4.79) is 0 S j3 S 2m. Therefore to use Eq. (4.80) to evaluate all the Ip integrals, the JP integrals must be evaluated over the range 0 S p G 2m + 4. The range of /3 for Ja in Eqs. (4.79) is 0 c p 4 2m + 2, which is a subset of the range needed for the I0 integrals. Therefore, all the Zp and JP integrals can be evaluated by evaluating Ja, 0 6 p S 2m + 4. Thus, for a dipole element (m = 2), only nine JP integrals would need to be evaluated. Perhaps the easiest way to evaluate the Jp (and Zp) integrals is to transform them to Fourier sine and cosine

14

D.S. Burnett, R.L. Holford I Comput. Methods Appl. Mech. Engrg. 164 (1998) 49-76

transforms with the change of variable x = kr - [, since Fourier routines may be found in many mathematics software packages. This is the approach used by the authors. are identical for every infinite element in a mesh because they are The radial integrals, RE!,, i = 1,2,3,4 independent of angular variables and are along identical radial paths, i.e. hyperbolas emanating from the same ellipsoid. Their computational cost is insignificant and they only need to be evaluated once for a given problem. Consequently, the numerical integration required to generate K” and M” for each infinite element involves only the evaluation of the 2-D angular integrals, making these 3-D elements as cheap to generate as 2-D elements. In addition, since the frequency dependence of the element is contained only in the radial integrals, generation of the infinite element mesh during a frequency sweep is essentially free after the first frequency. Finally, it is noted that in the limiting cases of prolate and oblate spheroidal coordinates, the above expressions for KI and ME are identical to those reported in [ 111.

5. Summary and conclusions We have presented the theoretical foundations for ellipsoidal acoustic infinite elements. A new system of ellipsoidal coordinates, comprising two angular and one radial coordinate, has been introduced. This system includes prolate and oblate spheroidal coordinates as limiting cases; hence, this element includes the authors’ previously reported spheroidal elements [ l,lO, 111 as special cases. Proof of convergence of the underlying ellipsoidal multipole expansion has been outlined and mathematical expressions for the stiffness and mass matrices have been derived and summarized with sufficient detail to enable their implementation into computer code. Numerical results for ellipsoidal elements used in ellipsoidal form are not presented here, but results in spheroidal form may be found in [ 1,111. Based on extensive theoretical and computational experience with both the boundary element method (BEM) and this infinite element method (IEM), we see the following advantages for this ellipsoidal infinite element technology. (1) When used in a prolate spheroidal form, this element has been found to be several hundred times faster than the BEM for the same accuracy, and this speed-up can increase much further depending on structural aspect ratio and frequency. This extraordinary computational efficiency is due in part (and primarily) to the local connectivity, which is characteristic of conventional FEA, and in part to the spheroidal geometry, which minimizes the number of finite-size fluid elements between the structure and the spheroid. Both features contribute significantly to decreasing the bandwidth (or wavefront) of the assembled equations. It is expected that when this element is used in ellipsoidal form, it will demonstrate even greater computational speeds because an ellipsoid can circumscribe most structural shapes more closely than a prolate spheroid, thereby reducing the size of the computational domain. (2) It is broadband. There are no critical frequencies as with the BEM, which require additional, complicated code to circumvent the problem. (3) It uses conventional FEA code architecture. Compared to the BEM, IEM code is much simpler and it is much easier to convert a structural code into a structural acoustics code. (4) This ellipsoidal infinite element employs basis functions in the outward (‘radial’) direction that are the leading terms of a new multipole expansion, which has been shown to represent the exact solution (of the exterior Helmholtz problem) in ellipsoidal coordinates. This should facilitate constructing a proof of convergence for the element, as it has already for the prolate spheroidal form [12], the latter proof being apparently the first such proof for any kind of infinite element. It should also facilitate the development of rigorous a posteriori error estimators. (5) This ellipsoidal infinite element is completely compatible with @-adaptive technology, since the FE representation is a polynomial in all directions and, as just noted, error estimators can be derived. (6) Element matrices are symmetric. This is due to the use of a conventional Bubnov-Gale&in formulation with a self-adjoint operator. (7) Element matrices are well conditioned. This is due to the separation of the 3-D integrals into products of 2-D and 1-D integrals, which, in turn, separates the oscillatory behavior over a semi-infinite interval from simple polynomials in the finite angular directions.

D.S. Burnett, R.L. Horford / Comput. Methods Appl. Mech. Engrg. 164 (1998) 49-76

75

(8) Element matrices are cheap to generate. The cost is approximately the same as for a 2-D element, again because of the separation of the 3-D integrals into 2-D and 1-D integrals. (9) This IEM technology can be extended to other forms of wave propagation in unbounded regions. Indeed, it has recently been extended to electromagnetism [32].

5. Note added in proof Since this article was accepted for publication, there has been some discussion among workers in the field of infinite elements questioning the use of the terminology ‘multipole expansion.’ This terminology was first introduced in [ 11, and has been used by us in [ 1 l] and in this paper. The consensus of these discussions is that our terminology is inappropriate, and that a preferred terminology is ‘radial expansion’ (as adopted in [lo]).

References [l] D.S. Burnett, A three-dimensional acoustic infinite element based on a prolate spheroidal multipole expansion, J. Acoust. Sot. Am. 96 (1994) 2798-2816. [2] D. Givoli, Numerical Methods for Problems in Infinite Domains (Elsevier, Studies in Applied Mech., No. 33, 1992). [3] The adjective ‘finite’ in ‘finite element’ means non-infinitesimal, not non-infinite. (See R.W. Clough, The finite element method in plane stress analysis, Proc. 2nd ASCE Conf. Electronic Computation, Pittsburgh, PA, (Sept, 1960) 345-378.) [4] P. Bettess, Infinite Elements (Penshaw, Sunderland, UK, 1992). [5] R.L. Ungless, An infinite finite element, M.A. SC. Thesis, Univ. of British Columbia, 1973. (Although the adjective ‘infinite finite’ is semantically correct, as explained in [3], all subsequent authors have dropped ‘finite’ after the word ‘infinite’.) [6] P. Bettess, Infinite elements, Int. J. Numer. Methods Engrg. 1I (1977) 53-64. [7] O.C. Zienkiewicz, K. Bando, P. Bettess, C. Emson and T.C. Chiam, Mapped infinite elements for exterior wave problems, Int. J. Numer. Methods Engrg. 21 (1985) 1229-1251. [8] R.J. Astley and W. Eversman, A note on the utility of a wave envelope approach in finite element duct transmission studies, J. Sound Vib. 76 (1981) 595-601. [9] L. Cremers and K.R. Fyfe, On the use of variable order infinite wave envelope elements for acoustic radiation and scattering, J. Acoust. Sot. Am. 97(4) (I 995) 2028-2040. [lo] R.L. Holford, A radial expansion for the acoustic field exterior to a prolate or oblate spheroid, J. Acoust. Sot. Am, submitted. [l I] D.S. Burnett and R.L. Holford, Prolate and oblate spheroidal acoustic infinite elements, Comput. Methods Appl. Mech. Engrg. 158 (1998) 117-142.. [12] L.F. Demkowicz and K. Gerdes, Convergence of the infinite element methods for the Helmholtz equation, Texas Inst. for Comp. and App. Math., Report 95-07, Univ. of Texas at Austin, May 1995. [13] K. Gerdes and L.F. Demkowicz, Solution of 3D-Laplace and Helmholtz equations in exterior domains using hp-infinite elements, Comput. Methods Appl. Mech. Engrg. 137 (1996)239-273. [14] L.F. Demkowicz and F. Ihlenburg, Convergence analysis for coupled finite/infinite element methods, Proc. Fourth U.S. Nat. Congress Comp. Mech., San Francisco, Aug. 1997. [15] K. Gerdes, The conjugated vs. the unconjugated infinite element method for the Helmholtz equation in exterior domains, Eidgenossische Technische Hochschule, Zurich, Switz., Research Report 96-l 1, July 1996. [I61 A.A. Oberai, P.M. Pinsky and M. Malhotra, A study on infinite elements and NRBCs for the exterior acoustics problem. Proc. Fourth U.S. Nat. Congress Comp. Mech., San Francisco, Aug. 1997. [I71 J.J. Shirron and I. Babuska, Infinite element methods for exterior Hehnholtz problems, Proc. Fourth U.S. Nat. Congress Comp. Mech., San Francisco, Aug. 1997. [ 181 R.J. Astley and J.P. Coyette, Mapped spheroidal wave-envelope elements for unbounded wave problems, Proc. Fourth U.S. Nat. Congress Comp. Mech., San Francisco, Aug. 1997. [I91 M. Malhotra, A.A. Oberai and PM. Pinsky, On the use of infinite elements and the DtN condition for solving large-scale problems in acoustics, Proc. Fourth U.S. Nat. Congress Comp. Mech., San Francisco, Aug. 1997. [20] J.L. Cipolla, Transient infinite elements for acoustics and shock, DE-Vol. 84-2, 1995 Design Eng. Tech. Conf.,Vol. 3 - Part B, ASME, 113-128 (1995). [21] R.J. Astley, A transient infinite element for multi-dimensional acoustic radiation, DE-Vol. 84-2, 1995 Design Eng. Tech. Conf., Vol. 3-Part B, ASME, 97-111 (1995). [22] P.M. Morse and H. Feshbach, Methods of Theoretical Physics, Chapter 5. (McGraw-Hill, New York, 1953). [23] Although it is not germane to this discussion, we may mention that, among all ellipsoids that surround a convex structure, K, there exists exactly one that has minimum volume. This ellipsoid is called the minimum ellipsoid or Luewner ellipsoid c~j’K. [See H. Buseman. The Geometry of Geodesics (Academic Press, New York, 1955). Sn. 16.1 [24] E.W. Hobson, The Theory of Spherical and Ellipsoidal Harmonics, Chapter 11.(Cambridge University Press, Cambridge, UK. 1931). [25] S. Wolfram, Mathematics: A System for Doing Mathematics by Computer, 2nd edition (Addison-Wesley, NY, 1991).

16

D.S. Burnett, R.L. Holford I Comput. Methods Appl. Mech. Engrg. 164 (1998) 49-76

[26] M. Marden, Geometry of Polynomials, Chapter 10. (American Mathematical Society, Providence, RI, 1966). [27] E.I. Jury, Theory and Application of the z-Transform Method, Chapter 3. (Wiley, NY, 1964). [28] L.F. Kallivokas, J. Bielak and R.C. MacCamy, A simple impedance-infinite element for the finite element solution of the three-dimensional wave equation in unbounded domains, Comput. Methods Appl. Mech. Engrg. 147 (1977) 235-262. [29] B.A. Szabo and I. Babuska, Finite Element Analysis (Wiley, New York, 1991). [30] D.S. Burnett, Finite Element Analysis: From Concepts to Applications (Addison-Wesley, Reading, Mass., 1987). [31] Ibid., p. 578. [32] D.S. Burnett and F.M. Labianca, A 3-D electromagnetic infinite element for wave propagation in unbounded domains, in preparation.