A comparison of approximate boundary conditions and infinite element methods for exterior Helmholtz problems

Computer methods in applied mechanics and engineerlng ELSEYIER

Comput.

Methods

Appl. Mech. Engrg.

164 (1998)

121- 139

A comparison of approximate boundary conditions and infinite element methods for exterior Helmholtz problems Joseph J. Shirron”‘“, Ivo BabGkab “Naval Research Laboratory, Code 7131, Washington, DC 203755320, USA bTICAM, The lJniversi@ of Texas at Austin, Austin, TX 78712, USA Received

3 November

1997

Abstract We compare the accuracy of a variety of approximate boundary conditions (ABCs) and two versions of an infinite element method (IEM) for approximating solutions of exterior Helmholtz problems. The infinite elements are of variable order in the radial direction and differ only in whether the test functions in the weak formulation are conjugated or not. We find that the accuracy of the ABCs decreases with increasing frequency, and that at higher frequencies their accuracy is several orders of magnitude lower than that obtained with the IEMs. We show that the convergence rate of the conjugated element is slower than the rate for the unconjugated element and that the latter element fails to converge in the far field. 01998 Elsevier Science S.A. All tights reserved.

1. Introduction The need to solve Helmholtz problems in unbounded domains arises in many areas of physics and engineering, for example in the analysis of acoustic scattering or radiation from a submerged body. When the dimensions of the body are small in comparison to the extent of the surrounding medium it is common to take the medium as unbounded. This assumption imposes difficulties on the development of numerical methods for solving such problems. One common method is to reformulate the exterior problem as a boundary integral equation over the surface of the body [28]. This reformulation has the advantage of automatically incorporating the Sommerfeld radiation condition [9,58] at infinity. Furthermore, when an approximate solution is computed using a boundary element method, only the surface of the body need be discretized, as opposed to the exterior volume. Hence, mesh construction is greatly simplified. In addition, this method involves only an unknown surface field, rather than a field defined on a volume, so the size of the linear system of equations for the discretized boundary integral equation is relatively small due to the reduction of the dimension of the problem. On the other hand, the coefficient matrix is dense and (usually) non-symmetric. This method also requires numerical techniques for the evaluation of oscillatory and (at least) weakly-singular integrals. Then, there is the difficulty in posing boundary integral equations that are uniquely solvable for all excitation frequencies. Several schemes to overcome this difficulty have been proposed, for example [22,26,27,47,54,56,60], but usually at the cost of increased algorithmic complexity and computational expense. The reader may also consult [2] for additional details and references on boundary element methods for the Helmholtz equation. Studies concerning the cost effectiveness of boundary element methods versus finite element methods (IEMs)

* Corresponding

author

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

122

J.J. Shirron,

I. Babufka

I Comput.

Methods

Appl.

Mech.

Engrg.

164 (1998)

121-139

for exterior Helmholtz problems have been conducted by Harari and Hughes [45]. They conclude, as Bettess [ 171 does for simpler problems, that finite element methods are more cost-effective for the test cases they report. Burnett [23] has also reached a similar conclusion. Our experience is in agreement with these researchers, so we shall limit our discussion from here on to FEMs for exterior problems. One natural way to use a FEM for the solution of a problem in an unbounded domain is simply to truncate the infinite exterior. A new problem on the resulting finite domain is then posed whose solution, ideally, should coincide with the restriction of the solution of the original problem. The main difficulty is in formulating a boundary condition on the truncation surface to replace the Sommerfeld radiation condition at infinity. For simple truncation surfaces, such as a circle in two dimensions or a sphere in three dimensions, the exact Dirichlet-to-Neumann (DtN) boundary condition is known. Keller and Givoli [42,50] have used the exact boundary condition with some success. However, the exact condition is given in an infinite series which must necessarily be truncated in any actual computation, and doing so introduces error. In addition, the exact boundary condition is nonlocal which, upon discretization by a FEM, results in the filling of dense blocks in an otherwise sparse system of linear equations. For this reason we restrict our attention to local approximate boundary conditions (ABCs) that preserve the sparsity and symmetry of the problem. In particular, we consider the conditions suggested by Bayliss and Turkel [12], Engquist and Majda [37,38] and Feng [48]. Research on the DtN is ongoing, and advances have been made in constructing local approximations to the nonlocal DtN map and in developing special high-order finite elements for use with high-order DtN conditions [32,41,55]. Other papers of interest include [44] and [4]. Exterior problems may also be solved by a combination of finite and infinite elements. In this case, instead of truncating the infinite domain one uses finite elements to discretize the domain out to a circle in Iw2 or a sphere (or spheroid) in [w3. The domain outside the circle or spheroid is discretized using infinite elements, which are elements having infinite extent in the radial direction. The infinite element shape functions are typically, though not always, patterned after the asymptotic expansion of the solution of the Helmholtz equation. The salient feature is their satisfaction of the Sommerfeld radiation condition. There is a fairly large body of literature devoted to the study of infinite elements, for example see [3,5-8,14-16,18-21,29-31,361. In this paper we shall report results for the infinite element formulation of [57] which is based directly on Burnett [24]. Our goal in this paper is to quantitatively assess the error in approximating solutions of exterior Helmholtz problems using a FEM together with ABCs or IEMs. We make this assessment for two simple model problems for which the exact solution is known: acoustic plane wave scattering from a rigid circle and a rigid sphere. In Section 2 we discuss the geometry for a problem in 2D and give its exact solution. We also develop a weak formulation of the problem on a truncated domain and give a condition for its unique solvability. Various approximate boundary conditions in 2D are discussed in Section 3, followed by two infinite element methods in Section 4. Section 5 contains a discussion of the results obtained for scattering by a circle. Details are given there concerning the construction of the meshes and our choice of error criterion. We tabulate results for the various ABCs and IEMs, including the size of the linear systems, their condition number and time for solution. At the end of the section we consider the quality of the approximation in the far field, as computed from the infinite elements. In Section 6 we introduce the model problem in 3D and give its exact solution. We consider only the ABCs of Bayliss and Turkel and show how they differ from the ones in 2D. Also, we discuss the modifications to the IEMs necessary in 3D. Section 7 follows with the results for the sphere, along with a brief discussion of the similarities to the results for the circle. A generalization of the infinite element methods to prolate spheroidal coordinates is made in Section 8 and a comparison of the two IEMs is made there for the case of a radiating hemispherically-capped cylinder. The paper ends with a summary and various comments in Section 9.

2. A model problem in Rz Here, we consider the problem of approximating wave. (One may also regard this as scattering cross-section.) This is a Neumann problem, but we r = {XE [w*: 1x1= u} be the circle of radius a

the wave scattered from a rigid circle insonified by a plane of a plane wave by an infinite cylinder with a circular obtain similar results for Dirichlet and Robin problems. Let and let R, = {x E I!%*: In\> a} denote its exterior. Polar

J.J. Shirron, I. Babushka / Comput. Methods Appl. Mech. Engrg. 164 (1998) 121-139

I??

coordinates (r, (3) are used throughout this section. A plane wave uinc(r, 0) = e”rcoaH with wave number k > 0 is incident on r from 0 = 0. The scattered wave K satisfies the Helmholtz equation Au+k2U=0

inn,,

the rigid (or sound-hard) au -=-an

condition

on r,

an

$

U

boundary

aUlnC

and the Sommerfeld

uniformly

(1)

(2)

radiation

condition

(r, 13)- iku(r, 0)

in 19.The incident

i”c(r, 6) =

wave has the expansion

2 i”J,(kr) m=--z

where _I, is the cylindrical series

= 0, [_51]

eime ,

Bessel function

of the first kind of order m. The exact solution is then given by the

r

uen(r,~9) = -

C

m=-z

i” j$$$ m

H,(kr) eime ,

where H,,, is the cylindrical Hankel function of the first kind of order m and the prime denotes differentiation with respect to argument. To approximate the solution of the scattering problem using ABCs the infinite exterior is truncated at a radius of R > a. Let r, = {XE R2 : 1x1= R} denote the truncation circle and let OR = {x E R2 : a < 1x1CR} denote the finite computational domain. See Fig. 1 for a depiction of the geometry. The solution of the original scattering problem is approximated by the solution u of the problem Au + k’u = 0

in OR ,

(6)

au -=-dn

on r,

(7)

aU’nC an

and

Fig. 1. Geometry

for the problem

in the truncated

domain.

J.J. Shirron, I. Babushka I Comput. Methods Appl. Mech. Engrg. 164 (1998) 121-139

124

au

j-yBu=O

on&.

The Robin boundary condition (8) replaces the Sommerfeld problem is to find a function u E H'(L&)such that

I

(Du . Du - k’uu) d0k -

I

vBu dT, =

r,

%

I

for all v E H’(f&). For the problem on the B : L2(rR) +L2(TR) must satisfy the inequality

1

radiation

condition.

A weak formulation

of this

auinc

vxdr

truncated

domain

to be

uniquely

solvable

the

operator

(10) for all u E L2(rR) and some constant C > 0, where the overbar denotes complex conjugation. As we mentioned earlier, the exact operator B (also called the Dirichlet-to-Neumann map) is nonlocal while the ABCs considered in this paper are all local. Consequently, the solution u will not agree with the exact solution ueX. Our intention is to quantify this error for various ABCs.

3. Some ABCs in R? The simplest ABC is obtained the finite radius R, i.e.

au

by imposing

the Sommerfeld

radiation

condition

not at infinity but, instead, at

on r, .

dn - iku = 0

(11)

Among the more widely used ABCs are those of Bayliss and Turkel [ 121 (see also [ 111). They define a hierarchy of operators

(12) that successively annihilate also Karp [49]), ikr

the first n terms in the far-field expansion

of Atkinson

[9] and Wilcox

[61,62] (see

zc F,#>

r”

(13)

’

that is, L&R, 0) = 0(R-n-“2). Their boundar y conditions follow quite naturally first-order Bayliss and Turkel boundary condition is simply $-(ik-&)u=O

on 4,

and their second-order

condition

by setting L,u(R, 0) = 0. The

(14)

is

(15)

(it-$)$++(2k’+T---&+-$-$)u=O.

Another sequence of local ABCs of increasing order was developed by Engquist and Majda [37,38]. Their method is based on the theory of pseudo-differential operators and its explanation is beyond the scope of this paper. We simply note that their first-order condition is identical to the first-order Bayliss and Turkel condition while their second-order condition is

at.4

an-

C

ik--+

1

2R

L(ik+i)$]u=O 2k2R2

on&.

(16)

J.J. Shirron, I. Babu&

I Comput. Methods Appl. Mech. Engrg.

164 11998)

121- 139

125

Still, another sequence of local ABCs was proposed by Feng [48]. Feng derived his conditions by truncating an asymptotic expansion of the exact Dirichlet-to-Neumann (DtN) map. His first-order condition coincides with the first-order condition of Bayliss and Turkel, and his second-order condition is

au --

(17)

an Feng’s third-order

au

-dn

ik--+

condition 1 2R

is

--!--[ik+i)(l+4-$)]u=O 8k2R2

on4.

(18)

4. Infinite element methods in Rz In the preceding section the original exterior domain was truncated and a new problem was posed on the finite domain, with a Robin condition (ABC) on the truncation boundary replacing the Sommerfeld radiation condition at infinity. In this section we consider infinite element methods in which the exterior is discretized with elements of infinite extent in the radial direction whose shape functions satisfy the Sommerfeld radiation condition. The surface r,, which denoted the truncation boundary in the previous section, will now denote the boundary between the finite element domain and the infinite element domain. As with the FEM/ABC methods, the solution u on the finite domain 0, is sought in the Hilbert space H’(,&). But now the solution in the infinite domain L?; is sought in terms of an ansatz suggested by the Atkinson/Wilcox expansion, i.e.

with the function through

p EH’(D),

Cu(r, 6)

where D = (0, 1) X (0, HIT).We define a function

U over the exterior

aSrcR,

where continuity at r = R is satisfied if u(R, 0) = ,u( 1,8) for 0 < 8 s 2~. A weak formulation of the problem is to find two functions u E H’(LnR) and ,LLE H’(D) such that

for all functions

domain

with u(R, .) = p( 1, 0)

V of the form

with u E H ‘(,flR), v E H’(D) and u(R, -) = v( 1, e). Here, 0: is the annular domain bounded by the circles of radius R and S. Note that the limits of the integrals in the first term of (21) cannot be evaluated independently. Separately, the limits of the integrals may be unbounded, but the limit of the combination of the area and line integrals is finite [24,57]. Another weak formulation of this problem can be obtained by simply replacing V in (21) by its complex conjugate. In this case it turns out that the limits of all the integrals are bounded and can be evaluated independently [57]. The weak formulations given above are by no means the only ones possible or the only ones of practical value. Variational formulations of the problem have been developed in weighted Sobolev spaces by Demkowicz and Gerdes [34,40], who have carried out a detailed analysis of the convergence and stability properties of such formulations as well as the formulation of Burnett [24]. The reader should also consult the report by Demkowicz

J.J. Shirron, I. Babushka I Comput. Methods Appl. Mech. Engrg. 164 (1998) 121-139

126

and Ihlenburg [35] which analyzes the convergence of the coupled finite/infinite element methods for the Helmholtz problem. The details of the infinite element shape functions are the same for either formulation. Denote the angular variable by 5 E [ - 1, l] and the reciprocal radius by the variable t = R/u (so in the domain 0; we have 0 < t S 1). Then, we construct the shape functions as tensor products of ‘radial’ and angular functions, (23) For the ‘radial’ functions

we use

T,(t) = 1 ,

(24)

T,(t)=PT(t)-

1

for j>O,

(25)

where PT is the shifted Legendre

polynomial

of degree j. For the angular functions

we use

N,(5) =;

(1 - 0,

(26)

N,(5)=&1

+I$),

(27)

Nk(t)=pk(6)-pk-2(t)

fork>1,

(28)

where Pk is the Legendre polynomial of degree k. Within an infinite element the angle 6’ depends on 5 in a way that is independent of the number of functions Nk retained in ,u. Rather, we implement a pq-version of the element in which the geometry is fixed and we may vary the orders p and q. (In the results that follow, the geometry is approximated by quadratic Lagrange interpolation of the nodal coordinates at element vertices and mid-sides.) In the finite domain flR we use the p-version of the finite element method [59]. It is then trivial to build the continuity condition u(R, .) = p( 1, *) into the approximation spaces. The field on an element in fi, with an edge on r, is represented by a linear combination of Nk for k = 0, . . . , p. The same is true for the adjacent infinite element at t = 1. So the only difficulty is in matching the coefficients hk to the coefficients of the finite-element edge.

5. Comparison of results in R2 Recall in the model problem a plane wave is scattered by a rigid circle of radius a. Here, we report results in which, for the ABCs, the exterior domain is truncated at a circle of radius R = a + A, where A = 2n/k is the acoustic wavelength. In other words, the computational domain is the annulus whose thickness is one wavelength. For the IEMs, the infinite elements model the domain 0; in which r 5 R. We present results for the non-dimensional frequencies ku = 1,5, 10. In each case the mesh for 0k is uniform with 16 elements in the radial direction and 16ka elements in the angular direction. This provides a resolution of 16 elements per wavelength in both the angular and radial directions. The finite elements are all cubic (p = 3) quadrilaterals. The infinite elements have cubic angular variation and are polynomials of degree q in the reciprocal radius t. For of r-l’* and re3’*, etc. q = 0 the radial decay is just r -“*, for q = 1 the decay is a linear superposition The discretization is sufficient so that nearly all of the error in the computed solutions can be attributed to inaccuracy of the ABC or IEM. Nevertheless, the mesh still limits the radial and angular resolution to which a solution can be computed. We shall use the phrase ‘threshold limit of the mesh’ to refer to the lowest possible error that can be achieved in a computation on that mesh. Our measure of error is the L, norm of the restriction of u - uex to r (i.e. the trace of u - uex on r),

((u- uexl(2 = =u

lr(u I

uex12 dT

ozn Ju(u, 0) - ~““(a, @)I*d0 .

(29)

J.J. Shirron, I. BabuSko I Comput. Methods Appl. Mech. Engrg. 164 (1998) 121-139

We deem the trace of u on r to be the most important be computed via the Helmholtz integral 1 u(x) = 5

$

(Y)Y(~, Y) - U(Y) s

?

computed

quantity.

The solution at any point in 0,

127

can

(31)

(n, y)] dT(y), x E n+ ,

where the free-space Green function fix, y) = -i/2H#lx - yl) and HA is the cylindrical Hankel function of the first kind or order zero. Obviously, the accuracy in computing u at any point x E 0, is limited only by the error in u on r since the normal derivative is known through (2). For this reason we are not concerned with assessing the error over the whole domain 0, or with assessing the error in the derivatives of u. (Note that this would not be the case if we were solving a Dirichlet problem. Then we would be interested in the error in the normal derivative of u on r). We also report the size of the system of linear equations that was solved in each case. The size is given as N X b, where N is the dimension of the system and b is the half-bandwidth. We have used the LAPACK [39] routines zgbtrf for the factorization and zgbtrs for the forward and backward substitutions. The condition number of each system is estimated with the routine zgbcon and is tabulated along with the time (in seconds) that it took a midrange workstation to compute the solution of the banded system of equations. We note that the systems of equations for the ABCs and the unconjugated infinite elements (UIEs) are complex symmetric (not Hermitian), while the systems for the conjugated infinite elements (CIEs) are non-symmetric. If LAPACK routines existed for the solution of complex symmetric banded systems, the solution times for the ABCs and CIEs could be roughly halved. Such routines do not exist, only ones for Hermetian positive definite systems are available. Unfortunately, the Helmholtz operator is not positive definite. However, it is not very difficult to develop a custom solver for complex symmetric, banded systems of equations (we have done so) and thereby exploit the symmetry to reduce the solution time and memory requirements by half. We do not report results obtained with our routines, for uniformity we use the LAPACK routines for general banded systems in each case. Table 1 shows the results obtained for the non-dimensional frequency ku = 1. Among the ABCs, the Sommerfeld (Somm) condition leads to the greatest error in trace, as expected. The second-order condition of Bayliss and Turkel (B&T-2) produces the least error. The third-order ABC of Feng (Feng-3) is also very effective. The results for the unconjugated infinite elements with q = 0, 1 (UIE-0 and UIE-1, respectively) are nearly identical. Results from Shirron [57] indicate exponential convergence of this element, so we conclude that the error obtained for UIE-1 represents the threshold for this mesh. Note that the condition number for UIE-1 is two orders of magnitude larger than for UIE-0. The trace error for the conjugated infinite elements also appears to be decreasing exponentially with increasing order q, with the threshold error reached by CIE-3. Note that, in contrast to the unconjugated case, the condition numbers for the conjugated elements are relatively insensitive to the order q. At this low frequency the system sizes are small and the equations are solved very quickly.

Table

I

Plane wave scattered by a circle of radius a at ka =

I

Method

llu - UCXII

System size

Cond. no.

CPU

Somm

8.695e - 2

1328 X 155

1.183e +4

3.3

B&T-l

6.846e - 3

1328 X 155

9.276e + 3

3.3

B&T-2

5.654e - 4

1328 X 155

7.337e + 3

3.3

E&M-2

3.036e - 3

1328 X 155

7.294e + 3

3.3

Feng-2

1.096e - 3

1328 X 155

8.104e + 3

3.3

Feng-3

5.863e - 4

1328 X 155

7.294e + 3

3.3

UIE-0

5.675e-4

1328 X 155

7.380e + 3

3.3

I

5.647e - 4

1376 X 151

9.31 le + 5

3.3

CIE-0 CIE- 1

9.916e - 2

1328 X 155

3.080e + 4

3.3

I .730e - 3

1376 X 151

2.258e + 4

3.3

CIE-2

5.663e - 4

1424 X 166

2.837e + 4

3.9

CIE-3

5.647e - 4

1472 X 172

2.757e f 4

4.9

UIE-

128

J.J. Shirron, 1. Babushka I Comput. Methods Appl. Mech. Engrg. 164 (1998) 121-139

Table 2 Plane wave scattered by a circle of radius a at ka = 5 Method

llu - UeXll

System size

Cond. no.

Somm B&T-l B&T-2 E&M-2 Feng-2 Feng-3

1.301e 9.284e 4.816e 3.905e 1.065e 4.178e -

6640 6640 6640 6640 6640 6640

3.571e + 4.17Oe + 3.337e + 3.341e + 3.554e + 3.340e +

WE-0 UIE- 1 UIE-2

5.146e - 3 6.664e - 5 6.439e - 5

CIE-0 CIE- 1 CIE-2 CIE-3 CIE-4 CIE-5

2.867e + 3.147e 3.662e 2.441e 1.188e 6.441e -

1 2 3 3 2 3

0 1 2 3 4 5

X X X X X X

324 324 324 324 324 324

CPU 3 3 3 3 3 3

72 12 80 80 81 82

6640 X 324 6880 x 340 7120 X 378

3.336e + 3 4.948e + 5 1.148e + 8

80 17 90

6640 6880 7120 7360 7600 7840

1.235e + 5 l.O90e+5 6.575e + 4 9.302e + 4 1.254e + 5 1.003et5

71 81 92 134 192 264

x 324 x 340 X 378 X 484 X 590 X 696

Table 2 shows the results obtained for ku = 5. Again, the Sommerfeld condition yields the largest error. The second-order ABCs of Bayliss and Turkel and Engquist and Majda produce errors on par with Feng-3. Exponential convergence with q is evident in the results for the unconjugated infinite elements. The trace error for UIE-2 is only slightly lower than that for UIE-1, leading us to conclude that the error for UIE-2 is the least error that can be obtained for this mesh. This conclusion is supported by the results for the conjugated elements. They also converge exponentially with order q (albeit more slowly) and they reach the same threshold. Just like the results for ka = 1, the condition numbers for the unconjugated elements are increasing exponentially with q while the conjugated elements show no growth trend. It is important to note that for the infinite elements the system sizes and bandwidths are increasing with q. Accordingly, the solution time increases, and the effect is not trivial. The solution time for CIE-5 is nearly four times that of CIE-0. (We note here that the variation in times for solving systems of the same size is due to differing run time loads on the computer.) Table 3 displays the results obtained for the high frequency ku = 10. Once more the Sommerfeld condition is the least accurate. However, now the difference between the other ABCs is considerably diminished. The trace errors for the first-, second- and third-order conditions are nearly identical, and they are several orders of Table 3 Plane wave scattered by a circle of radius a at ka = 10 Method

IIU- UeXll

System size

Cond. no.

CPU

Somm B&T-l B&T-2 E&M-2 Feng-2 Feng-3

2.204e 2.061e 2.190e 2.07Oe 2.793e 2.12Oe -

13280 X 13280 X 13280 X 13280 X 13280 X 13280 X

2.695e + 2.69Oe + 3.270e + 3.273e + 3.29le + 3.272e f

3 3 3 3 3 3

158 153 166 167 166 166

UIE-0 UIE- 1 UIE-2

2.224e - 2 3.058e - 4 7.03Oe - 5

13280 X 324 13760 x 340 14240 X 378

3.269e + 3 5.724e + 4 3.358e + 8

168 163 195

CIE-0 CIE-1 CIE-2 CIE-3 CIE-4 CIE-5 CIE-6

2.239e + 1.458e + 6.921e 1.53Oc 2.704e 3.524e 3.760e -

13280 X 13760 X 14240 X 14720 X 15200 X 15680 X 16160 x

5.485e + 5 3.306e + 5 2.035e + 5 1.433e + 5 1.923e + 5 1.92Oef5 2.134e + 5

154 176 215 322 423 576 803

1 2 2 2 2 2

0 0 1 1 2 3 4

324 324 324 324 324 324

324 340 378 484 590 696 802

J.J. Shirron, I. Babushka I Comput. Methods Appl. Mech. Engrg.

164 (1998) 121-139

129

magnitude higher than the errors obtained by the most accurate infinite elements. For the unconjugated elements, the trace errors decrease exponentially while the condition numbers increase at a similar rate. The errors also drop exponentially with order q for the conjugated elements, but the condition numbers are stable. Lack of core memory prevents us from running the case CIE-7, but the trend indicates that the trace error in this case would be limited by the error threshold for the mesh. Once more we call attention to the fact that the solution times for the infinite elements are increasing with order q. The solution time for CIE-6 is roughly five times that for CIE-0. An intriguing possibility presents itself when an exterior Helmholtz problem is solved with infinite elements, namely, whether the solution can be obtained accurately throughout the infinite domain. When the problem is solved using an ABC the field is known only over a finite computational domain, and one must resort to the Helmholtz integral to compute the field outside that domain. However, when the problem is solved with infinite elements the entire exterior is modeled and it is natural to ask about the quality of the solution in the infinite element domain. The answer is perhaps nonintuitive: an accurate solution can be obtained with conjugated elements (by taking q large enough) but an accurate solution is not always available with unconjugated elements (regardless of the value of q). We illustrate this for our model problem. But first let us introduce the function pex defined through (32) which represents the variation of the exact solution after factoring out the asymptotic oscillation and decay. This is precisely the ansatz assumed in (19), so we can directly compare the exact function pex with that computed using either version of the infinite elements. Figs. 2-4 are plots of the real part of ,u as a function of R/r at the angle 0 = 0 for the three non-dimensional frequencies ka = 1,5, 10, respectively. In each case the field in the conjugated element very accurately approximates the exact field. The field computed from the unconjugated element agrees well at r = R then diverges from the exact field as r+m. In Figs. 2 and 3 the conjugated and unconjugated elements are compared for the same order q, which may be misleading since the elements converge at different rates. In Fig. 4 the results are plotted for the conjugated elements of order 2 and 4; the curve for CIE4 is indistinguishable from the exact curve. The angle 0 = 0 corresponds to the direction of backscatter of the incident wave, but is otherwise not special. Similar results are found at other angles and also for the imaginary part of p. We conclude from these results that the unconjugated infinite elements serve to provide a highly accurate boundary condition on r,, but fail to provide an accurate approximation of the solution in 0; (that is, in the infinite element domain itself). With the conjugated infinite elements (and large enough q) a highly accurate approximation of the solution can be computed in all of a+ (that is, in both the finite element and the infinite element domains). On the other hand, the trace error on r for the unconjugated infinite elements converges much faster than for the conjugated elements. The inability of the UIE to provide accurate solutions out to the far field is not related to the approximability

-0.17,

,

Scaled ExteriorField: k&l I

,

.

-

I Exact

o-----oCl&1

-0201 0.0

0.2

’

0.4

0.6

h

0.8

FUr

Fig. 2. Real

part of scaled exterior solution at 0 = 0, ka = 1

J.J. Shirron, I. Babushka I Comput. Methods Appl. Mech. Engrg. 164 (1998) 121-139

130

sceled Exteriw Field:ka=5

,

0.80,

,

I

4.5,

,

Scaled ExteriorField:ka-10

1

,

(

P

-

‘9 ‘8

3.5 -

Exact

o-..-----o CIE_2

‘9

* - - 0 WE-2 0.80 -

b--d

i_/

,

0.0

0.2

0.8

0.4

,

T

0.8

1.o

I

0.5’

0.0

0.2

3

I

0.4

FUr

0.8

CIE-4

1

0.8

R/r

Fig. 3. Real part of scaled exterior solution at 0 = 0, ka = 5. Fig. 4. Real part of scaled exterior

solution at 0 = 0, ka = 10.

properties of the trial spaces. Recall that both IEMs are based on the same variational form of the exterior problem, the difference is only in the test spaces. At the infinite-dimensional level both IEMs are equivalent. In particular, the trial spaces are the same in each method, so the approximability properties are the same for the UIE and CIE methods. Our infinite element trial spaces are tensor products of eik’/r times polynomials in t = R/r by piecewise polynomials in the angular direction. Results from [57] demonstrate that the spaces V, = span{t”}z are well suited to approximating the radial variation of the exact (scaled) solution, with the error of the best approximation decreasing exponentially in L,(O, 1) with increasing q. The suitability of the trials spaces to approximating the exact solution is demonstrated by the convergence of the CIE solutions for all values of r 3 R. The failure of the unconjugated elements to accurately approximate the solution as r +a is a direct consequence of the instability of the bilinear form in (21). More precisely, the bilinear form A(U, V) = A( p, v) =f,”

(DU . DV- k2UV) doi

cc

CQ1 UV dT,

- ik J

(33)

with

(35) and ,u, v E H = H’(D) yN =

min

is unstable.

By this we mean that the inf-sup

constants

]A(K V)l

max

O#fiEH,,,O#vEHN

Ill4 II4

(36) '

for the bilinear form A(*, a) on the finite-dimensional subspaces HN C H are not bounded from below as N + ~0. (Here, )(*I]is induced by some inner product on the Hilbert space H, and N depends on the angular discretization, the element angular order p and ‘radial’ order q.) Let us examine the consequences of this. Let uN solve the problem Find uN E H, such that A@,, v,) = F(u,), where F is a bounded

linear functional

V uN E HN ,

(37)

on H. From BabuSka [lo] we have that if (38)

J.J. Shirron, I. Babushka I Comput. Methods Appl. Mech. Engrg. 164 (1998) 121-139

then the error is bounded

131

by

llUeX-UN(14(1+cCNIYN)IJUeX-WNII,

(39)

Vw,EH,,

where C, is a bound for A(., *) on H,,,, i.e. A+,,

u,) c C,.&,,$l/u,,l] ,

(40)

v u,,r, u, E H,v .

In [57] we show (for the L,(D) norm) that there is a sequence of functions up E HN (the best approximation of with element order q for sufficiently fine angular ueX) such that 11~~~- u:]] converges to zero exponentially discretization. On the other hand, we also show that the stability constant C,/y, grows exponentially with q, is increasing with q. Consequently, the UIE and that this growth is so rapid that the quantity C,,,llueX- uzl]/r, method is simply not convergent in the L, norm over (0, 1) X (0,2n). Nevertheless, we are able to prove convergence of the unconjugated elements on r, (i.e. at r = R) by a line of reasoning that depends on the symmetry of A(., .) see [57] for more detail. The instability of the bilinear form in the UIE method manifests itself in the ill-conditioning of the linear systems computed after discretization. This ill-conditioning has consequences for iterative solution techniques of the linear equations. Both theory and experience show that the convergence rate is either slowed or obviated for ill-conditioned systems. In view of the exponentially increasing condition numbers that we have calculated for the FEM/UIE systems, it may be inappropriate to attempt iterative solutions for larger values of q. See [52,53] for some experience by Malhotra and Pinsky with iterative solutions and their methods for dealing with the ill conditioning of the infinite element systems.

6. A model problem

in lR3

We now consider the problem of computing the field scattered by a rigid sphere insonified by a plane wave. For the remainder of this paper we redefine the surface r = {x E Iw3 : 1x1 = a} and the exterior domain 0, = {XE [w3: 1x1~ a}. In spherical coordinates, a plane wave uinc(r, 0) = elkr “I ’ with wave number k I=-0 is incident on r from 8 = 0. The scattered wave u satisfies the Helmholtz equation (1) and the Neumann boundary condition (2), as well as the Sommerfeld radiation condition !iEr uniformly

1

g(6)-iku(re^)

I

=O,

for all Jil = 1. The incident

(41) wave has the expansion

[l],

m

I.4

i”c(r, 0) = 2 (2n + l)i”P,(cos 11=il

0)j,(kr),

(42)

independent of the azimuthal angle 4. The function j, is the spherical Bessel function of the first kind of order IZ. The exact solution is given by the series r j

uex(r, 0) = - C (2n + l)i”P,(cos n=o where h, is the spherical

Hankel function

R Uex(r, 0) = elkcrpR) 7 pex

’W

0)“-h;(k,.) h,(kr) ’

,

(43)

of the first kind of order n. For later use we also define ,uex through forrZ=R,

which represents the variation of the exact solution after factoring out the asymptotic oscillation and decay behavior. For the ABCs the infinite exterior is truncated at a radius of R > a. We redefine 4 = {XE [w3: 1x1==R} to denote the truncation sphere and L$ = {XE [w3: a < 1x1CR} to denote the finite computational domain. For the IEMs r, is the interface between the finite element domain L& and the unbounded exterior L?; = {XE [w3 : /xl> R).

132

J.J. Shirron, I. Babushka I Comput. Methods Appl. Mech. Engrg.

164 (1998) 121-139

In three dimensions we consider only the ABCs of Sommerfeld, which is unchanged and Turkel. The Atkinson/ Wilcox expansion in spherical coordinates is

from (1 l), and Bayliss

(45) The family of Bayliss and Turkel boundary conditions, which again annihilate expansion (45), follow from setting L,,u = 0 on r,, where L”=,fj The first-order

(-$-ik+v),

successive

terms in the far-field

nal.

(46)

ABC is simply

&A

-_ an

on r,,

and the second-order

condition

(47)

is

2(ik-:)$+(T--$+2k’)u+&$(sinf3%)+&-=0

on&.

(48)

The weak formulation of the problem on the truncated domain is identical to (9) and is subject to the same solvability condition (10). The development of the infinite elements in three dimensions parallels the treatment in 2D, except the ansatz is now e with the function

,u E H’((0, 1) X (0, T) X (0,2n)).

We define the function

U over the exterior domain through

where continuity at r = R is satisfied if u(R, 8, 4) = p( 1, f3,4) for 0 =Z0 G n and 0 G 4 G 21r. The weak formulation of the problem is analogous to(21) with obvious changes to account for the extra dimension in 4. As before, the limits of the integrals cannot be evaluated independently since the surface integral cancels unbounded terms in the volume integral [24,57]. Just as in 2D, another weak formulation can be obtained by merely replacing the test functions with their complex conjugates. And in this case it also turns out that the limits of all integrals are bounded and can be evaluated independently.

7. Comparison of results in [w3 Our model problem in 3D is axisymmetric, i.e. the solution is independent of the azimuthal angle 4. For simplicity, we shall exploit this fact and explicitly remove the 4 dependence throughout. This means that the surface r reduces to the semi-circle of radius a, that rR reduces to the semi-circle of radius R, that the volume OR reduces to the area {x = (r, 0) : r S a =SR, 0 S 13S T} etc. Under this assumption the problem is essentially two-dimensional. As in the circular case, we consider the three non-dimensional frequencies ka = 1,5, 10. We construct meshes for 0, that are uniform with 16 elements in the radial direction and 8ka elements in the angular direction. Again, this provides a resolution of 16 elements per wavelength in both the polar angular and radial directions. The finite elements are all cubic (p = 3) quadrilaterals. The infinite elements have cubic variation in 19and are polynomials of degree q in the reciprocal radius t. Hence, for 9 = 0 the radial decay is just r-‘, for 9 = 1 the decay is in terms of r-’ and rP2, etc.

J.J. Shirron, I. Babushka ! Comput. Methods Appl. Mech. Engrg. 164 (1998) 121-139 Table 4 Plane wave scattered

by a sphere of radius a at ka =

133

I

Method

I/u - UrXll

System size

Cond. no.

CPU

Somm B&T- 1 B&T-2

3.521e - 2 2.27le - 3 4.255e - 4

713 x 86 713 x 86 713 x 86

1.184e+5 1.145e + 5 1.247e + 5

0.7 0.7 0.7

WE-0 UIE- I

4.262e - 4 4.255e - 4

713 x 86 738 X 79

1.245e + 5 2.410e + 6

0.7 0.6

CIE-0 CIE- I CIE-2

3.554e - 2 4.468e - 4 4.254e - 4

713 x 86 738 x 79 763 X 87

2.294e + 5 1.179e + 5 2.431e + 5

0.7 0.6 0.8

Table 4 has the results for plane wave scattering at ku = 1. The Sommerfeld condition yields the largest trace error. The first-order ABC of Bayliss and Turkel is an order of magnitude lower; the second-order ABC is two orders of magnitude lower. (In fact, it is limited by the least error possible for this mesh.) The trace error for the unconjugated infinite element with q = 0 appears to be at nearly the threshold limit, since the error for IJIE-1 is only slightly lower. (Though, note the twenty-fold increase in the condition number of UIE-1 over UIE-0.) The conjugated infinite element with q = 0 is equivalent to the Sommerfeld ABC. The error for CIE-1 is down by two orders of magnitude over CIE-0, and the error for CIE-2 is limited by the threshold error for the mesh. The results for plane wave scattering at ku = 5 are shown in Table 5. At this higher frequency the Sommerfeld and B-&T-l ABCs exhibit roughly the same accuracy. The B&T-2 condition is better than the two other ABCs by roughly a factor of 20. Exponential convergence is evident for the unconjugated elements, with the trace error for UIE-2 being limited by the threshold error. As in 2D, their condition numbers grow exponentially with order q. The error for the conjugated elements is also decreasing exponentially, but their condition numbers show no growth trend. The solution times for the infinite elements is increasing with q as more linear equations are added to the system and the bandwidth grows. In Table 6 we record the results for plane wave scattering at ku = 10. At the highest frequency there is little distinction between the Sommerfeld and the first-order Bayliss and Turkel conditions. The error for B&T-2 is only a factor of 10 smaller than the other ABCs, but still several orders of magnitude larger than the threshold error that can be reached using infinite elements. Once again, the exponential convergence of the infinite elements is evident (down to the threshold error), as is the exponential growth of the condition numbers for the unconjugated elements. For the spherical infinite elements we are also interested in the quality of the approximation for r Z= R, i.e. within the element itself. Figs. 5-7 are plots of the real part of the scaled exterior solution ,X as a function of R/r at the polar angle 0 = 0. (Recall that puex is the function that remains after factoring out from uex its Table 5 Plane wave scattered

by a sphere of radius a at ka = 5

Method

/la - aeall

System size

Cond. no.

CPU

Somm

B&T- 1 B&T-2

8.990e - 2 4.35Oe - 2 2.352e ~ 3

3369 X 140 3369 x 140 3369 x 140

2.352e + 4 2.345e + 4 3.414e + 4

8.2 8.2 8.4

UIE-0 UIE- 1 UIE-2

2.51 le - 3 4.713e - 5 4.638e - 5

3369 x 140 3490 x 174 3611 x 193

3.415e + 4 3.708e + 6 7.407e + 8

7.3 Il.4 13.7

CIE-0 CIE-I CIE-2 CIE-3 CIE4 CIE-5

7.429e 1.844e 1.795e 1.230e 7.147e 4.646e

3369 3490 3611 3732 3850 3974

7.267e 2.753e 5.640e 7.241e 1.360e I .249e

7.6 Il.4 14.6 21.7 28.5 37.9

-

1 1 2 3 5 5

x x x X x x

140 174 193 248 303 358

+ f + + + +

5 5 5 5 6 6

J.J. Shirron, I. Babushka I Comput. Methods Appl. Mech. Engrg. 164 (1998) 121-139

134 Table 6 Plane wave scattered Method

by a sphere of radius a at ka = 10 System size

IIU- Ue”ll

Cond. no.

CPU

Somm

1.367e -

1

6689 X 140

3.157e + 4

18

B&T- 1 B&T-2

1.156e - 1 1.214e - 2

6689 X 140 6689 X 140

3.157e + 4 3.193e + 4

18 18

UIE-0 UIE-1 UIE-2

1.233e - 2 1.774e - 4 3.938e - 5

6689 X 140 6930 X 174 7171 X 193

3.193e + 4 9.381e + 6 3.982e + 9

18 26 31

CIE-0 CIE- 1 CIE-2 CIE-3 CIE-4 CIE-5 CIE-6 CIE-7 CIE-8

1.243e 7.446e 3.149e 8.720e 1.505e 1.989e 2.016e 4.317e 3.908e

6689 6930 7171 7412 7653 7894 8135 8376 8617

2.923e 1.249e 1.181e 6.361e 8.349e 4.615e 2.103e 1.288e 6.150e

+ -

0

1 1 2 2 3 4 5 5

X 140 X 174 x 193 X 248 X 303 X 358 X 413 X 468 X 523

+ + + + + + + + +

6 6 6 5 5 5 6 6 5

18 25 29 42 56 73 94 120 148

Sealed Exterior Field: ka-1 -0.70

I

0.4

0.6 R/r

Fig. 5. Real part of scaled exterior solution at 0 = 0, ka = Fig. 6. Real part of scaled exterior

1.

solution at 8 = 0, ka = 5.

Sceld Extwlor Field: ka-10 0.0

R/r Fig. 7. Real part of scaled exterior solution at 0 = 0, ka = 10.

0.8

1.0

J.J. Shirron, I. Babushka I Cornput. Methods Appl. Mech. Engrg. 144 (1998) 121-139

135

asymptotic oscillation and decay.) These figures indicate that the situation is much the same as in 2D: the conjugated elements closely approximate the exact solution over R G r < m while the approximation diverges in the unconjugated elements as r -+ ~0. Note that the comparison of the conjugated and unconjugated elements for the same order 4 is misleading due to their different rates of convergence. In Fig. 7 we plot curves for both CIE-2 and CIE-4; the curve for CIE4 is nearly indistinguishable from the curve for the exact function.

8. Generalization

to prolate spheroidal

coordinates

One class of radiators and scatterers of interest to the Navy consists of bodies that are long and thin, like submarines. It would be extremely inefficient to use a spherical boundary to truncate the exterior domain for such a slender body. In that way, for even the closest fitting sphere, an enormous volume of fluid would need to be modeled. Upon discretization, this would lead to extremely large systems of linear equations with correspondingly large bandwidths. The optimal situation is instead to truncate the exterior with a prolate spheroid that is tailored to fit the body as closely as possible. This idea was first suggested by Burnett [24], who has recently extended the idea by proposing an ellipsoidal truncation boundary [25]. In this way the finite element domain is kept as small as practical. One may then apply an ABC developed for spheroidal boundaries (see [33] and [13] for conditions on an ellipse in 2D), apply the DtN in spheroidal coordinates [43], or use spheroidal infinite elements [24,57]. In either of these strategies the size of the linear system that results from discretization would be a small fraction of the size obtained from its spherical counterpart. As we mentioned before, the DtN condition is not local. A consequence is that on discretization dense blocks are filled in an otherwise sparse system of linear equations generated by a FEM. These dense blocks increase the storage requirements and solution times for such systems. For these reasons we decided against implementing the DtN. In principle, one can derive a hierarchy of ABCs in prolate spheroidal coordinates analogous to those of Bayliss and Turkel in spherical coordinates. However, we have chosen not to do so based on the results in the previous section. The indication there was that the accuracy of ABCs decreased markedly as the frequency increased. To improve the accuracy one could move the truncation boundary farther out, but that would increase the computational cost. Or one could use a higher order condition. However, that would necessitate using a FEM with better than Co continuity. Instead, we have chosen to develop infinite element methods in prolate spheroidal coordinates. Let 2a be the distance between the foci of a prolate spheroid aligned with its long axis coinciding with the z axis in Cartesian coordinates and its center at the origin. The prolate spheroidal coordinates (r, 13,4) are related to the Cartesian coordinates through x = J;z-d’

sin 6 cos 4 ,

(51)

4’ = \iv’ - a2 sin 13sin 4,

(52)

z=t-cost?,

(53)

for r 3 a, 0 4 0 C T and 0 G 4 < HIT. Recall that the infinite element shape functions in spherical coordinates were given by an ansatz suggested by the asymptotic expansion of Atkinson and Wilcox. This is also how we construct the shape functions in prolate spheroidal coordinates. The far-field expansion is now [24,46]

(54) which is valid outside the smallest prolate spheroid enclosing then of the form

the body. The infinite element shape functions

are

(55) where ,u tensor product of polynomials of degree q in t = R/r by angular shape functions of degree p. Note that this has exactly the same form as the spherical shape functions; all that has changed is the geometric meaning of

J.J.

136

Shirron,

I. Babushka

I Comput.

Methods

Fig. 8. Hemispherically-capped

Appl.

Mech.

Engrg.

164

(1998)

121-139

cylinder in a prolate spheroid.

the coordinates. Not surprisingly, the conjugated and unconjugated versions of the prolate spheroidal infinite element inherit convergence and stability behavior that is similar to that of their spherical counterparts. To illustrate the efficacy of the prolate spheroidal infinite elements for problems in the exterior of long, slender bodies we consider the case of a radiating hemispherically-capped cylinder. We shall specify the Neumann data on the body such that it radiates as an acoustic point source with pole at the origin. (In this case, if the finite/infinite element boundary were spherical, infinite elements of order q = 0 would be sufficient for approximating the exact solution.) This problem is axisymmetric so we explicitly remove the dependence of the solution on the azimuthal angle (6.We take the cylinder to have radius b = 1, length L = 10 and we truncate the el ment domain with a prolate spheroid of radius R = 6.9 and inter-focal length 2a such that S = nY= R2 - a2 = 2. Fig. 8 shows a coarse mesh (not the one used for computation) for the two-dimensional finite element domain. In the actual computations the domain is meshed quasi-uniformly with quadrilateral elements for which the longest edge is 6 = 0.1. We show results only at the frequency kb = 5, where the acoustic wavelength A = 1.26, so the mesh has roughly 12 elements per wavelength. The finite elements are cubic (p = 3) and the infinite elements have cubic angular variation and are of order q in the reciprocal radius. Table 7 records the results for the conjugated and unconjugated infinite elements of various orders. The first column is the error in the solution on the surface of the capped cylinder. For the UIE we note the reasonably fast convergence of the trace error up to q = 2, beyond which point the instability of the method prevents further exponential convergence. The instability is evident in the exponentially increasing condition numbers for the FEM/UIE systems. For the conjugated version we see very slow convergence (but still exponential) with increasing order q. Unlike the UIE method, the CIE method is stable (as evidenced by the condition numbers of the FEM /CIE systems) and the trace error can be driven down to the mesh threshold by taking ever larger orders q. As can be seen, the CIE converges to an error that is two orders of magnitude smaller than the error for the UIE. But this is achieved only at a substantial computational cost: the solution time for CIE-30 is more than 17 times longer than for UIE-2 and the memory requirement is 13 times greater. Nevertheless, if high accuracy is required only the CIE method can be counted on to deliver it. In defense of the UIE method, it is worth mentioning that for UIE-3 the maximum relative error in the solution over the surface of the capped cylinder is only 3.6%, which may be of sufficient accuracy for many purposes.

Table 1 Radiation from a 10: 1 hemispherically-capped Method

IIU- 0

cylinder at kb = 5

System size

Cond. no.

CPU

UIE-0

8.656e - 1

6175 X 78

3.156e + 04

12

UIE- 1

1.976e - 1

7226

5.255e - 2

1617 X 85

+ 07 + 10

14

UIE-2

9.827e 3.846e

UIE-3

5.201e - 2

8128 x 97

1.379e + 12

20

6115 X 78

5.835e + 6

12

3.214e + 6

26

3.855e + 6

50

CIE-0

1.016e + 1

CIE-5

3.753e

+ 0

X 81

9030 11285

17

CIE-10

1.196e+O

CIE-15

5.088e - 1

x 121 x 181 13540 X 241

CIE-20

3.236e - 2

15795 x 301

5.094e + 7

140

9.233e + 6

89

CIE-25

2.144.e - 3

18050 X 361

1.669e+l

209

CIE-30

5.367e

20305 X 421

1.594e + 8

297

- 4

J.J. Shirron, 1. BabuZka 1 Comput. Methods Appl. Mech. Engrg. 164 (1998) IZI- 139

137

9. Conclusions We have compared various approximate boundary conditions and infinite element methods for the exterior Helmholtz equation for two simple cases where exact solutions are available: plane wave scattering by a rigid circle and a rigid sphere. We measured the error in trace on f, since if ulr were computed accurately the field at any point x E R, could be computed with corresponding accuracy via the Helmholtz integral (3 1). Although the results discussed in this paper were for a Neumann boundary condition, the results for Dirichlet and Robin boundary conditions on r (not reported here) are qualitatively the same. The results in 3D were very similar to those obtained in 2D. Our meshes were constructed in a special way. At each frequency ka we placed the truncation surface 1; one wavelength out from f. We then constructed a uniform mesh with 16 elements in the radial direction and 8ku (in 3D) or 16ku (in 2D) elements in the angular direction. The finite elements were cubic; the infinite elements were cubic in the angular direction and of variable order q in the reciprocal radius. In general, we found that the trace errors for the ABCs increased with ku. For low ku the error from the better conditions was comparable to that obtained using infinite elements. In particular, the error from the second-order Bayliss and Turkel condition was limited by the least error possible for the mesh. For higher values of ku the distinctions between the ABCs diminished. The fact that a condition was first- or second-order had less bearing on its accuracy, the exception being the Sommerfeld condition which uniformly yielded the largest errors. At ku = 10 the errors from the ABCs were several orders of magnitude above the threshold error for the mesh, while the infinite elements of high enough order were so accurate as to be limited by the threshold. The ABCs did have the advantages of yielding the smallest systems of linear equations (which were complex symmetric) and the lowest condition numbers for those systems. For both versions of the infinite element the trace error decreased exponentially with increasing order q. However, the rate of convergence for the conjugated elements was considerably slower than for the unconjugated elements. Still, at each frequency it was possible to take q large enough that the error in the computed solution for the -conjugated elements was limited only by the threshold error for the mesh. Such convergence was not always possible for the unconjugated elements, at least not in prolate spheroidal coordinates. The condition numbers of the linear systems using conjugated elements were relatively insensitive to the order q, while they grew exponentially for the systems using unconjugated elements. For q > 0 the infinite element systems were larger and had higher bandwidths than the ABC systems, which meant longer solution times. For the conjugated elements in particular, the increase in the solution time was not trivial; in the worst cases it was a factor of 5 or more. When a solution is approximated using finite elements and an ABC, one obtains information only on a finite computational domain. To compute the solution outside that domain one must resort to the Helmholtz integral. By contrast, when a solution is sought using infinite elements, the computational domain is itself unbounded and an approximation is available in both the near and far fields. We have shown that the quality of the approximation differs greatly depending on which version of the infinite element was used. The approximation computed with conjugated elements converged to the exact solution for all R S r < ~0 as the order q was increased. However, the approximation obtained using unconjugated elements diverged from the exact solution as r + cc. Consequently, no reliable far-field information could be extracted from the unconjugated elements, while very accurate information could be had from the conjugated elements (but that computation would be more costly). The inability of the unconjugated infinite elements to converge to the exact solution over R 4 r c: m is a consequence of the instability of the bilinear form on which they are based. This instability precludes convergence of the UIE method and manifests itself in the ill conditioning of the FEM/UIE linear systems.

References [l] M. Abramowitz and IA. Stegun, Handbook of mathematical functions, U.S. Department of Commerce, 1972. [2] S. Amino, P.J. Harris and D.T. Wilton, Coupled Boundary and Finite Element Methods for the Solution of the Dynamic Fluid-Structure Interaction Problem (Springer-Verlag. Berlin, 1992). [3] R.J. Astley, Wave envelope and infinite elements for acoustical radiation, Int. J. Numer. Methods Fluids 3 (1983) 507-526.

138

J.J. Shirron, 1. Babushka I Comput. Methods Appl. Mech. Engrg.

164 (1998)

121-139

[4] R.J. Astley, FE mode-matching schemes for the exterior Helmholtz problem and their relationship to the FE-DtN approach, Comm. Numer. Methods Engrg. 12(4) (1996) 257-267. [5] R.J. Astley, J.P. Coyette and G.J. Macaulay, Mapped wave envelope elements for acoustical radiation and scattering, J. Sound Vib. 170 (1994) 97-l 18. [6] 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. [7] R.J. Astley and W. Eversman, Finite element formulations for acoustic radiation, J. Sound Vib. 88(1983) 47-64. [S] R.J. Astley and W. Eversman, Wave envelope and infinite element schemes for fan noise radiation from turbofan inlets. Am. Instit. Aeronaut. Astronaut. J. 22 (1984) 1719-1726. [9] FS? Atkinson, On Sommerlled’s ‘radiation condition’, The Philosophical Mag. 40 (1949) 645-651. [lo] I. Babushka, Error-bounds for finite element method, Numer. Math. 16 (1971) 322-333. [ 1 l] A. Bayliss, M. Gunzburger and E. Turkel, Solution of elliptic equations in exterior regions, SIAM J. Appl. Math. 42 (1982) 430-451. [12] A. Bayliss and E. Turkel, Radiation boundary conditions for wave-like equations, Comm. Pure Appl. Math. 33 (1980) 707-725. [13] G. Benporat and D. Givoli, Solution of unbounded domain problems using elliptic artificial boundaries, Comm. Numer. Methods Engrg.ll(9) (1995) 734-741. [14] P. Bettess and O.C. Zienkiewicz, Diffraction and refraction of surface waves using finite and infinite elements, Int. J. Numer. Methods Engrg. 11 (1977) 1271-1290. [15] P. Bettess, Infinite elements, Int. J. Numer. Methods Engrg. 11 (1977) 53-64. [16] P. Bettess, More on infinite elements, Int. J. Numer. Methods Engrg. 15 (1980) 1613-1626. [17] P. Bettess, Operation counts for boundary integral and finite element methods, Int. J. Numer. Methods Engrg. 17 (1981). [18] P. Bettess, A simple wave envelope element example, Comm. Appl. Numer. Methods 3 (1987) 77-80. [ 191 P. Bettess, Infinite Elements (Penshaw, Sunderland, UK, 1992). [20] P. Bettess and J.C. Bettess, Infinite elements for dynamic problems: Part 1, Engrg. Comput. 8 (1991) 99-124. [21] P. Bettess and J.C. Bettess, Infinite elements for dynamic problems: Part 2, Engrg. Comput. 8 (1991) 125-151. [22] H. Brakhage and P. Werner, Uber das Dirichletsche Aussenraumproblem fiir die Helmholtzsche Schwingungsgleichung, Arch. Math. 16 (1965) 325-329. [23] D.S. Burnett, 3-D structural acoustic modeling using a new high accuracy acoustic infinite element: A cost comparison with the BEM, Private communication of prepublication results, in preparation. [24] D.S. Burnett, A three-dimensional acoustic infinite element based on a prolate spheroidal multipole expansion, J. Acoust. Sot. Am, 96 (1994) 2798-2816. [25] D.S. Burnett and R.L. Holford, An ellipsoidal infinite element for 3D radiation and scattering, Comput. Methods Appl. Mech. Engrg., submitted. [26] A.J. Burton and G.F. Miller, The application of integral equation methods to the numerical solution of some exterior boundary value problems, Proc. Roy. Sot. Lond. A323 (1971) 201-210. [27] C.C. Chien, H. Rajiyah and S.N. Atluri, An effective method for solving the hypersingular integral equations in 3-D acoustics, J. Acoust. Sot. Am. 88 (1990) 918-937. [28] D. Colton and R. Kress, Integral equation methods in scattering theory, in: Pure and Applied Mathematics (John Wiley & Sons, Inc., New York, 1983). [29] L. Cremers and K.R. Fyfe, Numerical methods for solving acoustic radiation problems, Can. Acoust. 19(4) (1991) 29-30. [30] L. Cremers, K.R. Fyfe and J.P. Coyette, A variable order infinite acoustic wave envelope element, J. Sound Vib. 171 (1994) 483-508. [3 11 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 (1995) 2028-2040. [32] I. Patlashenko, D. Givoli and J.B. Keller, High-order boundary conditions and finite elements for infinite domains, Comput. Methods Appl. Mech. Engrg. 143 (1997) 13-39. [33] A.F. Peterson, D.B. Meade, W. Slade and K.J. Webb, Comparison of local radiation boundary conditions for the scalar Helmholtz equation with general boundary shapes, IEEE Trans. Antennas Propagation 43( 1) (1995) 6-10. [34] L. Demkowicz and K. Gerdes, Convergence of the infinite element methods for the Helmholtz equation, Technical Report 95-07, The Texas Institute for Computational and Applied Mathematics, University of Texas at Austin, May 1995. [35] L. Demkowicz and F. Ihlenburg, Analysis of a coupled finite-infinite element method for exterior Helmholtz problems, Technical Report 96-l 1, The Texas Institute for Computational and Applied Mathematics, University of Texas at Austin, November 1996. [36] J.A. Eaton and B.A. Regan, Application of the infinite element method to acoustic scattering problems, AIAA J. 34( 1) (1996) 29-34. [37] B. Engquist and A. Majda, Absorbing boundary conditions for the numerical simulation of waves, Math. Comput. 31 (1977) 629-65 1. [38] B. Engquist and A. Majda, Radiation boundary conditions for acoustic and elastic wave calculations, Comm. Pure Appl. Math. 32 (1979) 313-357. [39] E. Anderson et al., LAPACK User’s Guide. SIAM, Philadelphia, 1992. [40] K. Gerdes and L. Demkowicz, Solution of 3D Laplace and Helmholtz equations in exterior domains using hp-infinite elements, Comput. Methods Appl. Mech. Engrg. 137(3-4) (1997) 239-273. [41] D. Givoli and J.B. Keller, Special finite elements for use with high-order boundary conditions, Comput. Methods Appl. Mech. Engrg. 119(3-4) (1994) 199-213. [42] D. Givoli, Non-reflecting boundary conditions, J. Comput. Phys. 94 (1991) l-29. [43] M.J. Grote and J.B. Keller, On nonreflecting boundary conditions, J. Comput. Phys. 122(2) (1995) 231-243. [44] I. Harari and D. Avraham, High-order finite element methods for acoustic problems, J. Comput. Acoust. 5( 1) (1997) 33-5 1. [45] I. Harari and T.J.R. Hughes, A cost comparison of boundary element and finite element methods for problems of time-harmonic acoustics, Comput. Methods Appl. Mech. Engrg. 97 (1992) 77-102.

J.J. Shirron, I. Babufka I Comput. Methods Appl. Mech. Engrg. [46] R.L. Holford, A multipole expansion [47] [48] [49] [50] [51] [52] [53] [54] [55]

164 (1998) IZI-139

for the acoustic field exterior to a prolate or oblate spheroid, AT&T Bell Laboratories,

139 Whippany,

NJ., J. Acoust. Sot. Am., to be submitted. D.S. Jones, Integral equations for the exterior acoustic problem, Quart. J. Mech. Appl. Math. 27 (1974) 129- 142. F. Kang, Finite element method and natural boundary reduction, in: Proc. Int. Congress of Mathematicians, Warsaw. Poland (1983) 1439-1453. S.V. Karp, A convergent ‘far-field’ expansion for two-dimensional radiation functions, Comm. Pure Appl. Math. 14 ( 1961 ;I 427-434. J.B. Keller and D. Givoli, Exact non-reflecting boundary conditions, J. Comput. Phys. 82 (1989) 172-192. N.N. Lebedev, Special Functions and Their Applications (Dover Publications, Inc., 1972). A.A. Oberai, M. Malhotra and PM. Pinsky, On the uyse of infinite elements and the DtN condition for solving large-scale problems in acoustics, to be submitted. R.W. Freund, M. Malhotra and PM. Pinsky, Iterative solution of multiple radiation and scattering problems in structural acoustics using a block quasi-minimal residual algorithm, Comput. Methods Appl. Mech. Engrg. l46( l-2) (1997) 173-196. 0.1. Panich, On the question of the solvability of the exterior boundary-value problems for the wave equation and Maxwell’s equations, Russian Math. Surveys 20 (1965) 221-226. I. Patlashenko and D. Givoli, Non-reflecting finite element schemes for three-dimensional acoustic waves, J. Comput. Acoust. 5(l) (1997) 955115.

[56] H.A. Schenk. Improved integral formulation for acoustic radiation problems, J. Acoust. Sot. Am. 44 (1968) 41-58. [57] J.J. Shirron, Solution of exterior Helmholtz problems using finite and infinite elements, Ph.D. Thesis, Department Mathematics, University of Maryland, 1995. [S8] A. Sommertled, Die Greensche funktion der Schwingungsgleichung, Jber. Deutsch. Math. Verein. 21 (1912) 309-353. [59] B. Szabo and I. BabuSka, Finite Element Analysis (John Wiley & Sons, Inc., New York, 1991). [60] F. Ursell, On the exterior problems of acoustics, Proc. Cambridge Philosophical Sot. 74 ( 1973) 117- 125. [6l] C.H. Wilcox, An expansion theorem for electromagnetic fields, Comm. Pure Appl. Math. 9 (1956) 115- 134. [62] C.H. Wilcox, A generalization of theorems of Rellich and Atkinson, Proc. Am. Math. Sot. 7 (1956) 271-276.

of Applied