- Email: [email protected]
Higher-order boundary infinite elements
Computer methods in applied mechanics and englneerlng Comput.
ELSEWER
Methods
Appl. Mech. Engrg.
164 (1998)
107-l 19
Higher-order boundary infinite elements Isaac Harari”‘“, Rami Shalom”, Paul E. Barboneb “Department of Solid Mechanics, Materials & Structures, Tel-Aviv University, 69978 Ramat Aviv, Israel ‘Department of Aerospace & Mechanical Engineering, Boston University, 110 Cummington Street, Boston, MA 02215. USA Received
23 July 1997; revised 23 November
1997
Abstract Higher-order approximations are implemented in a novel approach to infinite element formulations for exterior problems of timeharmonic acoustics. This approach is based on a functional which provides a general framework for domain-based computation of exterior problems. Two prominent features of this formulation simplify the task of discretization: the lack of integration over the unbounded domain, and weak enforcement of continuity between finite elements and infinite elements. Consequently, the infinite elements mesh the interface only and need not match the finite elements on the interface. Various infinite element approximations for two-dimensional configurations with circular interfaces are presented. Numerical results demonstrate the good performance of this scheme. 01998 Elsevier Science S.A. All rights reserved.
1. Introduction Boundary-value problems in exterior domains pose a unique challenge to computation, since the unbounded region is inappropriate for direct discretization. Numerical solutions to exterior problems for the Helmholtz equation, which governs problems of time-harmonic acoustic radiation and scattering, have been sought frequently by methods that are based on integral representations (see, e.g. the survey in [26]), relating quantities on the physical boundary. An alternative strategy is domain-based computation. Finite element methods are the numerical technique of choice for numerous classes of boundary-value problems, and are being used increasingly in structural acoustics (see e.g. [ 1,4,13,18,25,28,29] and the survey [ 161). In order to use domain-based discretization in an exterior problem, a bounded computational domain is formed by introducing an artificial boundary. Correct far-field behavior is then enforced either by specifying proper boundary conditions on this boundary, or by assuming interpolation with suitable behavior in the complement of the computational domain [Ill. The former approach is associated with non-reflecting boundary conditions, and in particular, the DtN method. This method was conceived as a general procedure for exterior boundary-value problems by Givoli and Keller [ 121, and is related to earlier work in acoustics [8,24]. The DtN method was recently shown [2] to be related to mode matching [7]. On the other hand, the term ‘infinite elements’ refers to a class of methods which is based on interpolation with suitable behavior in the complement of the computational domain [3,5,6,10,3 11. Global DtN boundary conditions can be very accurate, but all of the degrees of freedom on the artificial boundary are coupled, potentially increasing the cost of computation. In contrast, infinite elements retain the element-based data structure of finite elements, thus preserving the bandedness of the discrete equations, but typically at the cost of including extra unknowns in the formulation arising from discretization of the unbounded domain. In the following, a novel approach to the development of infinite element formulations for exterior problems
* Corresponding
author.
0045-7825 /98/$19.00 0 1998 Elsevier Science S.A. All rights reserved PII: SOO4S-7825(98)00049-B
I. Harari et al. I Comput. Methods Appl. Mech. Engrg. 164 (1998) 107-119
108
of time-harmonic acoustics [15] is employed. In Section 2 the exterior acoustics problem is decomposed by an artificial boundary into inner and outer fields. The inner field will be approximated by finite elements and the outer field by infinite elements. A simple functional weakly enforces continuity of the two fields. The functional is defined on the bounded inner domain and interface only, motivating the term boundary infinite elements. This formulation is exact, free of constraint degrees of freedom, and can be used to recover many well-known formulations developed previously [14], including the DtN boundary conditions. Criticism [9,20,21,23] of a similar approach [22] is examined in [15], leading to the conclusion that for our purposes (namely, enforcing flux continuity across the boundary) this procedure can lead to an effective and accurate numerical method. The weak form is modified in Section 3 to accommodate piecewise smooth representations. Linear and higher-order infinite element approximations are presented. Numerical results that validate the performance of the infinite elements are presented in Section 4.
2. Exterior boundary-value problem of acoustics Let $J?C lRd be a d-dimensional unbounded region. The boundary of 3, denoted by r, is internal and assumed piecewise smooth (Fig. 1). The outward unit vector normal to r’ is denoted by n. We assume that r admits the partition
r =<
U 4, where 4 II 4 = 0.
2.1. Strong form We consider a boundary-value problem related to acoustic radiation and scattering governed by the Helmholtz equation: find U: ,G%! + C, the spatial component of the acoustic pressure or velocity potential, such that -Zu=f u=g
au
z=
in%?
(1)
oncr ikh
liir’d-“‘2
(2)
on & ($-iku)=O
Here _% := Au + k2u is the Helmholtz operator, A is the Laplace operator and k E C is the wave number, Im k > 0; au/ dn : = Vu . n is the normal derivative and V is the gradient operator; i = J-1 is the imaginary unit; r is the distance from the origin; and f: % + C, g: r, -+ C and h: r, + C are the prescribed data.
Fig. 1. An unbounded
region with an internal boundary
1. Harari et al. I Comput. Methods Appl. Mech. l$grg.
164 (1998) 107-I 19
109
Eq. (4) is the Sommefeld radiation condition and allows only outgoing waves at infinity, The radiation condition requires that energy flux at infinity be positive, thereby guaranteeing that the solution to the boundary-value problem (l)-(4) is unique (see [27, pp. 296-2991 and [30, pp. 55-601 and references therein). Appropriate representation of this condition is crucial to the reliability of any numerical formulation of the problem (l)-(4). We now outline the development of a simple variational formulation that couples two parts of an acoustics problem in a partitioned domain by weakly enforcing continuity at the interface [ 151. 2.2. Partitioned problem The unbounded domain 3 is partitioned by a smooth artificial boundary IY,into a bounded inner domain 0’ and its unbounded outer complement a0 (Fig. 2), expressed analytically as
cB==niuno
(5)
where 0’ fl0’ = 0 and r, =z 02. The solution of the original boundary-value problem (l)-(4) field u’, where uiln,, = 0 and ~‘(~i = 0 so that u=
is decomposed into the inner field ui and outer
u 1’
on 0’ u0 on 0”
Continuity across f, is enforced weakly. We assume that f vanishes in the outer domain: .f=O
in 0’
2.3. Variational statement
(7) with weak continuity
The variational form of the boundary-value problem (l)-(4) is stated in terms of sets of trial solutions, Y’ and Y” for ui and u’, respectively. For Neumann problems (4 = Q)), Y’ = H’(L?). Otherwise, Dirichlet boundary conditions on < must also be satisfied by functions in 9’. The treatment of Y” is less conventional, namely u”(~~“=Oin~“,~tm_r(d-1)/2
($-ikilo)
=O>
Continuity across the artificial boundary is enforced weakly. Thus, we seek u’ E Y’ and u0 E Y0 that render II(u) stationary, where
Fig. 2. The partitioned
domain.
110
I. Harari et al. I Comput. Methods Appl. Mech. Engrg. 164 (1998) 107-119
1 Il(u)=zcc(u’,ui)-L(ui)+
a(w, u) =
I
0’
0 (
3,
1 Ui _ _2 Uo r, )
(VW *Vu - wk*u) da
(10)
(11) (12) Here, a(., *) and (e, *)TRare symmetric
bilinear
forms which are not inner products.
In this formulation there is integration only over the bounded domain and the artificial boundary. Thus, the outer field u” may be viewed as specifying boundary conditions on the artificial boundary r,, thereby defining a boundary-value problem for the inner field ui in the bounded domain ni. Continuity across the artificial boundary is weakly enforced by energy flux-like terms, reminiscent of the radiation condition that is being replaced. Specifically, the Euler-Lagrange equations of (9) are -2~’
=f
in 0’
ikh
on J,
(13)
dU’
z= U’ =
u*
on r,
aui au” on r, dn’ an”
-=--
(14) (15)
(16)
These provide satisfaction of the differential equation in Oi and boundary conditions on 4, and enforce continuity of the functions and their normal derivatives across the artificial boundary (n’ denotes the outward normal with respect to 0’). Previous investigations of the well-posedness of bounded-domain formulations for exterior problems [ 17,191 indicate that enforcing continuity of the functions and their normal derivatives is fitting. Admissible outer fields are sufficiently regular for II, namely, au”l&z” EH-“2(G) [15]. The stationary point of the functional (9) is obtained by setting its first variation equal to zero. This is equivalent to:
Here, WI E vi and w” E “Y” are the arbitrary variations of ui and u”, respectively. Accordingly, for Neumann problems “v^’= HI(@); otherwise, homogeneous counterparts of Dirichlet boundary conditions on c must also be satisfied by functions in “Yi. In the outer field “Y””= Y”. The formulation is symmetric since V w” E v” and VuOEYO
(18) The weighting functions in the weak form (17) are not conjugated, in accordance with the complex-valued functional (9). A similar weak form with conjugated weighting functions, as employed by the wave envelope procedure [3], can be derived directly from the Euler-Lagrange equations (13)-(16).
I. Hurari et al. I Comput. Methods Appl. Mech. Engrg. 164 (1998) 107-119
111
3. Piecewise outer field representations We now account for possible discontinuities across infinite element boundaries when u” is approximated by piecewise smooth functions. Therefore, we consider a partition of the outer domain LP into nsec nonoverlapping sectors R” with boundaries r”, s = 1, . . . , nbec. For example, Fig. 3 shows the partition in a two-dimensional configuration with a circular interface. We denote the union of sector interiors by
fin” =
u
0.’
(19)
r=l The interior
sector boundaries
are
(20) where 80” is the boundary of no. One side of rynt is designated as ‘ + ’ and the other side is ‘ - ’ (Fig. 4). Allowing rp,,, we define the jump across rynt as:
for discontinuities
across
(21) Here x+ and x are obtained by approaching ry,, from the + and - sides, respectively. Similarly, n + and n _ are unit normals to ry,, pointing out of the + and - sides, respectively. The jump is invariant to the arbitrary + and - designations.
3.1. Continuous, For continuous,
piecewise-smooth
functions
piecewise-smooth
Y” = ( u” 1U”E Co(C),
representations
of the outer field the set of trial solutions
234”= 0 in ho, l& r ‘d-i”2 (~_iku’)
Fig. 3. Partition
in a two-dimensional
configuration
=o>
with a circular
Fig. 4. Two sides of an interior sector boundary.
is modified
(22)
interface.
I. Harari et al. I Comput. Methods Appl. Mech. Engrg. 164 (1998) 107-119
112
and Y” = 9’“. The weak form must now account for -possible jumps in derivatives across ry,,,. Thus we seek ui E Y’ and u” E 9’” such that
(23)
The Euler-Lagrange
I-1 au” dn
= 0
equations are (13)-(16)
and
on rynt
(24)
enforcing continuity of normal derivatives across interior sector boundaries. In this case du”ldn” E H-“*(T” fl r,), s = 1, . . . , n set, as required on each sector. The formulation retains its symmetry since tl w” E V” and Vu”EYO
(25) 3.2. Injinite element approximation The formulation (23) is a generalization of (17) that accounts for possible discontinuities across infinite element boundaries. We may now employ infinite element basis functions that exactly satisfy the Helmholtz equation within each sector and are continuous across infinite element boundaries. This approximation is employed directly with Eqs. (23). In most cases, the circumferential behavior of these functions on the artificial boundary differs from piecewise polynomial variation of finite element shape functions. Weak enforcement of continuity across the artificial boundary is then performed naturally within the framework of our formulation (23). Consequently, the infinite element mesh need not be the restriction of the finite element mesh to the artificial boundary. This flexibility may offer a significant advantage in practical applications. To fix ideas consider the two-dimensional configuration with a circular interface shown in Fig. 3. The outer domain 0“ is partitioned into non-overlapping sectors 0S=]R,w[X]8,,8,+,[,
s=l,...,
nsec.
(26)
Here, % < !, + , and 4, Let+ , = 8, + 2n. A typical sector is shown in Fig. 5. The infinite element approximation is based on interpolation of values on the arc r = R, 0, < 8 < 13,+, .
3.2.1. Linear infinite elements The basic two-noded infinite element is based on the representation u”=N,u”(R,~,)+N,u”(R,~,+,),
0,<0
(27)
where N, =
N2 =
Hy’(kr)
CT+, - 8
Hb”(kR)
Af?,
Hy’(kr)
0 - 0,
Hb’ ‘(kR)
A@,
(28)
(29)
are shape functions that satisfy the homogeneous Helmholtz equation within the sector, and AO,= OS+,- 0,. These functions exhibit outgoing oscillatory behavior in the radial direction and linear circumferential variation. The expressions required for infinite element arrays are simply
I. Harari et al. I Comput. Methods Appl. Mech. Engrg. I64 (1998) 107-119
Fig. 5. A typical sector in a two-dimensional
configuration
with a circular
113
interface.
(30)
N,V’C0) =
e- e,
(31)
Aey
and f-$
(R, 0) = -k
$&R,t?)=
-k
H;“‘(kR)
$+, - c9
Hb’ ‘(kR)
(32)
A@7
Hb”‘(kR)
28 - 0 H;’ ‘(kR) 4
(33)
The prime denotes differentiation with respect to the argument. The element array that is obtained from the term -(&v”/&r”, uO)~, in the second equation in (23). weakly enforcing continuity of the inner and outer fields across the interface, is proportional to the standard mass matrix of a linear line element kR
H;“‘(kR)-4 H;' “(kR)
6
2
1
[ 1
2
1
(34)
The element array that is obtained from the jump term (w”, [~3u”ldn]),~“~ in the second equation in (23), weakly enforcing derivative continuity of the outer field across the interior sector boundaries, is proportional to the standard stiffness of a linear element
z$-; The coefficient
-2
(35)
cyOOis a special case of a term that arises in higher-order
ffi HI'‘(kr)H’,“(kr) CY ylr = i R Hl”(kR)H;‘(kR)
(jr T
approximation
(36)
Integration of the Hankel functions is performed analytically, prior to computation. In subsequent implementations the unequal-order coefficient need not be integrated since it is evaluated by symmetry considerations, namely that the sum of (34) and (35) is symmetric.
I. Harari et al. I Comput. Methods Appl. Mech. Engrg. 164 (1998) 107-I 19
114
The evaluation of the integrals, when needed, is based on properties of U,(Z), solutions of the homogeneous Bessel equation of order v LJ?$u,= 0
(37)
%?,u, := z(zui>’ + (2’ - V2)U,
(38)
where
It follows from U,.s2l”U” -
u,~wup =0
(39)
that for Hankel functions, as solutions of the Bessel equation, ff v+ =p
kR
Hl”‘(kR)
p2 - v2
H;‘(kR)
HI”‘(kR) -
Hl”(kR)
(40)
1
In practice, the differentiated terms are eliminated by standard formulas for Bessel functions to simplify the implementation. Equal-order coefficients are evaluated by approximating the limit. (41)
at small values of V.Numerical experience indicates that this limit converges without difficulty for arithmetic of sufficient precision. The dependence of cyOO, the coefficient needed for the linear elements, on kR is shown in Fig. 6. The coefficient is well behaved and monotonically decreasing, indicating the diminishing importance of the jump terms for a given problem, as the location of the artificial boundary is moved outwards. As previously noted, the infinite element mesh need not match the finite element mesh on the interface r,. If linear infinite elements coincide with linear finite elements then the coupling terms in (23) are simply
real imaginary
Fig. 6. The coefficient
required
______________
for two-noded
infinite elements.
I. Harari et al. I Comput. Methods Appl. Mech. Engrg. 164 (1998) 107-119
proportional to (34). Otherwise, more elaborate flexibility offered by the formulation.
treatment
of the coupling
terms is required
115
for the added
3.2.2. High-order injinite elements Several extensions to higher-order approximation are possible. The basic interpolation in (27) can be enriched by a bubble Hx’(kr) sin(v,(8 - 19,)) where V, = nnlAt9,. The added coefficient in the approximation is a generalized degree of freedom rather than a nodal value. Regarding radial integration, only CZ,,,,is required as in the case of the linear element. Since 4, “I is associated with a bubble function in the circumferential direction it need not be evaluated. The matrix equations are statically condensed back to the size of the original two-noded element so that no additional computational effort is required in the solution. Alternatively, the outer field may be based on a truncation of the series representation
u” =
2
,?=0
In particular,
Hli,‘(k r )( a, COS(V,(~ - k),)) + b, sin(v,,(o - 0,)))
a three-degree-of-freedom H;“(kr)
1 N’ = 5
(
H;“(kR)
combines
the interpolatory
shape functions
HI,: ‘(kr) +
Hb’‘(kr) Hb”(kR)
element
(42)
Hj,; ‘(kR)
COO, (0 - es ))
(43)
cos(v, (0 - e,>>
(44)
HL: ‘(kr)
- H1,:‘(kR)
with the bubble Hy '(kr) sin(v, (0 - 19,))~again associated with a generalized degree of freedom. Due to the orthogonality of these functions, terms multiplying equal-order radial integrals aPI* vanish, so that no radial integration whatsoever is required. Again, static condensation of internal degrees of freedom is employed.
4. Numerical results: Radiation from a sector of a cylinder Numerical results are presented for exterior problems bounded internally by infinite circular cylinders of radius a. Soft boundary conditions are specified on the wet surface to represent a pressure-release cylinder. We consider the non-uniform radiation from an infinite circular cylinder with a constant inhomogeneous value on an arc (-a < ~9< cu) and vanishing elsewhere, so that there are two points of discontinuity in the boundary data. The normalized analytical solution to this problem for a cylinder of radius a is
(45) The prime on the sum indicates that the first term is halved. For low wave numbers this solution is relatively uniform in the circumferential direction. The directionality of the solution grows as the wave number is increased, and the solution becomes attenuated at the side of the cylinder opposite the radiating element. A circular interface r, is located at R = 2.5~ Finite elements are used in the inner domain LJi and infinite elements in the outer domain no. In the following examples the mesh of infinite elements (linear unless otherwise noted) coincides with the finite element mesh on the interface (although this is not required by the formulation). The inner domain ni is discretized by 3 X 36 bilinear quadrilateral finite elements (Fig. 7). We examine the problem with a geometrically non-dimensionalized wave number ka = 1 (the wavelength is about three times the diameter of the cylinder and four times the width of the domain), and a resolution of approximately 13 nodal points per wavelength in the finite element mesh.
4.1. Linear elements We select (Y= 20”. Fig. 8 shows the real part of the analytical and numerical solutions. The low-amplitude oscillations of the analytical solution in the vicinity of the wet surface are merely an artifact of the truncated
I. Harari et al. I Comput. Methods Appl. Mech. Engrg. 164 (1998) 107-I 19
116
Fig. 7. The computational domain quadrilateral finite elements.
exterior
to a cylinder
of radius a, with a circular
interface
at R = 2Sa, discretized
by 3 X 36 linear
series representation (60 terms) of the discontinuity in the boundary data, and are not relevant to the validation of the numerical results. The numerical solution captures the essential physics of the problem, without visible reflection from the interface. The error relative to the interpolation of the analytic solution, measured in the energy norm, is
(46) Fig. 9 depicts the solutions along the interface. The performance of the infinite elements is again quite good.
Fig. 8. Radiation denote numerical
from a sector of a cylinder, solution).
Fig. 9. Radiation
ka = 1, - 13 points per wavelength
from a sector of a cylinder,
(dashed contours denote analytical
ka = 1, -13 points per wavelength,
along the interface
solution, solid contours R = 2.5a.
I. Harari et al. I Comput. Methods Appl. Mech. Engrg. 164 (1998) 107- I19
Fig. 10. Radiation from a sector of a cylinder, ka = 1, - 13 points per wavelength, analytical solution, solid contours denote numerical solution). Fig.
bubble-enriched
11. Radiation from a sector of a cylinder, ka = 1, -13 points per wavelength,
117
infinite elements (dashed contours denote
bubble-enriched
infinite elements,
along the interface
R = 2Sa.
4.2. Bubble-enriched
elements
Here, elements enriched by the sine bubble are employed. enriched numerical solutions. The relative error
bh- 4, =
1o
9o/ .
,
,,
Fig. 10 shows the real part of the analytical
and
(47)
0
Iu’ I I
‘.._ :., ‘.., ‘... ‘\
j
;: ,.:
/’
,:’
j’ ,I’ :
t
,%.,
Fig. 12. Radiation from a sector of a cylinder, ka = 1, - 13 points per wavelength, analytical solution, solid contours denote numerical solution). Fig. 13. Radiation R = 2Sa.
from a sector of a cylinder,
ka = 1, -13
trigonometric
points per wavelength,
infinite elements (dashed contours
trigonometric
infinite elements,
denote
along the interface
118
I. Harari et al. I Comput. Methods Appl. Mech. Engrg. 164 (1998) 107-119
shows improvement over the linear infinite element. Fig. 11 depicts the solutions along the interface. The improvement over the linear infinite elements is apparent. 4.3. Trigonometric-based
elements
Here, three-degree-of-freedom elements based on a truncation of the trigonometric series (42) are employed. Fig. 12 shows the real part of the analytical and numerical solutions. The relative error (48) is a degradation from the linear infinite element. Fig. 13 depicts the solutions along the interface.
5. Conclusions
In this work, higher-order approximations are implemented in a novel method for developing infinite elements for exterior problems of time-harmonic acoustics. This technology is based on a simplified functional (9) for the partitioned problem, weakly enforcing continuity at the interface without additional degrees of freedom. The infinite elements mesh only the outer boundary of the finite element domain and need not match the finite elements on the interface. The implementation of these infinite elements is straightforward. Radial integration, when required, is performed analytically. Two approaches to higher-order approximation are presented: one is based on enrichment of low-order elements, and the other on the truncation of a trigonometric series. The second approach is free of radial integration. Infinite elements with three degrees of freedom developed by both approaches are examined. In both cases, static condensation is employed so that the computational effort required in the solution of the higher-order elements is no greater than for lower-order elements. Numerical results validate the improved performance of the bubble-enriched elements and indicate a need to employ more terms in the trigonometric series-based elements.
Acknowledgment
This research was supported by the U.S. Office of Naval Research. The authors wish to thank Parama Barai for calculation of error norms presented herein.
References [l] N.N. Abboud and P.M. Pinsky, Finite element dispersion analysis for the three-dimensional second-order scalar wave equation, Int. J. Numer. Methods Engrg. 35(6) (1992) 1183-1218. [2] 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. [3] R.J. Astley, G.J. Macaulay and J.-P. Coyette, Mapped wave envelope elements for acoustical radiation and scattering, J. Sound Vib. 170(l) (1994) 97-118. [4] I. Babulka, F. Ihlenburg, E.T. Paik and S.A. Sauter, A generalized finite element method for solving the Helmholtz equation in two dimensions with minimal pollution, Comput. Methods Appl. Mech. Engrg. 128(3-4) (1995) 325-359. [5] P. Bettess, Infinite elements, Int. J. Numer. Methods Engrg. 11(l) (1977) 53-64. [6] D.S. Burnett, A three-dimensional acoustic infinite element based on a prolate spheroidal multipole expansion, J. Acoust. Sot. Am. 96(5) (1994) 2798-2816. [7] H.S. Chen and CC. Mei, Oscillations and wave forces in a man-made harbor in the open sea, in: Tenth Symposium on Naval Hydrodynamics, pp. 573-594, Cambridge, Massachusetts, 1974. Office of Naval Research, Department of the Navy. [8] K. Feng, Asymptotic radiation conditions for reduced wave equation, J. Comput. Math. 2(2) (1984) 130-138. [9] R.H. Gallagher, E. Haugeneder and H.A. Mang, Limitations on methods of handling incompatible basis functions in boundary value problems, Trans. Am. Nucl. Sot. 33 (1979) 326-327. [IO] 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) (1996) 239-273.
1. Harari et al. I Comput. Methods Appl. Mech. Engrg. 164 (1998) 107-119 [l I] D. Givoli, Numerical methods for problems Publishing Co., Amsterdam, 1992).
in infinite domains,
in: Vol. 33 Studies
in Applied
Mechanics
119 (Elsevier
Scientific
[l2] D. Givoli and J.B. Keller, A finite element method for large domains, Comput. Methods Appl. Mech. Engrg. 76( 1) (1989) 41-66. [I31 K. Grosh and PM. Pinsky, Complex wave-number dispersion analysis of Galerkin and Galerkin least squares methods for fluid-loaded plates, Comput. Methods Appl. Mech. Engrg. I 13( 1-2) (1994) 67-98. [I41 I. Harari, P.E. Barbone and J.M. Montgomery Finite element formulations for exterior problems: Application to hybrid methods, non-reflecting boundary conditions and infinite elements, Int. J. Numer. Methods Engrg. 40( 15) (1997) 2791-2805. [15] I. Harari, P.E. Barbone, M. Slavutin and R. Shalom, Boundary infinite elements for the Helmholtz equation in exterior domains, Int. J. Numer. Methods Engrg. 41(6) (1998) 1105-I 131. [I61 I. Harari, methods [ 171 I. Harari domains,
K. Grosh, T.J.R. Hughes, M. Malhotra, PM. Pinsky, J.R. Stewart and L.L. Thompson, Recent developments in finite element for structural acoustics, Arch. Comput. Methods Engrg. 3(2-3) (1996) 131-309. and T.J.R. Hughes, Analysis of continuous formulations underlying the computation of time-harmonic acoustics in exterior Comput. Methods Appl. Mech. Engrg. 97(l) (1992) 103-124.
[ 181 I. Harari and T.J.R. Hughes, Galerkinlleast-squares [19] [ZO] [2l] [22] [23] [24]
finite element methods for the reduced wave equation with nonreflecting boundary conditions in unbounded domains, Comput. Methods Appl. Mech. Engrg. 98(3) (1992) 41 l-454. 1. Harari and T.J.R. Hughes, Studies of domain-based formulations for computing exterior problems of acoustics, Int. J. Numer. Methods Engrg. 37( 17) (1994) 2935-2950. E. Haugeneder and H.A. Mang, On an improper modification of a variational principle for finite element plate analysis, Z. Angew. Math. Mech. 59( 1I) (1979) 637-640. E. Haugeneder and H.A. Mang, Admissible and inadmissible simplifications of variational methods in finite element anal,ysis, in: S. Nemat-Nasser, ed., Variational Methods in the Mechanics of Solids, (Pergamon Press, Oxford, 1980) 187-198. S.W. Key, A specialization of Jones’ generalization of the direct-stiffness method of structural analysis, AIAA J. 9(5) (197111 984-985. H.A. Mang and R.H. Gallagher, A critical assessment of the simplified hybrid displacement method, Int. J. Numer. Methods Engrg. 1l(1) (1977) 145-167. M. Masmoudi, Numerical solution for exterior problems, Numer. Math. 51(l) (1987) 87-101.
[25] A. Oberai and P.M. Pinsky, A multiscale finite element method for the Helmholtz equation, Comput. Methods Appl. Mech. Engrg., 154(3-4) (1998) 281-297. [25] RI? Shaw, Integral equation methods in acoustics, in: C.A. Brebbia, ed., Boundary Elements X, Vol. 4 (Computational Mechanics Publications, Southampton, 1988) 22 I-244. [27] I. Stakgold, Boundary Value Problems of Mathematical Physics, Vol. II (Macmillan, New York, 1968). [28] J.R. Stewart and T.J.R. Hughes, h-adaptive finite element computation of time-harmonic exterior acoustics problems in two dimensions, Comput. Methods Appl. Mech. Engrg. 146( l-2) (1997) 65-89. [29] L.L. Thompson and P.M. Pinsky, A Galerkin least-squares finite element method for the two-dimensional Helmholtz equation, Int. J. Numer. Methods Engrg. 38(3) (1995) 371-397. [30] C.H. Wilcox, Scattering Theory for the d’Alembert Equation in Exterior Domains (Springer-Verlag. Berlin, 1975). [31] 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(7) (1985) 1229-1251.
ELSEWER
Methods
Appl. Mech. Engrg.
164 (1998)
107-l 19
Higher-order boundary infinite elements Isaac Harari”‘“, Rami Shalom”, Paul E. Barboneb “Department of Solid Mechanics, Materials & Structures, Tel-Aviv University, 69978 Ramat Aviv, Israel ‘Department of Aerospace & Mechanical Engineering, Boston University, 110 Cummington Street, Boston, MA 02215. USA Received
23 July 1997; revised 23 November
1997
Abstract Higher-order approximations are implemented in a novel approach to infinite element formulations for exterior problems of timeharmonic acoustics. This approach is based on a functional which provides a general framework for domain-based computation of exterior problems. Two prominent features of this formulation simplify the task of discretization: the lack of integration over the unbounded domain, and weak enforcement of continuity between finite elements and infinite elements. Consequently, the infinite elements mesh the interface only and need not match the finite elements on the interface. Various infinite element approximations for two-dimensional configurations with circular interfaces are presented. Numerical results demonstrate the good performance of this scheme. 01998 Elsevier Science S.A. All rights reserved.
1. Introduction Boundary-value problems in exterior domains pose a unique challenge to computation, since the unbounded region is inappropriate for direct discretization. Numerical solutions to exterior problems for the Helmholtz equation, which governs problems of time-harmonic acoustic radiation and scattering, have been sought frequently by methods that are based on integral representations (see, e.g. the survey in [26]), relating quantities on the physical boundary. An alternative strategy is domain-based computation. Finite element methods are the numerical technique of choice for numerous classes of boundary-value problems, and are being used increasingly in structural acoustics (see e.g. [ 1,4,13,18,25,28,29] and the survey [ 161). In order to use domain-based discretization in an exterior problem, a bounded computational domain is formed by introducing an artificial boundary. Correct far-field behavior is then enforced either by specifying proper boundary conditions on this boundary, or by assuming interpolation with suitable behavior in the complement of the computational domain [Ill. The former approach is associated with non-reflecting boundary conditions, and in particular, the DtN method. This method was conceived as a general procedure for exterior boundary-value problems by Givoli and Keller [ 121, and is related to earlier work in acoustics [8,24]. The DtN method was recently shown [2] to be related to mode matching [7]. On the other hand, the term ‘infinite elements’ refers to a class of methods which is based on interpolation with suitable behavior in the complement of the computational domain [3,5,6,10,3 11. Global DtN boundary conditions can be very accurate, but all of the degrees of freedom on the artificial boundary are coupled, potentially increasing the cost of computation. In contrast, infinite elements retain the element-based data structure of finite elements, thus preserving the bandedness of the discrete equations, but typically at the cost of including extra unknowns in the formulation arising from discretization of the unbounded domain. In the following, a novel approach to the development of infinite element formulations for exterior problems
* Corresponding
author.
0045-7825 /98/$19.00 0 1998 Elsevier Science S.A. All rights reserved PII: SOO4S-7825(98)00049-B
I. Harari et al. I Comput. Methods Appl. Mech. Engrg. 164 (1998) 107-119
108
of time-harmonic acoustics [15] is employed. In Section 2 the exterior acoustics problem is decomposed by an artificial boundary into inner and outer fields. The inner field will be approximated by finite elements and the outer field by infinite elements. A simple functional weakly enforces continuity of the two fields. The functional is defined on the bounded inner domain and interface only, motivating the term boundary infinite elements. This formulation is exact, free of constraint degrees of freedom, and can be used to recover many well-known formulations developed previously [14], including the DtN boundary conditions. Criticism [9,20,21,23] of a similar approach [22] is examined in [15], leading to the conclusion that for our purposes (namely, enforcing flux continuity across the boundary) this procedure can lead to an effective and accurate numerical method. The weak form is modified in Section 3 to accommodate piecewise smooth representations. Linear and higher-order infinite element approximations are presented. Numerical results that validate the performance of the infinite elements are presented in Section 4.
2. Exterior boundary-value problem of acoustics Let $J?C lRd be a d-dimensional unbounded region. The boundary of 3, denoted by r, is internal and assumed piecewise smooth (Fig. 1). The outward unit vector normal to r’ is denoted by n. We assume that r admits the partition
r =<
U 4, where 4 II 4 = 0.
2.1. Strong form We consider a boundary-value problem related to acoustic radiation and scattering governed by the Helmholtz equation: find U: ,G%! + C, the spatial component of the acoustic pressure or velocity potential, such that -Zu=f u=g
au
z=
in%?
(1)
oncr ikh
liir’d-“‘2
(2)
on & ($-iku)=O
Here _% := Au + k2u is the Helmholtz operator, A is the Laplace operator and k E C is the wave number, Im k > 0; au/ dn : = Vu . n is the normal derivative and V is the gradient operator; i = J-1 is the imaginary unit; r is the distance from the origin; and f: % + C, g: r, -+ C and h: r, + C are the prescribed data.
Fig. 1. An unbounded
region with an internal boundary
1. Harari et al. I Comput. Methods Appl. Mech. l$grg.
164 (1998) 107-I 19
109
Eq. (4) is the Sommefeld radiation condition and allows only outgoing waves at infinity, The radiation condition requires that energy flux at infinity be positive, thereby guaranteeing that the solution to the boundary-value problem (l)-(4) is unique (see [27, pp. 296-2991 and [30, pp. 55-601 and references therein). Appropriate representation of this condition is crucial to the reliability of any numerical formulation of the problem (l)-(4). We now outline the development of a simple variational formulation that couples two parts of an acoustics problem in a partitioned domain by weakly enforcing continuity at the interface [ 151. 2.2. Partitioned problem The unbounded domain 3 is partitioned by a smooth artificial boundary IY,into a bounded inner domain 0’ and its unbounded outer complement a0 (Fig. 2), expressed analytically as
cB==niuno
(5)
where 0’ fl0’ = 0 and r, =z 02. The solution of the original boundary-value problem (l)-(4) field u’, where uiln,, = 0 and ~‘(~i = 0 so that u=
is decomposed into the inner field ui and outer
u 1’
on 0’ u0 on 0”
Continuity across f, is enforced weakly. We assume that f vanishes in the outer domain: .f=O
in 0’
2.3. Variational statement
(7) with weak continuity
The variational form of the boundary-value problem (l)-(4) is stated in terms of sets of trial solutions, Y’ and Y” for ui and u’, respectively. For Neumann problems (4 = Q)), Y’ = H’(L?). Otherwise, Dirichlet boundary conditions on < must also be satisfied by functions in 9’. The treatment of Y” is less conventional, namely u”(~~“=Oin~“,~tm_r(d-1)/2
($-ikilo)
=O>
Continuity across the artificial boundary is enforced weakly. Thus, we seek u’ E Y’ and u0 E Y0 that render II(u) stationary, where
Fig. 2. The partitioned
domain.
110
I. Harari et al. I Comput. Methods Appl. Mech. Engrg. 164 (1998) 107-119
1 Il(u)=zcc(u’,ui)-L(ui)+
a(w, u) =
I
0’
0 (
3,
1 Ui _ _2 Uo r, )
(VW *Vu - wk*u) da
(10)
(11) (12) Here, a(., *) and (e, *)TRare symmetric
bilinear
forms which are not inner products.
In this formulation there is integration only over the bounded domain and the artificial boundary. Thus, the outer field u” may be viewed as specifying boundary conditions on the artificial boundary r,, thereby defining a boundary-value problem for the inner field ui in the bounded domain ni. Continuity across the artificial boundary is weakly enforced by energy flux-like terms, reminiscent of the radiation condition that is being replaced. Specifically, the Euler-Lagrange equations of (9) are -2~’
=f
in 0’
ikh
on J,
(13)
dU’
z= U’ =
u*
on r,
aui au” on r, dn’ an”
-=--
(14) (15)
(16)
These provide satisfaction of the differential equation in Oi and boundary conditions on 4, and enforce continuity of the functions and their normal derivatives across the artificial boundary (n’ denotes the outward normal with respect to 0’). Previous investigations of the well-posedness of bounded-domain formulations for exterior problems [ 17,191 indicate that enforcing continuity of the functions and their normal derivatives is fitting. Admissible outer fields are sufficiently regular for II, namely, au”l&z” EH-“2(G) [15]. The stationary point of the functional (9) is obtained by setting its first variation equal to zero. This is equivalent to:
Here, WI E vi and w” E “Y” are the arbitrary variations of ui and u”, respectively. Accordingly, for Neumann problems “v^’= HI(@); otherwise, homogeneous counterparts of Dirichlet boundary conditions on c must also be satisfied by functions in “Yi. In the outer field “Y””= Y”. The formulation is symmetric since V w” E v” and VuOEYO
(18) The weighting functions in the weak form (17) are not conjugated, in accordance with the complex-valued functional (9). A similar weak form with conjugated weighting functions, as employed by the wave envelope procedure [3], can be derived directly from the Euler-Lagrange equations (13)-(16).
I. Hurari et al. I Comput. Methods Appl. Mech. Engrg. 164 (1998) 107-119
111
3. Piecewise outer field representations We now account for possible discontinuities across infinite element boundaries when u” is approximated by piecewise smooth functions. Therefore, we consider a partition of the outer domain LP into nsec nonoverlapping sectors R” with boundaries r”, s = 1, . . . , nbec. For example, Fig. 3 shows the partition in a two-dimensional configuration with a circular interface. We denote the union of sector interiors by
fin” =
u
0.’
(19)
r=l The interior
sector boundaries
are
(20) where 80” is the boundary of no. One side of rynt is designated as ‘ + ’ and the other side is ‘ - ’ (Fig. 4). Allowing rp,,, we define the jump across rynt as:
for discontinuities
across
(21) Here x+ and x are obtained by approaching ry,, from the + and - sides, respectively. Similarly, n + and n _ are unit normals to ry,, pointing out of the + and - sides, respectively. The jump is invariant to the arbitrary + and - designations.
3.1. Continuous, For continuous,
piecewise-smooth
functions
piecewise-smooth
Y” = ( u” 1U”E Co(C),
representations
of the outer field the set of trial solutions
234”= 0 in ho, l& r ‘d-i”2 (~_iku’)
Fig. 3. Partition
in a two-dimensional
configuration
=o>
with a circular
Fig. 4. Two sides of an interior sector boundary.
is modified
(22)
interface.
I. Harari et al. I Comput. Methods Appl. Mech. Engrg. 164 (1998) 107-119
112
and Y” = 9’“. The weak form must now account for -possible jumps in derivatives across ry,,,. Thus we seek ui E Y’ and u” E 9’” such that
(23)
The Euler-Lagrange
I-1 au” dn
= 0
equations are (13)-(16)
and
on rynt
(24)
enforcing continuity of normal derivatives across interior sector boundaries. In this case du”ldn” E H-“*(T” fl r,), s = 1, . . . , n set, as required on each sector. The formulation retains its symmetry since tl w” E V” and Vu”EYO
(25) 3.2. Injinite element approximation The formulation (23) is a generalization of (17) that accounts for possible discontinuities across infinite element boundaries. We may now employ infinite element basis functions that exactly satisfy the Helmholtz equation within each sector and are continuous across infinite element boundaries. This approximation is employed directly with Eqs. (23). In most cases, the circumferential behavior of these functions on the artificial boundary differs from piecewise polynomial variation of finite element shape functions. Weak enforcement of continuity across the artificial boundary is then performed naturally within the framework of our formulation (23). Consequently, the infinite element mesh need not be the restriction of the finite element mesh to the artificial boundary. This flexibility may offer a significant advantage in practical applications. To fix ideas consider the two-dimensional configuration with a circular interface shown in Fig. 3. The outer domain 0“ is partitioned into non-overlapping sectors 0S=]R,w[X]8,,8,+,[,
s=l,...,
nsec.
(26)
Here, % < !, + , and 4, Let+ , = 8, + 2n. A typical sector is shown in Fig. 5. The infinite element approximation is based on interpolation of values on the arc r = R, 0, < 8 < 13,+, .
3.2.1. Linear infinite elements The basic two-noded infinite element is based on the representation u”=N,u”(R,~,)+N,u”(R,~,+,),
0,<0
(27)
where N, =
N2 =
Hy’(kr)
CT+, - 8
Hb”(kR)
Af?,
Hy’(kr)
0 - 0,
Hb’ ‘(kR)
A@,
(28)
(29)
are shape functions that satisfy the homogeneous Helmholtz equation within the sector, and AO,= OS+,- 0,. These functions exhibit outgoing oscillatory behavior in the radial direction and linear circumferential variation. The expressions required for infinite element arrays are simply
I. Harari et al. I Comput. Methods Appl. Mech. Engrg. I64 (1998) 107-119
Fig. 5. A typical sector in a two-dimensional
configuration
with a circular
113
interface.
(30)
N,V’C0) =
e- e,
(31)
Aey
and f-$
(R, 0) = -k
$&R,t?)=
-k
H;“‘(kR)
$+, - c9
Hb’ ‘(kR)
(32)
A@7
Hb”‘(kR)
28 - 0 H;’ ‘(kR) 4
(33)
The prime denotes differentiation with respect to the argument. The element array that is obtained from the term -(&v”/&r”, uO)~, in the second equation in (23). weakly enforcing continuity of the inner and outer fields across the interface, is proportional to the standard mass matrix of a linear line element kR
H;“‘(kR)-4 H;' “(kR)
6
2
1
[ 1
2
1
(34)
The element array that is obtained from the jump term (w”, [~3u”ldn]),~“~ in the second equation in (23), weakly enforcing derivative continuity of the outer field across the interior sector boundaries, is proportional to the standard stiffness of a linear element
z$-; The coefficient
-2
(35)
cyOOis a special case of a term that arises in higher-order
ffi HI'‘(kr)H’,“(kr) CY ylr = i R Hl”(kR)H;‘(kR)
(jr T
approximation
(36)
Integration of the Hankel functions is performed analytically, prior to computation. In subsequent implementations the unequal-order coefficient need not be integrated since it is evaluated by symmetry considerations, namely that the sum of (34) and (35) is symmetric.
I. Harari et al. I Comput. Methods Appl. Mech. Engrg. 164 (1998) 107-I 19
114
The evaluation of the integrals, when needed, is based on properties of U,(Z), solutions of the homogeneous Bessel equation of order v LJ?$u,= 0
(37)
%?,u, := z(zui>’ + (2’ - V2)U,
(38)
where
It follows from U,.s2l”U” -
u,~wup =0
(39)
that for Hankel functions, as solutions of the Bessel equation, ff v+ =p
kR
Hl”‘(kR)
p2 - v2
H;‘(kR)
HI”‘(kR) -
Hl”(kR)
(40)
1
In practice, the differentiated terms are eliminated by standard formulas for Bessel functions to simplify the implementation. Equal-order coefficients are evaluated by approximating the limit. (41)
at small values of V.Numerical experience indicates that this limit converges without difficulty for arithmetic of sufficient precision. The dependence of cyOO, the coefficient needed for the linear elements, on kR is shown in Fig. 6. The coefficient is well behaved and monotonically decreasing, indicating the diminishing importance of the jump terms for a given problem, as the location of the artificial boundary is moved outwards. As previously noted, the infinite element mesh need not match the finite element mesh on the interface r,. If linear infinite elements coincide with linear finite elements then the coupling terms in (23) are simply
real imaginary
Fig. 6. The coefficient
required
______________
for two-noded
infinite elements.
I. Harari et al. I Comput. Methods Appl. Mech. Engrg. 164 (1998) 107-119
proportional to (34). Otherwise, more elaborate flexibility offered by the formulation.
treatment
of the coupling
terms is required
115
for the added
3.2.2. High-order injinite elements Several extensions to higher-order approximation are possible. The basic interpolation in (27) can be enriched by a bubble Hx’(kr) sin(v,(8 - 19,)) where V, = nnlAt9,. The added coefficient in the approximation is a generalized degree of freedom rather than a nodal value. Regarding radial integration, only CZ,,,,is required as in the case of the linear element. Since 4, “I is associated with a bubble function in the circumferential direction it need not be evaluated. The matrix equations are statically condensed back to the size of the original two-noded element so that no additional computational effort is required in the solution. Alternatively, the outer field may be based on a truncation of the series representation
u” =
2
,?=0
In particular,
Hli,‘(k r )( a, COS(V,(~ - k),)) + b, sin(v,,(o - 0,)))
a three-degree-of-freedom H;“(kr)
1 N’ = 5
(
H;“(kR)
combines
the interpolatory
shape functions
HI,: ‘(kr) +
Hb’‘(kr) Hb”(kR)
element
(42)
Hj,; ‘(kR)
COO, (0 - es ))
(43)
cos(v, (0 - e,>>
(44)
HL: ‘(kr)
- H1,:‘(kR)
with the bubble Hy '(kr) sin(v, (0 - 19,))~again associated with a generalized degree of freedom. Due to the orthogonality of these functions, terms multiplying equal-order radial integrals aPI* vanish, so that no radial integration whatsoever is required. Again, static condensation of internal degrees of freedom is employed.
4. Numerical results: Radiation from a sector of a cylinder Numerical results are presented for exterior problems bounded internally by infinite circular cylinders of radius a. Soft boundary conditions are specified on the wet surface to represent a pressure-release cylinder. We consider the non-uniform radiation from an infinite circular cylinder with a constant inhomogeneous value on an arc (-a < ~9< cu) and vanishing elsewhere, so that there are two points of discontinuity in the boundary data. The normalized analytical solution to this problem for a cylinder of radius a is
(45) The prime on the sum indicates that the first term is halved. For low wave numbers this solution is relatively uniform in the circumferential direction. The directionality of the solution grows as the wave number is increased, and the solution becomes attenuated at the side of the cylinder opposite the radiating element. A circular interface r, is located at R = 2.5~ Finite elements are used in the inner domain LJi and infinite elements in the outer domain no. In the following examples the mesh of infinite elements (linear unless otherwise noted) coincides with the finite element mesh on the interface (although this is not required by the formulation). The inner domain ni is discretized by 3 X 36 bilinear quadrilateral finite elements (Fig. 7). We examine the problem with a geometrically non-dimensionalized wave number ka = 1 (the wavelength is about three times the diameter of the cylinder and four times the width of the domain), and a resolution of approximately 13 nodal points per wavelength in the finite element mesh.
4.1. Linear elements We select (Y= 20”. Fig. 8 shows the real part of the analytical and numerical solutions. The low-amplitude oscillations of the analytical solution in the vicinity of the wet surface are merely an artifact of the truncated
I. Harari et al. I Comput. Methods Appl. Mech. Engrg. 164 (1998) 107-I 19
116
Fig. 7. The computational domain quadrilateral finite elements.
exterior
to a cylinder
of radius a, with a circular
interface
at R = 2Sa, discretized
by 3 X 36 linear
series representation (60 terms) of the discontinuity in the boundary data, and are not relevant to the validation of the numerical results. The numerical solution captures the essential physics of the problem, without visible reflection from the interface. The error relative to the interpolation of the analytic solution, measured in the energy norm, is
(46) Fig. 9 depicts the solutions along the interface. The performance of the infinite elements is again quite good.
Fig. 8. Radiation denote numerical
from a sector of a cylinder, solution).
Fig. 9. Radiation
ka = 1, - 13 points per wavelength
from a sector of a cylinder,
(dashed contours denote analytical
ka = 1, -13 points per wavelength,
along the interface
solution, solid contours R = 2.5a.
I. Harari et al. I Comput. Methods Appl. Mech. Engrg. 164 (1998) 107- I19
Fig. 10. Radiation from a sector of a cylinder, ka = 1, - 13 points per wavelength, analytical solution, solid contours denote numerical solution). Fig.
bubble-enriched
11. Radiation from a sector of a cylinder, ka = 1, -13 points per wavelength,
117
infinite elements (dashed contours denote
bubble-enriched
infinite elements,
along the interface
R = 2Sa.
4.2. Bubble-enriched
elements
Here, elements enriched by the sine bubble are employed. enriched numerical solutions. The relative error
bh- 4, =
1o
9o/ .
,
,,
Fig. 10 shows the real part of the analytical
and
(47)
0
Iu’ I I
‘.._ :., ‘.., ‘... ‘\
j
;: ,.:
/’
,:’
j’ ,I’ :
t
,%.,
Fig. 12. Radiation from a sector of a cylinder, ka = 1, - 13 points per wavelength, analytical solution, solid contours denote numerical solution). Fig. 13. Radiation R = 2Sa.
from a sector of a cylinder,
ka = 1, -13
trigonometric
points per wavelength,
infinite elements (dashed contours
trigonometric
infinite elements,
denote
along the interface
118
I. Harari et al. I Comput. Methods Appl. Mech. Engrg. 164 (1998) 107-119
shows improvement over the linear infinite element. Fig. 11 depicts the solutions along the interface. The improvement over the linear infinite elements is apparent. 4.3. Trigonometric-based
elements
Here, three-degree-of-freedom elements based on a truncation of the trigonometric series (42) are employed. Fig. 12 shows the real part of the analytical and numerical solutions. The relative error (48) is a degradation from the linear infinite element. Fig. 13 depicts the solutions along the interface.
5. Conclusions
In this work, higher-order approximations are implemented in a novel method for developing infinite elements for exterior problems of time-harmonic acoustics. This technology is based on a simplified functional (9) for the partitioned problem, weakly enforcing continuity at the interface without additional degrees of freedom. The infinite elements mesh only the outer boundary of the finite element domain and need not match the finite elements on the interface. The implementation of these infinite elements is straightforward. Radial integration, when required, is performed analytically. Two approaches to higher-order approximation are presented: one is based on enrichment of low-order elements, and the other on the truncation of a trigonometric series. The second approach is free of radial integration. Infinite elements with three degrees of freedom developed by both approaches are examined. In both cases, static condensation is employed so that the computational effort required in the solution of the higher-order elements is no greater than for lower-order elements. Numerical results validate the improved performance of the bubble-enriched elements and indicate a need to employ more terms in the trigonometric series-based elements.
Acknowledgment
This research was supported by the U.S. Office of Naval Research. The authors wish to thank Parama Barai for calculation of error norms presented herein.
References [l] N.N. Abboud and P.M. Pinsky, Finite element dispersion analysis for the three-dimensional second-order scalar wave equation, Int. J. Numer. Methods Engrg. 35(6) (1992) 1183-1218. [2] 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. [3] R.J. Astley, G.J. Macaulay and J.-P. Coyette, Mapped wave envelope elements for acoustical radiation and scattering, J. Sound Vib. 170(l) (1994) 97-118. [4] I. Babulka, F. Ihlenburg, E.T. Paik and S.A. Sauter, A generalized finite element method for solving the Helmholtz equation in two dimensions with minimal pollution, Comput. Methods Appl. Mech. Engrg. 128(3-4) (1995) 325-359. [5] P. Bettess, Infinite elements, Int. J. Numer. Methods Engrg. 11(l) (1977) 53-64. [6] D.S. Burnett, A three-dimensional acoustic infinite element based on a prolate spheroidal multipole expansion, J. Acoust. Sot. Am. 96(5) (1994) 2798-2816. [7] H.S. Chen and CC. Mei, Oscillations and wave forces in a man-made harbor in the open sea, in: Tenth Symposium on Naval Hydrodynamics, pp. 573-594, Cambridge, Massachusetts, 1974. Office of Naval Research, Department of the Navy. [8] K. Feng, Asymptotic radiation conditions for reduced wave equation, J. Comput. Math. 2(2) (1984) 130-138. [9] R.H. Gallagher, E. Haugeneder and H.A. Mang, Limitations on methods of handling incompatible basis functions in boundary value problems, Trans. Am. Nucl. Sot. 33 (1979) 326-327. [IO] 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) (1996) 239-273.
1. Harari et al. I Comput. Methods Appl. Mech. Engrg. 164 (1998) 107-119 [l I] D. Givoli, Numerical methods for problems Publishing Co., Amsterdam, 1992).
in infinite domains,
in: Vol. 33 Studies
in Applied
Mechanics
119 (Elsevier
Scientific
[l2] D. Givoli and J.B. Keller, A finite element method for large domains, Comput. Methods Appl. Mech. Engrg. 76( 1) (1989) 41-66. [I31 K. Grosh and PM. Pinsky, Complex wave-number dispersion analysis of Galerkin and Galerkin least squares methods for fluid-loaded plates, Comput. Methods Appl. Mech. Engrg. I 13( 1-2) (1994) 67-98. [I41 I. Harari, P.E. Barbone and J.M. Montgomery Finite element formulations for exterior problems: Application to hybrid methods, non-reflecting boundary conditions and infinite elements, Int. J. Numer. Methods Engrg. 40( 15) (1997) 2791-2805. [15] I. Harari, P.E. Barbone, M. Slavutin and R. Shalom, Boundary infinite elements for the Helmholtz equation in exterior domains, Int. J. Numer. Methods Engrg. 41(6) (1998) 1105-I 131. [I61 I. Harari, methods [ 171 I. Harari domains,
K. Grosh, T.J.R. Hughes, M. Malhotra, PM. Pinsky, J.R. Stewart and L.L. Thompson, Recent developments in finite element for structural acoustics, Arch. Comput. Methods Engrg. 3(2-3) (1996) 131-309. and T.J.R. Hughes, Analysis of continuous formulations underlying the computation of time-harmonic acoustics in exterior Comput. Methods Appl. Mech. Engrg. 97(l) (1992) 103-124.
[ 181 I. Harari and T.J.R. Hughes, Galerkinlleast-squares [19] [ZO] [2l] [22] [23] [24]
finite element methods for the reduced wave equation with nonreflecting boundary conditions in unbounded domains, Comput. Methods Appl. Mech. Engrg. 98(3) (1992) 41 l-454. 1. Harari and T.J.R. Hughes, Studies of domain-based formulations for computing exterior problems of acoustics, Int. J. Numer. Methods Engrg. 37( 17) (1994) 2935-2950. E. Haugeneder and H.A. Mang, On an improper modification of a variational principle for finite element plate analysis, Z. Angew. Math. Mech. 59( 1I) (1979) 637-640. E. Haugeneder and H.A. Mang, Admissible and inadmissible simplifications of variational methods in finite element anal,ysis, in: S. Nemat-Nasser, ed., Variational Methods in the Mechanics of Solids, (Pergamon Press, Oxford, 1980) 187-198. S.W. Key, A specialization of Jones’ generalization of the direct-stiffness method of structural analysis, AIAA J. 9(5) (197111 984-985. H.A. Mang and R.H. Gallagher, A critical assessment of the simplified hybrid displacement method, Int. J. Numer. Methods Engrg. 1l(1) (1977) 145-167. M. Masmoudi, Numerical solution for exterior problems, Numer. Math. 51(l) (1987) 87-101.
[25] A. Oberai and P.M. Pinsky, A multiscale finite element method for the Helmholtz equation, Comput. Methods Appl. Mech. Engrg., 154(3-4) (1998) 281-297. [25] RI? Shaw, Integral equation methods in acoustics, in: C.A. Brebbia, ed., Boundary Elements X, Vol. 4 (Computational Mechanics Publications, Southampton, 1988) 22 I-244. [27] I. Stakgold, Boundary Value Problems of Mathematical Physics, Vol. II (Macmillan, New York, 1968). [28] J.R. Stewart and T.J.R. Hughes, h-adaptive finite element computation of time-harmonic exterior acoustics problems in two dimensions, Comput. Methods Appl. Mech. Engrg. 146( l-2) (1997) 65-89. [29] L.L. Thompson and P.M. Pinsky, A Galerkin least-squares finite element method for the two-dimensional Helmholtz equation, Int. J. Numer. Methods Engrg. 38(3) (1995) 371-397. [30] C.H. Wilcox, Scattering Theory for the d’Alembert Equation in Exterior Domains (Springer-Verlag. Berlin, 1975). [31] 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(7) (1985) 1229-1251.