3-D Acoustic Infinite Element Based on a Prolate Spheroidal Multipole Expansion

In the field of computational structural acoustics, the problem of efficiently modeling the acoustic field in large exterior domains has remained a difficult challenge for over a quarter century. Many techniques have evolved. Those historically receiving greatest attention in the literature are: the boundary integral equation method or, as it is often called, the boundary element method (BEM), which is based on the surface Helmholtz integral equation applied to the surface of the structure; and infinite elements, which are based on the Helmholtz differential equation (the reduced wave equation), applied to semi-infinite sectors of the domain that are exterior to an artificial boundary surrounding the structure, with finite elements between the structure and the boundary.

The BEM is the method of choice for most researchers and code developers. Its primary advantage is that it is based on a mathematically rigorous formulation, namely, the Helmholtz integral representation, that (i) satisfies the Sommerfield radiation condition and (ii) represents the exact solution throughout the exterior domain. In addition, its reduction of the dimensionality of the problem by one has long been thought to be a significant computational advantage. However, we have found that this is actually disadvantageous because of the much greater bandwidth of the equations, due to the inherent global connectivity of the method.

By contrast, infinite elements were hitherto never intended to represent the exact solution, or even an approximate solution, within the element itself. Based on ad hoc physical assumptions or asymptotic theories, they were only meant to provide an approximate non-reflecting condition on an artificial boundary surrounding the structure, thereby enabling a sufficiently accurate solution to be obtained in the finite region interior to the boundary. The solution exterior to the boundary is then usually obtained from the exterior Helmholtz integral, using the pressure and velocity data computed on the surface of the structure.

Two types of infinite elements have previously been used for acoustic applications: the exponential decay and the "mapped" element.

The exponential decay element approximates the spatial decay of the acoustic pressure amplitude by an exponential, p.varies.e.sup.-.gamma.r e.sup.-ikr, where .gamma. is an empirically adjusted positive number and r is any coordinate that extends to infinity. Because this decay function is inconsistent with acoustical theory, the accuracy rapidly deteriorates away from the artificial boundary.

The mapped element is based on the asymptotic form of the field for large r, using the leading terms of the lower order spherical Hankel functions, namely, ##EQU1## Because of the mapping, the element was intended to be usable in any shape and orientation, i.e., not necessarily aligned with spherical coordinate surfaces, thereby permitting the "radial" sides of different elements to emanate from different spherical origins.

However, we have found that for this representation to converge, (i) the elements must lie outside the smallest sphere circumscribing the structure, and (ii) the sides of all elements in a mesh must conform to radial lines emanating from the origin of the single r coordinate in the expansion. Since the fluid region between the circumscribing sphere and the structure must be meshed with finite-size acoustic elements, the total acoustic mesh would become very large and therefore computationally very inefficient for structures of large aspect ratio. Thus, the element is practical only for "chunky" structures, i.e., those that can be closely circumscribed by a sphere. This is a serious limitation for many applications.

An additional problem of the mapped element is that element generation is relatively expensive. Moreover, generation of the matrices for the mapped element requires inversion of an ill-conditioned matrix. This causes quadrature problems at very high frequencies.

The spheroidal infinite element described in this paper breaks with tradition in that it is based on a mathematically rigorous formulation, namely, a new multipole expansion (a power series in 1/r, where r is a prolate spheroidal radius) that (i) satisfies the Sommerfeld radiation condition and (ii) represents the exact solution throughout the domain exterior to a minimal circumscribing prolate spheroid. When the leading terms of the expansion, a polynomial in 1/r multiplied by an oscillatory factor, are incorporated in an infinite element, accuracy may be increased arbitrarily throughout the element by increasing the number of terms in the polynomial, analogous to p-extension in conventional polynomial-based finite elements.

The spheroidal element improves on the mapped element by removing all of the aforementioned practical limitations. It (i) eliminates the 3-D mapping that mixes together the (troublesome) infinite direction with the (well-behaved) finite directions, (ii) generalizes the circumscribing sphere to a circumscribing spheroid in order to accommodate structures of both large and small aspect ratios and (iii) takes advantage of separable coordinates to reduce the 3-D integrals to products of much more cheaply computed 2-D and 1-D integrals, which simultaneously eliminates the ill conditioning.

We believe that the inventive technique is probably the most efficient technique now available for large computational structural acoustics problems. In addition, the technique used for this acoustic element can be applied to other types of wave propagation such as electromagnetism and elastodynamics.

Our new element is readily incorporated in any structural finite element code. We have, in fact, implemented the inventive technique and incorporated this new element in such a code, which is referred to below as our exemplary IEM code. We have tested this code on a variety of problems to validate its accuracy, as discussed below.

Finite element codes, in general, are well known and need not be described here in detail. Numerous books are available on this subject, including D. S. Burnett, Finite Element Analysis, Addison-Wesley, Reading, Mass., 1987 (reprinted 1988), which is hereby incorporated by reference.

FIG. 1 is a flowchart illustrating the basic steps in the use and operation of a finite element code.

FIG. 2 depicts an arbitrary, prolate structure together with its smallest circumscribing sphere S, having radius r.sub.0.

FIGS. 3a-3d depict the closed coordinate surfaces of four orthogonal coordinate systems. These surfaces are: (a) sphere; (b) prolate spheroid; (c) oblate spheroid; and (d) ellipsoid.

FIG. 4 depicts a set of confocal ellipses and hyperbolas. When rotated about the z-axis, these curves form the spheroidal and hyperboloidal surfaces, respectively, of the prolate spheroidal coordinate system.

FIG. 5 illustrates the prolate spheroidal coordinates r, .theta., and .phi., and related geometric parameters. In the figure, the prolate radius r is defined as ##EQU2##

FIG. 6 depicts an arbitrary, prolate structure and a minimal circumscribing prolate spheroid S, with prolate radius r.sub.0.

FIG. 7 depicts an illustrative prolate spheroidal infinite element.

FIG. 8 depicts a typical mesh using infinite elements according to the invention.

FIG. 9 shows the base of an illustrative infinite element, lying on the infinite element spheroid. The dashed lines in the figure are .theta., .phi.-coordinate lines.

FIG. 10 is a 2D slice through a 3D mesh for a typical structural acoustics problem.

FIG. 11 shows the geometric discontinuity at an interface between adjacent infinite and finite-size elements.

FIG. 12 illustrates a local node numbering convention for the particular case of a quadratic Lagrange quadrilateral (n=9) in the angular directions and a quadrupole (m=3) in the radial direction (N=m.times.n=27).

FIG. 13 compares the performance of a BEM code to that of the exemplary IEM code for four different meshes for a structural acoustic problem. The IEM fluid mesh, which is topologically identical in the radial direction for all four meshes, was adjusted to yield approximately the same rms error as the BEM on the surface of the structure.

FIG. 14 is a graph of the normal velocity along the meridian of a cylindrical thin shell with rigid, flat, end caps. The cylinder has an aspect ratio of 6.7:1, and is excited by a ring force. The product ka is equal to 1. The solid curve represents the analytical solution, and the discrete points represent the infinite element calculation.

FIG. 15 is a graph of the far-field radiation pattern of a cylindrical thin shell having elastic, hemispherical end caps. The cylinder has an aspect ratio of 6:1, and is excited by an axial point force. The product ka is equal to 1.9. As in FIG. 14, the solid curve represents the analytical solution, and the discrete points represent the infinite element calculation.

The finite element method (FEM) is a method, performed with the help of a computer, for predicting the behavior of a physical system by obtaining numerical solutions to mathematical equations that describe the system and its loading conditions. The use of the FEM may be thought of as comprising three phases: preprocessing, solution, and postprocessing. These phases are now discussed in further detail with reference to FIG. 1.

In the preprocessing phase, the physical domain of the problem is partitioned into a pattern of subdomains of simple geometry, referred to as "elements". The resulting pattern is referred to as a "mesh". In addition, problem data such as physical properties, loads, and boundary conditions are specified. This procedure is represented as step 10 of the figure.

The solution phase comprises steps 20-50 in the figure. In the first of these four steps, numerical values are computed for the coefficients in the "element matrix equation" associated with each element in the mesh. The element matrix equation, the derivation of which is indicated in box 15 of the figure, is a set of numerically computable mathematical formulas that are derived theoretically and implemented into the computer code that forms the FEM program.

During the use of the FEM program, the code for these formulas is accessed by each element in the mesh, and numerical values are computed for the coefficients in the formulas using the geometric and physical data associated with each element.

We are providing, for the first time, a set of such formulas (i.e., the element matrix equation) that is specific to a prolate spheriodal acoustic infinite element.

The procedure used in the derivation of the element matrix equation, which is described in detail below, embodies the following general ideas. The unknown field variable, for which the finite-element analysis is seeking a solution (such as the acoustic pressure), is represented approximately within each element as a finite sum of known functions, referred to as "shape" functions. These shape functions are usually chosen to be polynomials. There are unknown parameters in this representation, referred to as "degrees of freedom (DOF)", that become the new unknowns which the finite-element analysis will find values for. The DOF are often the values that the unknown field variable takes at specific points in the element, referred to as "nodes". The purpose of this representation is that when values for the DOF are subsequently computed in step 50, a solution will then be known everywhere, continuously, throughout each element. This is possible because at that stage, both the shape functions and the parameters will be known, and these, together, define the complete solution.

The representation of the unknown field variable in terms of shape functions is then inserted into the governing physics equations (which are, typically, differential or integral equations) that express the physical laws to which the physical system is subject. These calculus equations reduce to a system of linear algebraic equations which, when written in matrix form, is the element matrix equation. The expressions for the coefficients in this matrix equation are manipulated using calculus and algebra until they are in a form that can be numerically evaluated by the computer. These numerical coefficients are the output of step 20.

In step 30, the element matrix equations from all the elements are combined together (i.e., they are said to be "assembled") into one large "system matrix equation." The matrices associated with elements that touch each other in the mesh will partially overlap, thereby establishing continuity in the field variable from element to element. Since the overlap is partial, the system matrix grows larger and larger as element equations are assembled, resulting in one very large system matrix.

The system matrix equation is then modified (step 40) to take account of the boundary and loading conditions. The system matrix equation is then solved (step 50), using conventional techniques of numerical analysis. Although there may be thousands, or even hundreds of thousands, of unknown DOF, the solution of this matrix equation is generally a very tractable problem. That is because the matrix elements tend to assume non-zero values only within a relatively narrow band along the matrix diagonal. A well-known measure of the width of this band (and thus, of the time required for solution) is the "rms half-bandwidth" B.sub.rms.

In the postprocessing phase, the solution to the system matrix equation is displayed in an appropriate, meaningful form (step 60). In this phase, other useful information may be derived from the solution and likewise displayed. The following section presents a detailed description of the derivation of the element matrix equation for the inventive prolate spheroidal acoustic infinite element.

Consider an arbitrary structure immersed in an infinite homogeneous fluid and vibrating at constant frequency .omega. (FIG. 2). Let S be a sphere of minimum radius r.sub.0 that just circumscribes the structure, as shown in FIG. 2. It is well known that the scattered and/or radiated pressure p exterior to S can be represented by the following multipole expansion in spherical coordinates r,.theta.,.phi.: ##EQU3## where k is acoustic wavenumber, the F.sub.n are smooth, differentiable functions, and the series converges in r>r.sub.0.

Since Eq. (1) is valid only outside a circumscribing sphere, then infinite elements based on Eq. (1) should only lie outside a circumscribing sphere. The region between the structure and the sphere must therefore be filled with finite-size acoustic elements. This limits the usefulness of such infinite elements to structures that are "chunky", i.e., that fill up most of the space inside a sphere. Long and/or flat structures would need an excessive number of finite-size elements (or DOF for p-type elements) to fill the space. Therefore a different multipole expansion is needed, one that is appropriate for non-chunky structures. Since multipole expansions are limited to a specific coordinate system (because they provide an explicit functional form for the "radial" coordinate that goes to infinity), the first step is to select an appropriate coordinate system.

The closed coordinate surfaces relating, respectively, to four orthogonal coordinate systems are shown in FIG. 3. They are all ellipsoids, i.e., quadric surfaces for which the three orthogonal principal cross-sections are ellipses, including circles as a degenerate form. The most general case, involving three different elliptical cross-sections, is case (d). The others are given special names: spheres, prolate spheroids and oblate spheroids. Prolate and oblate spheroids are generated by rotating ellipses about their major and minor axes, respectively.

Prolate spheroidal coordinates can closely circumscribe any structure for which two dimensions are similar and the third comparable or greater, i.e., structures that vary from "chunky" to "long". The applications of interest herein involve long (or chunky) structures, for which the appropriate coordinate system is prolate spheroidal. Therefore, described below is the development of a prolate spheroidal infinite element.

1. Prolate Spheroidal Coordinates

FIG. 4 shows a family of confocal ellipses and hyperbolas, which are orthogonal. When rotated about the major axis of the ellipses (the z-axis), the curves become quadric surfaces, called, respectively, prolate spheroids and two-sheeted hyperboloids. Combined with planes through the z-axis, these orthogonal surfaces are the coordinate surfaces comprising the prolate spheroidal system.

FIG. 5 defines the parameters relevant to the prolate spheroidal coordinates r, .theta., .phi.: ##EQU4##

Surfaces of constant r are confocal spheroids, constant-.theta. surfaces are confocal two-sheeted hyperboloids, and constant-.phi. surfaces are planes containing the z-axis.

The transformation from Cartesian to prolate spheroidal coordinates is as follows: ##EQU6##

2. A Prolate Spheroidal Multipole Expansion

Consider the same structure as in FIG. 2, but now let S be a prolate spheroid of minimal prolate radius r.sub.0 that just circumscribes the structure (FIG. 6). The scattered and/or radiated pressure p exterior to S can be represented by the following multipole expansion in prolate spheroidal coordinates r, .theta., .phi.: ##EQU7## with properties similar to Eq. (1). The series in Eq. (6) converges absolutely and uniformly in r,.theta. and .phi. in any region r.gtoreq.r.sub.0 +.epsilon.>r.sub.0. There is a six term recursion formula that defines the functions G.sub.n (.theta.,.phi.;k) in terms of G.sub.0 (.theta.,.phi.;k).

Eq. 6 is a generalization of Eq. (1), including the latter as a special case. Whereas Eq. (5) is appropriate for a wide range of shapes, varying from chunky to arbitrarily long, Eq. (1) is only appropriate for chunky shapes. Stated geometrically, prolate spheroids have a wide range of shapes but spheres have only one.

The appropriateness of Eq. (6) for long structures can be seen by considering the difference in radial response along the major and minor axes of the prolate spheroidal coordinates. Along the major axis the prolate radial coordinate, r, is the distance from the origin, a; hence, r is the same as the spherical radius. In any other direction r is not the distance from the origin. In particular, along the minor axis b is the distance from the origin. Since r=.sqroot.b.sup.2 +f.sup.2 , then close to the structure, where b/f<<1, dr/db.congruent.b/f. Thus, r is changing more slowly than b. This can be appreciated by reexamining FIG. 4; since spheroids are surfaces of constant r, a 1/r.sup.n decay is spread over greater radial distances at .theta.=.pi./2 than at .theta.=0. In short, the spatial rate of decay in the radial direction in Eq. (6) is like spherical spreading along the major axis, but close to the side of a long structure the decay is slower, similar to cylindrical (2-D) spreading, as one might expect. In the far field, a, b, and r approach equality, so the decay approaches spherical spreading in all directions.

The infinite element is shown in FIG. 7, and a 2-D cross-section of a typical 3-D mesh outside an irregular, 3-D body is shown in FIG. 8. One face of the element, the "base", must attach to, and conform to the shape of, a prolate spheroidal surface of prolate radius r.sub.1 surrounding the structure, called the "infinite element spheroid". The base may be a curvilinear quadrilateral (shown here) or triangle; it need not conform to the .theta.,.phi. coordinate lines (FIG. 9). This permits one to construct a completely general 2-D mesh on the spheroid, comprising curvilinear quadilaterals and triangles of any shape and orientation. Any type of FE representation may be used over the base, e.g., Lagrange, serendipity or hierarchic polynomials of any degree. (The quadratic Lagrange nodal pattern shown here, and in later figures, is for illustrative purposes.)

The infinite element spheroid uniquely determines two loci, which in turn define a system of confocal spheroids and hyperboloids (cf. FIG. 4) that are used in the construction of the rest of the element. The side faces of the element are the loci of confocal hyperbolas emanating from the sides of the base. A multipole of order ##EQU8## requires m layers of nodes that are on confocal spheroids of prolate radii r.sub.1, r.sub.2, . . . , r.sub.m. The nodal pattern is identical on each spheroid. The value m=2 corresponds to a dipole, m=3 to a quadrupole (as shown in FIG. 7), m=4 to an octupole, etc. A "variable-multipole" element contains a variable number of layers of nodes, permitting the analyst to refine the mesh radially in the infinite region by increasing the multipole order, m, analogous to p-extension for finite-size elements.

Finally, there is an outer spheroidal face, S, at prolate radius r, that recedes to infinity in the theoretical development. Thus, in deriving the element matrix equation, the element begins with a finite size, r.sub.1 .ltoreq.r.ltoreq.r, so that the Sommerfeld radiation condition can be applied to the outer face, and then the limit is taken as r.fwdarw..infin..

Conventional .xi.,.eta. coordinates (illustrated in FIG. 9 for a quadrilateral) are defined over the element cross-section, i.e., on spheroids, by the mapping ##EQU9## where n is the number of nodes in the quadrilateral or triangle, .theta..sub..nu.,.phi..sub..nu. are the prolate spheroidal angular coordinates of the .nu..sup.th node, and .chi..sub..nu..sup.a (.xi.,.eta.) are interpolation polynomials. (Alternatively, blending functions could be used if the elements are hierarchic.)

Since the base, intermediate nodal layers, and outer face conform to confocal spheroids and the side faces are the loci of confocal hyperbolas, the element is a right cylinder in r,.theta.,.phi.-space (or r,.xi.,.eta.-space) and the 3-D integrals separate into well-conditioned 2-D angular integrals and 1-D infinite integrals. An additional benefit is that the above mapping is limited to just the two finite angular dimensions, i.e., the .theta.,.phi. coordinates, over spheroidal surfaces. The radial coordinate uses a separate mapping. This avoids the numerical ill-conditioning associated with 3-D mappings that mix together the infinite dimension with the two finite dimensions.

The use of these infinite elements is illustrated in FIG. 10, which shows a typical mesh for a structural acoustic problem. It is a cross-section of a 3-D mesh of the structure and the surrounding infinite fluid region. The structure is a circular cylindrical shell with rounded conical end caps and containing some internal structures. There are two layers of finite-size acoustic elements between the structure and the infinite element spheroid and then a layer of infinite elements, shown with dashed lines, outside the spheroid. The dashed lines are all hyperbolas. The first layer of finite-size elements is generated by projecting all the nodes on the shell's outer surface the same distance outward along normals to the shell, creating a surface parallel to the shell. The second layer projects the nodes on the parallel surface outward, along normals, to the spheroid. This is a simple procedure for automatic mesh generation for convex surfaces; concave surfaces would require a more sophisticated 3-D mesh generator.

There may or may not be a discontinuity (in .xi.,.eta. coordinates and therefore in the dependent variable, pressure) between the infinite elements and adjacent finite-size elements along the infinite element spheroid (FIG. 11), depending on how the geometry for the finite-size elements is generated. If these finite-size elements use a conventional polynomial mapping based on the global Cartesian coordinates of the nodes, then the elements will not conform exactly to the shape of the spheroid. However, the mapping in Eq. (7) for the infinite element defines a set of .xi.,.eta. coordinates that lie exactly on the spheroid. Thus there is a "sliver" of space, comprising gaps and overlaps, between the finite and infinite elements, and the .xi.,.eta. coordinates in the two elements do not exactly coincide on their interface. This is not a problem because the error due to this geometric discontinuity converges to zero, as the mesh is h or p-refined, at the same rate as does the (well-known) error at the boundary of a curved domain when using conventional polynomial-mapped elements. Indeed, simple numerical calculations show that the maximum spatial separation between identical .xi.,.eta. values on the faces of adjacent finite and infinite elements is typically several orders of magnitude less than internodal distances. (Using blending functions for the geometry of the finite-size elements would, of course, eliminate the discontinuity.)

Time-harmonic (e.sup.i.omega.t) behavior is governed by the 3-D Helmholtz equation,

where k is the wavenumber (=.omega./c), c is sound speed (=.sqroot.B/.rho.), B is bulk modulus, .rho. is density, and .rho. is the complex-valued amplitude of scattered and/or radiated pressure: ##EQU10##

To ensure uniqueness of the solution, the pressure must satisfy the Sommerfeld radiation condition at the outer "boundary" at infinity: ##EQU11## where the lower case o, read "little o", means "faster than". The r in Eq. (10) is a spherical r. However, prolate spheroidal coordinates approach spherical coordinates as r.fwdarw..infin., so Eq. (10) can be used with prolate spheroidal coordinates in the limit as r.fwdarw..infin..

When Eq. (10) is applied to the outer face of the element S, at prolate radius r, as it recedes to infinity, Eq. (10) can be satisfied by the proper choice of shape functions.

1. General Expression; DOF Numbering

The scattered and/or radiated pressure is represented as follows, ##EQU12## Here .psi..sub..nu..sup.a (.xi.,.eta.) are "angular" shape functions that interpolate p over spheroidal surfaces confocal to the infinite element surface, and .psi..sub..mu..sup.r (r) are "radial" shape functions that interpolate p along hyperbolic rays. Interelement C.sup.0 -continuity is established by enforcing the interpolation property:

The currently preferred local node numbering convention is to begin with the nodes on the base of the element and then proceed radially outward, a layer of nodes at a time. This is summarized in Table I and illustrated in FIG. 12.

2. Angular Shape Functions

The functions .psi..sub..nu..sup.a (.xi.,.eta.) are conventional 2-D polynomials (serendipity, Lagrange or hierarchic).

For example, for the quadratic Lagrange quadrilateral elements depicted in FIGS. 7, 9 and 12, the angular shape functions are ##EQU13##

If the functions .psi..sub..nu..sup.a (.xi.,.eta.) are also used for the mapping functions .chi..sub..nu..sup.a (.xi.,.eta.) in Eq. (7), then the element is "isoparametric" in the angular directions.

3. Radial Shape Functions

The functions .psi..sub..mu..sup.r (r) use a truncated form of the radial part of the multipole expansion in Eq. (6), namely, an mth order multipole expansion: ##EQU14## The phase factor e.sup.ikr.sbsp..mu. does not need to be included; if omitted, it will simply appear in each h.sub..mu..mu.' in Eq. (20) below. The factors k.sup..mu.' in the denominators are also not necessary; they are included only to make the h.sub..mu..mu.' nondimensional.

The coefficients h.mu..mu.' are determined by the requirement of interelement C.sup.0 -continuity. Applying Eq. (14) to Eq. (16) yields m sets of m linear algebraic equations:

where

and [I] is the identity matrix. Therefore,

This procedure defines m layers of nodes (with n nodes on each layer) lying on spheroids of prolate radii r.sub.1,r.sub.2, . . . ,r.sub.m (cf. FIG. 12).

To illustrate, consider a dipole element (m=2), ##EQU15## Inverting a 2.times.2 [S] matrix yields ##EQU16##

The procedure in Eqs. (16)-(20) is the one that has been used to date. However, a hierarchic formulation would have the usual advantages of ease of mesh refinement (by p-extension), improved numerical conditioning and elimination of all nodes exterior to the infinite element spheroid. The last advantage is especially important because it would eliminate the need to locate nodes on hyperbolic trajectories, which is not a standard technique in FE mesh generators. To convert to a hierarchic formulation, the angular directions would employ the standard 2-D hierarchic shape functions for quadrilaterals and triangles. The radial direction would use the mapping .zeta.=1-2r.sub.1 /r, which linearly maps the interval 1/r.epsilon.[1/r.sub.1,0) to the interval .zeta..epsilon.[-1,1), and then employ the standard 1-D hierarchic shape functions in .zeta., excluding the linear function that is unity at infinity, i.e., at .zeta.=1.

4. Stiffness, Mass, and Radiation Damping Matrices, and the Element Matrix Equation

Further mathematical details of the derivation of these expressions may be found in the Appendix attached hereinbelow as well as in the article by D. S. Burnett, "A 3-D Acoustic Infinite Element Based on a Prolate Spheroidal Multipole Expansion," J. Acoust. Soc. Am. 96 (4), October, 1994, which is hereby incorporated by reference.

1. Formal Expressions

The matrices for the element matrix equation are derived by starting with a finite-size element, i.e., with the outer face on a spheroid of prolate radius r (cf. FIG. 7), and then taking the limit as r.fwdarw..infin.. Thus, applying the Galerkin-weighted residual method to Eq. (8) over a single element yields ##EQU17## using k.sup.2 =.omega..sup.2 .rho./B. The first integral is converted into a surface integral and another volume integral using the identity .psi..sub.i .gradient..sup.2 p=.gradient..(.psi..sub.1 .gradient.p)-.gradient..psi..sub.i..gradient.p and the divergence theorem. Substituting Eq. (11) into the volume integrals yields the following element matrix equation:

where the stiffness matrix [K], mass matrix [M], and pressure gradient vector {F} are, respectively, ##EQU18##

The surface integral for {F} is over the entire boundary of the element, S.sup.(e). It is split into two integrals: one over the outer face S.sup.(e) and the other over the remaining faces S.sup.(e) -S.sup.(e), ##EQU19##

Consider the first integral in Eq. (26). As r.fwdarw..infin., prolate spheroidal coordinates approach spherical coordinates so .differential.p/.differential.n.fwdarw..differential.p/.differential.r and dS.fwdarw.r.sup.2 sin .theta.d.theta.d.phi., where r, .theta.,.phi. are spherical coordinates. To evaluate .differential.p/.differential.r as r.fwdarw..infin., substitute Eqs. (11), (12), and (16) into Eq. (10), which yields a stronger form of the Sommerfeld condition, ##EQU20## where the upper-case O, read "big oh," means "at least as fast as." Substitute .differential.p/.differential.r in Eq. (27) for .differential.p/.differential.n in the first integral. The O(1/r.sup.2) term makes no contribution to the integral because .psi..sub.i is O(1/r) and dS.varies.r.sup.2 as r.fwdarw..infin.. In the remaining term ikp, substitute Eq. (11) for p.

In the second integral in Eq. (26), substitute the balance of linear momentum for .differential.p/.differential.n, viz., .differential.p/.differential.n=.omega..sup.2 .rho.u.sub.n, where u.sub.n is the amplitude of the normal component of particle displacement.

Making these substitutions, Eq. (26) becomes

where ##EQU21##

The [C] matrix is the radiation damping matrix, representing radiation loss to infinity. The {D} vector permits specification of u.sub.n on the side or bottom faces of the element.

Substituting Eq. (28) into Eq. (24) yields, for the element matrix equation,

where, summarizing, the formal expressions for the stiffness, mass, and radiation damping matrices are, respectively, ##EQU22## The D.sub.i are zero in virtually all practical applications.

2. Transformation of Integrals; Final Expressions

The remaining mathematics transform the integrals in Eq. (31) to expressions that can be numerically evaluated. (The D.sub.i integrals are ignored because they are zero in virtually all practical applications.) Mathematical details may be found in the Appendix; following is a brief description of the principal steps.

Transform the Integrals in Eq. (31) to Prolate Spheroidal Coordinates r,.theta.,.phi..

The differential volume and surface elements are dV=J.sub..nu. drd.theta.d.phi. and dS=J.sub.s d.theta.d.phi., where J.sub..nu. is the volume Jacobian and J.sub.s is the surface Jacobian. J.sub..nu. and J.sub.s, as well as the gradient operator, V, in K.sub.ij in Eq. (31), are expressed in prolate spheroidal coordinates. This separates the 3-D integrals for K.sub.ij and M.sub.ij into products of 2-D "angular" integrals over .theta.,.phi. and 1-D "radial" integrals over r, and the 2-D integral for C.sub.ij into the product of a 2-D angular integral and a function of r.

Develop Final Expressions for the Angular Integrals

Transform the .theta.,.phi. coordinates to local .xi.,.eta. coordinates using the coordinate mapping in Eq. (7) [cf. FIG. 9]. The resulting well-defined integrals can be numerically integrated using Gauss rules in the conventional FE manner.

Develop Final Expressions for the Radial Integrals for K.sub.ij and M.sub.ij and Radial Function for C.sub.ij

Substitute Eq. (16) into each of the radial integrals and radial function and perform various algebraic operations. Some of the integrals become well-defined Fourier sine or cosine transforms, which can be evaluated by standard algorithms available in many mathematics software packages. The other integrals, as well as the radial function, result in undefined oscillatory terms, which are treated in the next step.

Form Final Expression for Element Matrix Equation

All the above expressions, including both the well-defined integrals and the undefined oscillatory terms, are substituted into Eq. (30), the element matrix equation. The parenthetical expression for the three matrices becomes, ##EQU23## where [K.sup..infin. ] and [M.sup..infin. ] comprise all the well-defined integrals. The undefined oscillatory terms all cancel, leaving as the final form of the element matrix equation,

where the expressions for K.sub.ij.sup..infin. and M.sub.ij.sup..infin. are as follows (D.sub.i are zero in virtually all practical applications): ##EQU24## where B is bulk modulus, .rho. is density, r.sub.1 is the prolate radius of the infinite element spheroid, and .epsilon..sub.1 (=f/r.sub.1) is the eccentricity of the infinite element spheroid.

The angular integrals, A.sub..nu.'.nu..sup.(i), i=1, . . . ,4, .nu.',.nu.=1, . . . ,n (cf. Table I), are ##EQU25## All four integrals in Eq. (35) can be numerically integrated using standard Gauss rules since the integrands are smooth and bounded (including the 1/sin .theta. term because J.infin..theta. as .theta..fwdarw.0).

The Jacobian, J, is computed from the coordinate mapping in Eq. (7), ##EQU26## and the derivatives of .psi..sub..nu..sup.a are evaluated in the conventional finite-element manner, ##EQU27## The angle .theta., for the functions sin .theta. and cos .theta., is also computed from Eq. (7).

The radial integrals, R.sub..mu.'.mu..sup.(i), i=1, . . . ,5, .mu.',.mu.=1, . . . ,m (cf. Table I), are ##EQU28## which are Fourier sine and cosine transforms that can be evaluated by standard algorithms available in many mathematics software packages, ##EQU29## which can be evaluated with the same software used for I.sub..beta., ##EQU30## and, from Eq. (44), a.sub..mu..alpha. =0 for .alpha.>m+1. ##EQU31## and h.sub..mu..alpha. are determined by Eq. (20).

The tensor products A.sub..nu.'.nu. R.sub..mu.'.mu. in Eq. (34) are N.times.N matrices that are constructed by multiplying each term in the m.times.m [R] matrix by the entire n.times.n [A] matrix, as shown in Eq. (47). Above and to the left of the matrix are tables, reproduced from Table I, showing the relationship of the indices .mu. and .nu. to j, and .mu.' and .nu.' to i: ##STR1##

Eq. (33), which is the element matrix equation, and its supporting Eqs. (34)-(47) are the equations cited in box 20 of FIG. 1. These equations are sufficient to implement this invention. They may be coded into any scientific programming language; FORTRAN is currently the most popular choice.

Some practical observations in implementing the above equations into software are as follows. The radial integrals in Eq. (39) are identical for every infinite element in a mesh (because they are independent of angular variables and are along identical prolate radial paths, i.e., confocal hyperbolas emanating from the same spheroid), so they only need to be evaluated once for a given problem; their computational cost is totally insignificant. Hence, the numerical integration required to generate [K.sup..infin. ] and [M.sup..infin. ] for each infinite element involves only the evaluation of the 2-D angular integrals, making these 3-D elements as cheap to generate as 2-D elements. In addition, since the frequency dependence of the element is contained only in the radial integrals, element generation during a frequency sweep is essentially free after the first frequency.

To demonstrate the power of the inventive method, we calculated a structural acoustics problem using both the exemplary IEM code and, for comparison, a conventional Boundary Element Method (BEM) code. The structure was a long (aspect ratio of 10), cylindrical, thin, steel shell with tapered end caps and no internal structure, surrounded by water. The shell was excited by an obliquely incident plane wave at ka=1. (The symbol a represents the radius of the cylindrical shell.)

FIG. 13 shows the total elapsed solution time for a sequence of four progressively refined meshes for each of the two comparative methods. The solution time represents matrix generation plus equation solving, but does not include exterior field evaluation.

Since this system has three orthogonal planes of symmetry, it was sufficient for meshes 1, 2, and 3 to model a single octant. Consequently, it was possible to treat the obliquely incident wave by performing four analyses of the octant, each analysis having different symmetric or antisymmetric boundary conditions on the symmetry planes, and then superposing the results of the analyses.

Mesh 4 doubled mesh 3, thus covering a quadrant of the structure. This would be needed for a structure having only two planes of symmetry. This had the effect of roughly tripling the BEM cost, but it did not increase the IEM cost. The number of degrees of freedom (DOF) per analysis doubled for both methods while the number of analyses was halved. However, the rms half bandwidth B.sub.rms, which almost doubled for the BEM, remained almost constant for the IEM.

The times reported in FIG. 13 are for the maximum possible number of symmetric and antisymmetric analyses for the most general asymmetric excitation, i.e., 8.times.the run time for one-octant analysis, and 4.times.the run time for one-quadrant analysis. Run times for meshes 3 and 4 for the BEM are estimated times; estimates of times for meshes in a refinement sequence have consistently been accurate to within a few percent.

It is evident from the figure that the BEM is not substantially more expensive than the IEM for problems having less than a few thousand DOF. However, the IEM rapidly becomes more economical as the problems grow larger.

It should be noted that in order to practice the IEM, it is necessary to generate a 3-D mesh in the fluid surrounding the structure. This is a relatively simple matter for structures with convex outer surfaces. For example, the exemplary code typically requires only two or three commands to generate the entire fluid mesh for a structure having such a surface. We simply project outward the 2-D mesh that is on the fluid-contacting surface of the structure. For more general structural geometries, 3-D meshes can be generated by a wide variety of commercially available mesh generation programs.

We have tested the accuracy of the IEM against theoretical analyses and other structural acoustics programs such as our BEM program (which uses the same structural code as our IEM program). We have run a wide variety of scattering and radiation problems on both the IEM and BEM codes, and have achieved very close agreement in all cases. In the comparative IEM runs, we terminated the mesh refinements in the fluid when pointwise differences in acoustic pressure and structural velocity (relative to the BEM) were, typically, about 1%. We refer to this level of agreement as "graphical accuracy." For more rigorous studies, we terminated when the rms difference of all nodal values of pressure and normal velocity on the fluid-contacting surface of the structure was, typically, between 0.1% and 1%.

For the problem of thin cylindrical shells described above, the following fluid meshes provided errors of 0.1%-1%: For scattering problems, two layers of fluid elements, quadratic in the radial direction, surround the structure, as depicted in FIG. 10. The outer surface of the first layer is parallel to the structural surface and lies at a distance of 3.lambda./32, where .lambda. is the characteristic wavelength in the fluid at the excitation frequency .omega.. The outer surface of the second layer is the infinite element spheroid, and it is separated from the inner parallel surface by 13.lambda./32 on both the semi-major and semi-minor axes. Hence, the infinite element spheroid is, on average, .lambda./2 from the structure. The infinite elements use only a dipole, with r.sub.2 =2r.sub.1. (The significance of the distances r.sub.1, r.sub.2, etc. is readily understood by reference to FIG. 7.)

Of course, quadrupoles, using r.sub.3 =3r.sub.1, and octupoles, using r.sub.4 =4r.sub.1, may be expected to yield greater accuracy than dipoles. In general, the accuracy of the infinite element appears to be independent of the radial location of the layers of multipole nodes. In this regard, we found that model problems were insensitive, to at least five significant figures, to a wide range of values of r.sub.1, r.sub.2, . . . ,r.sub.m.

The meshes used for scattering problems can also be used to solve radiation problems excited by point or line forces on the shell, as depicted, for example, in FIGS. 14 and 15. However, this class of radiation problems calls for somewhat greater radial refinement within the infinite element spheroid near the places where the force is applied. We found that this requirement is satisfied, and an even smaller error is everywhere obtained, by replacing the two layers of radially quadratic elements by a single layer of radially quartic elements. (Both meshes use the same number of DOF.) Our practical experience suggests that the most efficient fluid mesh for this class of radiation problems, in the 0.1%-1% error range, is a single layer of radially cubic elements, with quartics near point or line forces, combined with dipole infinite elements.

1. Formal Expressions

The element matrices are derived by starting with a finite-size element, i.e., with the outer face on a spheroid of prolate radius r, and then taking the limit as r.fwdarw..infin.. Thus, applying the Galerkin weighted residual method to Eq. (8) over a single element yields ##EQU32## using k.sup.2 =.omega..sup.2 .rho./B. The first integral is converted into a surface integral and another volume integral using the identity .psi..sub.i .gradient..sup.2 p=.gradient..(.psi..sub.i .gradient.p)-.gradient..psi..sub.i. .gradient.p and the divergence theorem. Substituting Eq. (11A) into the volume integrals yields the following element matrix equation:

where the stiffness matrix, [K], mass matrix, [M], and pressure gradient vector, {F }, are, respectively, ##EQU33##

The surface integral for {F} is over the entire boundary of the element, S.sup.(e). It is split into two integrals: one over the outer face, S.sup.(e), and the other over the remaining faces, S.sup.(e) -S.sup.(e), ##EQU34##

Consider the first integral in Eq. (26A). As r.fwdarw..infin., prolate spheroidal coordinates approach spherical coordinates so .differential.p/.differential.n.fwdarw..differential.p/.differential.r and dS.fwdarw.r.sup.2 sin .theta.d.theta.d.phi., where r,.theta.,.phi. are spherical coordinates. To evaluate .differential.p/.differential.r as r.fwdarw..infin., substitute Eqs. (11A), (12A) and (16A) into Eq. (10A), which yields a stronger form of the Sommerfeld condition, ##EQU35## where the upper case o, read "big o", means "at least as fast as". Substitute .differential.p/.differential.r in Eq. (27A) for .differential.p/.differential.n in the first integral. The O(1/r.sup.2) term makes no contribution to the integral because .psi..sub.i is O(1/r) and dS.varies.r.sup.2 as r.fwdarw..infin.. In the remaining term, ikp, substitute Eq. (11A) for p.

In the second integral in Eq. (26A), substitute the balance of linear momentum for .differential.p/.differential.n, viz., .differential.p/.differential.n=.omega..sup.2 .rho.u.sub.n, where u.sub.n is the amplitude of the normal component of particle displacement.

Making these substitutions, Eq. (26A) becomes ##EQU36##

The [C] matrix is the radiation damping matrix, representing radiation loss to infinity.

The {D} vector permits specification of u.sub.n on the side or bottom faces of the element. These faces are either on the boundary of, or interior to, the domain.

Face on boundary of domain: A nonzero value corresponds to a moving boundary with a known normal displacement (or velocity) amplitude. Along the side faces, a nonzero u.sub.n must decay fast enough in the radial direction for the integral in Eq. (29A) to be finite.

A zero value corresponds to a rigid boundary, i.e., motion permitted parallel to the boundary but not perpendicular to it. Perhaps the most common occurrence of a rigid infinite boundary is when there is a plane of symmetry through the structure and fluid, permitting the analyst to model only the structure and fluid-on one side of the plane, i.e., an infinite half-space problem. A plane of symmetry occurs when all passive system properties (geometry, constraints, material properties and, if present, hydrostatic prestresses) are symmetric with respect to a plane. If the excitation is also symmetric, or just the symmetric component of an asymmetric excitation is being considered, then the system response must possess the same symmetry. Therefore, vector quantities must be a mirror reflection across the plane. In the case of particle displacement, this means that the normal component must vanish, in order to preserve material continuity.

Face in interior of domain: Every interior face is common to two adjacent infinite or finite-size acoustic elements. Assembly of the {D} vectors from two adjacent elements yields, on the interface, a sum of two surface integrals with opposite sign because the outward normals are in opposite directions. The pair of integrals therefore represents the difference of u.sub.n across the interface. A nonzero difference corresponds to a jump discontinuity in the normal displacement, i.e., a separation of the fluid. A zero difference corresponds to preservation of material continuity, the usual condition in practical applications.

Substituting Eq. (28A) into Eq. (24A) yields for the element matrix equation:

where, summarizing, the stiffness, mass and radiation damping matrices are, respectively, ##EQU37## The remainder of Section IIID transforms the integrals in Eq. (31A) to expressions that can be numerically evaluated. (The author's IEM code has so far only used zero values for D.sub.i in the infinite elements in all applications.)

2. Transformation of Integrals to Prolate Spheroidal Coordinates

The differential volume and surface elements are dV=J.sub..nu. drd.theta.d.phi. and dS=J.sub.S d.theta.d.phi., where J.sub..nu. is the volume Jacobian and J.sub.S is the surface Jacobian (on spheroids). Since the faces of the infinite element conform to spheroids (.theta.,.phi.-surfaces) and hyperbolas (r-lines), the integration limits for the volume integrals can be separated into angular limits and radial limits. Thus, Eq. (31A) becomes ##EQU38## where .sigma..sup.(e) is the "spheroidal cross section" of the element, i.e., any confocal spheroidal surface inside the element and bounded by the side faces. There is only one spheroidal cross section in r,.theta.,.phi.-space because the .theta.,.phi. coordinates of the boundary of the cross section are independent of r. Hence, S.sup.(e) is equivalent to .sigma..sup.(e) in the surface integration for C.sub.ij. The transformation of d.sub.V, d.sub.S and .gradient. to prolate spheroidal coordinates is as follows.

Volume Jacobian ##EQU39## where d=xi+yj+zk is the position vector and a.sub.r,a.sub..theta.,a.sub..phi. are unit vectors, with dimensions of length, in the r,.theta.,.phi. directions, respectively. For orthogonal coordinates, a.sub.r.a.sub..theta. .times.a.sub..phi. =1. The scale factors h.sub.r, h.sub..theta. and h.sub..phi. are given by ##EQU40## Using Eq. (5A) yields ##EQU41## and Eq. (33A) becomes

Surface Jacobian

J.sub.S is needed in C.sub.ij only on the outer face, i.e., on the spheroid at r=r, ##EQU42## and it is needed only in the limit as r.fwdarw..infin. (while f remains constant). Hence, ##EQU43## which is the same as the surface Jacobian in spherical coordinates because prolate spheroidal coordinates.fwdarw.spherical coordinates as r.fwdarw..infin..

Gradient ##EQU44## The gradient is needed only in K.sub.ij, where it occurs as V.psi..sub.i..gradient..psi..sub.j J.sub.V. Using Eqs. (35A), (36A) and (39A) yields ##EQU45##

Eqs. (40A), (36A) and (38A) can now be substituted into K.sub.ij, M.sub.ij and .sub.ij, respectively, in Eq. (32A). In addition, Eq. (12) is substituted for .psi..sub.i and .psi..sub.j, introducing primes to distinguish the indices .nu. and g that are associated with the index i:

Performing these substitutions, and noting that the angular shape functions are implicitly functions of .theta.,.phi. [because of the coordinate mapping in Eq. (7)], all the 3-D integrals separate into products of 2-D "angular" integrals, A.sub..nu.'.nu..sup.(i), and 1-D "radial" integrals, R.sub..mu.'.mu..sup.(i). The 2-D surface integral for C.sub.ij separates in a similar manner. The resulting expressions are as follows. ##EQU46## The symbols R.sub..mu.'.mu..sup.(1) and R.sub..mu.'.mu..sup.(5) include an overbar because a term will later be removed from each integral [cf. Eqs. (58A) and (73A)] and the remaining terms will be labeled R.sub..mu.'.mu..sup.(1) and R.sub..mu.'.mu..sup.(5) so that all five radial integrals will have uniform symbols (no overbars) in the final expressions [Eq. (82A)].

The tensor products A.sub..nu.'.nu. R.sub..mu.'.mu. are N.times.N matrices that are constructed by multiplying each term in the m.times.m [R] matrix by the entire n.times.n [A] matrix, as shown in Eq. (46A). Above and to the left of the matrix are tables, reproduced from Table II, showing the relationship of the indices .mu. and .nu. to j, and .mu.' and .nu.' to i. ##STR2##

3. Angular Integrals

Because of the coordinate mapping in Eq. (7), d.theta.d.phi.=Jd.xi.d.eta., where the Jacobian is ##EQU47## Hence Eqs. (43A) become, ##EQU48##

All four integrals can be numerically integrated using Gauss rules since the integrands are smooth and bounded (including the 1/sin .theta. term because J.varies..theta. as .theta..fwdarw.0). The derivatives of .psi..sub..nu..sup.a are evaluated in the conventional manner, ##EQU49## Indeed, the entire evaluation of the angular integrals follows conventional finite element analysis procedures.

4. Radial Integrals

R.sub..mu.'.mu..sup.(1)

From Eq. (16A), ##EQU50## and, from Eq. (52A), a.sub..mu..alpha. =0 for .alpha.>m+1.

Substituting Eq. (53A) into the first of Eqs. (44A) yields ##EQU51## A factor k is introduced into the integrals (and cancelled by the 1/k in L.sub..mu.'.mu.) to make them dimensionless. The b.sub.2 integral is separated from the others because the upper limit yields an undefined term (which will later be cancelled by similar terms from M.sub.ij and C.sub.ij). Formally, we shall write ##EQU52## The other integrals are well defined, so the limiting process only requires changing the upper limit from r to .infin.. Eq. (55A) becomes ##EQU53##

The I.sub..beta. integrals can be evaluated two ways.

(i) Sine and Cosine Integrals ##EQU54## where Si and Ci are the well known sine and cosine integrals. For .beta.>1, integration of I.sub..beta. by parts yields the recursion relation ##EQU55##

(ii) Fourier Sine and Cosine Transforms

Making the change of variable

Eq. (61A) becomes ##EQU56##

The author used the second approach in his IEM code because of the availability of a Fourier sine and cosine transform routine, which was also used for the more complicated J.sub..mu. integral for R.sub..mu.'.mu..sup.(4) below.

R.sub.m'm.sup.(2)

From Eq. (16), ##EQU57## and h.sub..mu..alpha. =0 for .alpha.>m. Substituting Eq. (66A) into the second of Eqs. (44A) and letting the upper limit be .infin., since all integrals are well defined, yields ##EQU58##

R.sub.m'm.sup.(3)

Substituting Eq. (53A) into the third of Eqs. (44A) and letting the upper limit be .infin., since all integrals are well defined, yields ##EQU59##

R.sub.m'm.sup.(4)

Substituting Eq. (66A) into the fourth of Eqs. (44A) and letting the upper limit be .infin., since all integrals are well defined, yields ##EQU60## The denominator does not vanish because [(kr).sup.2 -.epsilon..sub.1.sup.2 .zeta..sup.2 ]=k.sup.2 (r.sup.2 -f.sup.2) and r.ltoreq.r.sub.1 >f. Using Eq. (64A) converts J.sub..mu. to Fourier sine and cosine transforms, ##EQU61## which can be evaluated with the same routine used for the I.sub..beta. integrals in Eq. (65A).

R.sub.m'm.sup.(5)

Substituting Eq. (66A) into the fifth of Eqs. (44A) and proceeding in the same manner as for R.sub..mu.'.mu..sup.(1) (viz. separating the c.sub.2 integral from the others because the upper limit produces an undefined oscillatory term) yields ##EQU62##

5. Radial Functions on Outer Face

The radiation damping matrix, C.sub.ij, in Eq. (42A) includes the term ##EQU63## Using Eq. (66A), ##EQU64##

Taking the limit as r.fwdarw..infin. eliminates all the c.sub..beta. terms except c.sub.2, which contains an undefined oscillatory term: ##EQU65##

6. Vanishing of Undefined Oscillatory Terms

The expressions for all the radial integrals [Eqs. (58A), (59A), (68A), (69A), (70A), (73A) and (74A)] and the radial functions on the outer face [Eq. (76A)] can now be substituted into the stiffness, mass and damping matrices in Eq. (42A). Adding those three matrices together as in the element matrix equation, Eq. (30A), yields

where [K.sup..infin. ] and [M.sup..infin. ] comprise all the well-defined integrals and are summarized in the next section. The "undefined oscillatory terms" comprise all terms containing the expression lim e.sup..sup.-i2kr. There are three such occurrences: Eqs. (58A), (73A) and (76A), which contribute to [K], [C] and [M], respectively. Thus, ##EQU66## using .omega.=ck and c.sup.2 =B/.rho.. From Eqs. (52A) and (54A), b.sub.2 =a.sub..mu.'1 a.sub..mu.1 =-h.sub..mu.'1 h.sub..mu.1, and from Eq. (67A), c.sub.2 =h.sub..mu.'1 h.sub..mu.1 =-b.sub.2. Hence,

In words, the stiffness and mass matrices consist of two types of terms: those that are independent of the location of the outer face (the "well-defined" integrals) and those that do depend on its location and therefore oscillate as the face recedes to infinity. The damping matrix, which represents application of the Sommerfeld radiation condition to the outer face, is completely oscillatory. Eq. (79A) states that the radiation condition exactly cancels the oscillatory terms in the stiffness and mass matrices.

7. Final Expressions

From Eqs. (77A) and (79A) it follows that the final form of the element matrix equation for the infinite acoustic element is ##EQU67## The vector {D} is given in Eq. (29A); it is zero in almost all practical applications. The angular integrals, A.sub..nu.'.nu..sup.(i), i=1,2,3,4, are given in Eq. (48A). The radial integrals are ##EQU68## and L.sub..mu.'.mu., I.sub..beta., J.sub..beta., b.sub..beta., c.sub..beta., and .zeta. are given by Eqs. (56A), (61A), (71A), (54A), (67A) and (60A), respectively. Eqs. (54A) and (67A), in turn, need h.sub..mu..alpha., which are determined by Eq. (20).

The radial integrals in Eq. (82A) are identical for every infinite element in a mesh (because they are independent of angular variables and are along identical prolate radial paths, i.e., hyperbolas emanating from the same spheroid), so they only need to be evaluated once for a given problem; their computational cost is totally insignificant. Hence, the numerical integration required to generate [K.sup..infin. ] and [M.sup..infin. ] for each infinite element involves only the evaluation of the 2-D angular integrals, making these 3-D elements as cheap to generate as 2-D elements. In addition, since the frequency dependence of the element is contained only in the radial integrals, element generation during a frequency sweep is essentially free after the first frequency.