This section presents the formulation of theinfinite element method(IEM), a technique for handling unbounded computational acoustic problems in the finite element framework. In the first part of this section the multipole expansion and the Atkinson – Wilcox theorem are discussed. Then, the infinite element formulation utilizing the Astley – Leis mapped approach is explained.
4.3.1 Multipole expansion
The solution of the Helmholtz equation in three dimensions in closed and open domains can be expanded in terms of separable solutions in spherical coordinatesx(r, θ, ϕ)as
p(x, ω) =
withPnmrepresenting Legendre polynomials,h(1)n andh(2)n denoting Hankel functions of the first and second kind, respectively, andAnm,Bnm,Cnm,Dnmbeing constants. The Hankel functions can be expaneded as When only outwardly propagating waves are assumed, which means that the Sommerfeld radiation conditon (4.19) is fulfilled, the second summation in equation (4.20) can be omitted and
4.3. INFINITE ELEMENTS 43
Figure 4.1.Illustration of the Atkinson – Wilcox criterion. PointAsatisfies the criterion, while pointBdoes not. After Fig. 7.4 of [13].
the expansions of the Hankel functions can be rearranged to give
multipole
withGn(θ, ϕ, ω)denoting a directivity function associated with thenth inverse power ofr. This form is called amultipole expansion. The expansion can also be interpreted in terms of near and far field contributions. The leading term, associated withr−1, determines the far field directivity, while the remaining terms contribute to the near and intermediate field. If expression (4.22) is truncated—as it will be the case when trial solutions are defined for infinite elements—the num-ber of terms retained determines the extent to which the truncated expansion is able to resolve near field effects.
The multipole Atkinson –
Wilcox theorem expansion (4.22) is expected to hold in the far field, but the extent of its validity
in the near field is not obvious. This question is answered by the Atkinson – Wilcox theorem [145].
It states that the sound field at any point that lies entirely outside a circumscribing sphere, which itself encloses all radiating and scattering sources, can be written as a multipole expansion of type (4.22), and that this expansion is absolutely and uniformly convergent.
The principle is illustrated in Figure 4.1, with the point A lying outside a circumscribing sphere of radiusaand therefore satisfies the criterion, whereas point B does not. The Atkin-son – Wilcox criterion was first derived for a spherical coordinate system [145] (Figure 4.1(a)), and a similar result was found later [36] for spheroidal or ellipsoidal coordinate systems (Fig-ure 4.1(b)). The inner region, in which the multipole expansion cannot be used, can be greatly reduced by using a spheroidal rather than a spherical surface when modeling oblong objects.
4.3.2 Mesh, mapping, shape and test functions
When a mesh infinite
element mesh is created using infinite elements, the infinite elements are attached to the artificial boundaryΓ∞, as depicted in Figure 4.2. The whole unbounded domain is subdivided into a finite number of elements of infinite size, so that (4.14) holds. The infinite elements must be defined so that their inner faces overlap with the outer faces of the finite elements located atΓ∞. The region occupied by the infinite elements is denoted byΩ∞in the following.
Ω∞
Figure 4.3.Coordinate mapping for infinite elements
There are a large number of formulations of infinite elements, as summarized in [13]. The infi-nite elements can be defined over the global coordinate system (unmappedorseparable formulation, see e.g. [36, 131]), or can be defined over common finite parent domains in the local coordinate system (mapped formulation). In the following our discussion is restricted to the Astley – Leis con-jugated formulation with mapped elements.
The basis of mapping infinite element formulations is the coordinate transform, which maps the element from a standard, finite domain into an infinite region of the physical domain. The mapping is carried out similar to equation (4.15), and is depicted in Figure 4.3. This finite to infinite transformation is performed by the proper definition of geometry shape functionsL(ξ) in the standard parent domain Oe. It can be seen that the simple one-dimensional mapping function mapped fromξ= 0. The multiplicators ofx1andx2in equation (4.23) are the geometrical shape functions of a simple line infinite element.
The mapping can directly be extended into two or three dimensions, as illustrated in Fig-ure 4.3 for a quadrilateral element. The element geometry is defined by themapping nodes, while thepressure nodes, corresponding to the DOFs of the solution, can be defined independently from the mapping nodes. The interpolating shape functions for the pressure field are defined in a different manner compared to standard finite elements: they are chosen so that they inherently follow the oscillating behavior of the solution [13, 15, 31]. This is achieved by inserting a phase term into the standard shape functions, represented by thephase functionµ(x), thus the approxi-mation of the pressure, in analogy with equation (4.8) becomes
pressure
withSidenoting the standard shape function on the inner face of the element,nf the number of mapping nodes on this face, andai representing the distance of these mapping nodes from an imaginary source point located insideΩ(see Figure 4.3). The functiona(η, ζ) =P
iSi(η, ζ)aithat
4.3. INFINITE ELEMENTS 45 results of the interpolating summation in equation (4.25) is continuous inside the element and in the whole domainΩ∞.
The other part of the pressure interpolation functions Pi can also be expressed in the local coordinate system. Astleyet al.[15] define these shape functions by using the standard finite element shape functions in theη and ζ directions and the Lagrange polynomialsLmj in the ξ direction as
Pi(ξ) =Si(η, ζ)1−ξ
2 Lmj (ξ). (4.26)
The Lagrange polynomials are defined as
Lagrange
withξk denoting theξcoordinate of thekth pressure node of the infinite element. The number of differentξlocations for the pressure nodes of the element is denoted here bymand is also called as theradial orderof the element. One of the pressure nodes is located on theΓ∞surface, atξ=−1, to ensure the direct coupling of the FE and IE parts of the mesh. The other pressure nodes are distributed in an equidistant manner over the−1≤ξ <1region.
Such a choice of shape functions provides the following properties 1. Ni(ξ) = 0in all nodes of the element, except for nodei.
2. Ni(ξ)≡Si(η, ζ)on the face of the element.
3. Ni(ξ)is a polynomial of ordermmultiplied by the phase terme−jk(r−a).
The first two properties ensure a compatible matching with the inner finite element mesh. From the third property the radial behavior ofNican be assessed. With keepingηandζconstant, the local behavior along theξ-axis gives the global behavior
Ni(x)∼hα1 The absence of the constant (orderr0) term in equation (4.28) is due to the multiplicator(1−ξ)/2 in the definition of the pressure shape functions, see eq. (4.26). It can be seen that the above behavior follows the expansion for outgoing waves defined in equations (4.20) and (4.22) trun-cated to themth order. Thus, by defining an infinite element of radial orderm, an mth order approximation of the multipole expansion can be achieved.
4.3.3 The discretized form
In the Astley – Leis conjugated formulation the test functions (also called weighting functions) are defined as the complex conjugate of the pressure interpolation functionsNi, multiplied by the Astley – Leis weightD(x). The test functionsφi(x)therefore become
Astley –
The weightD(ξ)cancels out the∝ (1−ξ)−2 term arising in the Jacobian from the derivatives
∂/∂ξLi, see equation (4.23). The gradients of the shape and weighting functions defined by equa-tions (4.26) and (4.29) are obtained as
∇Ni=∇Pie−jkµ−jk∇µPie−jkµ (4.30)
∇φi=∇DPiejkµ+D∇Piejkµ+ jk∇µPiejkµ. (4.31)
It can be seen that the complex multiplicators of the weight and shape functions,ejkµande−jkµ, respectively, cancel out in all terms of the productsφiNjand∇φi· ∇Nj.
Substituting the definitions from the previous section into the discrtized weak form, the fol-lowing form of infinite element matrices
infinite
The final matrix equation of the discretized finite–infinite element problem becomes
Kp0+ jωCp0−ω2Mp0= jωAv. (4.33) Here the matrices K, C, and M already incorporate the finite element contributions defined in equation (4.13). The damping matrix C has only infinite element contributions. All three matrices are real valued and sparse; however, the Astley – Leis weighting introduces asymmetry.
In equation (4.33) the notationp0 is applied since the nodal values of the solution do not exactly correspond to the complex pressure amplitudes at the given nodes in case of infinite elements.
An issue matrix
conditioning
that needs to be taken into account in the numerical solution of the matrix equation (4.33) is the condition number of the system matrices; as using infinite elements of high radial orders can result inill-conditioned matrices. This issue is addressed in a number of publications, and is usually solved by choosing special orthogonal polynomials for the pressure interpolation functionsPi(x), see references [31, 51, 52, 74]. It is also worth mentioning that under special conditions (see [14]), similar to (4.12), equation (4.33) can be transformed and solved in the time domain in case of transient problems.