Wave boundary elements: a theoretical overview presenting applications in scattering of short waves

Abstract

It is well known that the use of conventional discrete numerical methods of analysis (FEM and BEM) in the solution of Helmholtz and elastodynamic wave problems is limited by an upper bound on frequency. The current work addresses this problem by incorporating the underlying wave behaviour of the solution into the formulation of a boundary element, using ideas arising from the Partition of Unity finite element methods. The resulting ‘wave boundary elements’ have been found to provide highly accurate solutions (10 digit accuracy in comparison with analytical solutions is not uncommon). Moreover, excellent results are presented for models in which each element may span many full wavelengths. It has been found that the wave boundary elements have a requirement to use only around 2.5 degrees of freedom per wavelength, instead of the 8–10 degrees of freedom per wavelength required by conventional direct collocation elements, extending the supported frequency range for any given computational resources by a factor of three for 2D problems, or by a factor of 10–15 for 3D problems. This is expected to have a significant impact on the range of simulations available to engineers working in acoustic simulation. This paper presents an outline of the formulation, a description of the most important considerations for numerical implementation, and a range of application examples.

Introduction

The use of numerical methods has become widespread in the analysis of a wide variety of wave problems. We focus here on the analysis of the propagation of acoustic and elastodynamic waves, by radiation or scattering, through isotropic media. Starting with the acoustic wave case, this is governed in the frequency domain by the well-known Helmholtz equation

$$\text{∇}{}^2\text{φ+k}{}^2\text{φ=0}$$

where φ is a potential function whose complex nature provides information on both the magnitude and phase of oscillation. The wave number, k, is given by

$$\text{k=}\text{2}\text{π}\text{λ}$$

where λ is the wavelength, and therefore the wave number is seen as a mathematically convenient way of expressing the frequency, ω, under consideration. The potential φ expresses the spatial variation of the wave, and may be simply multiplied by e^−iωt to arrive at the time-varying field.

In order to present the new work in wave boundary elements, it is necessary first to present a brief overview of the mathematical formulation of the conventional direct collocation boundary element method [1]. Here, the potential and its normal derivative on the boundary are expressed in a weak form integral equation

$$\text{φ(x)}\text{2}\text{+}\text{∫}\text{Γ}\text{∂}\text{G(x,y)}\text{∂}\text{n}\text{φ(y)}\text{d}\text{Γ}{}_{\mathrm{y}}\text{=}\text{∫}\text{Γ}\text{∂}\text{φ(y)}\text{∂}\text{n}\text{G(x,y)}\text{d}\text{Γ}{}_{\mathrm{y}}$$

where x and y are, respectively, the locations of a fictitious source and a ‘field’ point, both on the boundary, Γ, to the domain. Γ is assumed here to be smooth. The term G(x,y) is the so-called fundamental solution, or free-space Green's function, describing the problem. For 2D problems this is given by

$$\text{G(x,y)=}\text{i}\text{4}\text{H}{}_0\text{(kr)}$$

where r is the distance from point x to point y, i is the imaginary square root of −1 and H₀ is a Hankel function of the first kind of order zero. The fundamental solution is simpler in 3D, being given by

$$\text{G(x,y)=}\text{e}{}^{\mathrm{<mtext>i</mtext><mtext>kr</mtext>}}\text{4}\text{π}\text{r}$$

Use of these fundamental solutions automatically satisfies the Sommerfeld radiation condition at infinity. Eq. (3) may be integrated numerically by subdividing the boundary into elements, over each of which the potential may be expressed in terms of nodal values according to conventional interpolation using Lagrangian shape functions, i.e.

$$\text{φ=}\text{N}{}^{\mathrm{<mtext>T</mtext>}}\text{Φ}$$

where $\text{N}$ is a vector containing the shape functions and Φ is a vector containing the nodal potentials. The superscript ‘T’ denotes transpose of a matrix or vector. The derivatives of the nodal potentials may be similarly interpolated. By placing the point x at each node in turn, a system of equations is produced in the nodal potentials that may be solved using a direct or iterative solver.

For scattering problems involving a known incident wave, the potential is often expressed as the summation of the incident potential, φⁱ, and the scattered potential, φ^s, and the incident potential is added to the right hand side of Eq. (3). Furthermore, for exterior radiation or scattering problems, some numerical consideration needs to be given to the elimination of the irregular frequencies, corresponding to the eigenvalues of an acoustic cavity bounded externally by Γ, that give rise to a non-uniqueness problem in the solution. This may be done using the approach of Burton and Miller [2] or using the CHIEF algorithm of Schenck [3] (or one of its derivatives).

Although this approach has found widespread use in acoustic analysis in industry, the method is limited in the frequency range that may be considered. The limit is imposed by constraints on computational resources and run time. The limit arises from the subdivision of the boundary, Γ, into elements, and the assumption embodied in Eq. (6) that over each element the potential varies according to some polynomial, often quadratic. This requires that each element be sufficiently small, in comparison with the wavelength, λ, that a quadratic polynomial is capable of capturing the required portion of the sinusoidal waveform. It is generally accepted that, in order to provide results of reasonable engineering accuracy, four quadratic elements (or 8–10 nodes) are required per full wavelength.

For problems involving physically large boundaries and/or short wavelengths, this constraint quickly causes the required computational resources to become substantially larger than those available today, and has the effect of imposing an upper bound on the frequency that may, in practice, be considered for any problem. For higher frequencies an alternative approach, such as Statistical Energy Analysis (SEA), is commonly employed. However, there is a lower frequency bound to the applicability of SEA, below which its assumption of high modal density becomes invalid. This leaves engineers in some disciplines with an awkward mid-range of frequencies in which it is difficult to obtain reliable computational predictions. One such field is the acoustics of automotive passenger compartments, in which the mid-frequency range may be found in the hundreds of Hz, which is an important part of the audible range.

Approaches that attempt to address this issue have recently started to emerge. The Fast Multipole Method (FMM) of Rokhlin [4] has been adopted by a number of researchers [5], [6] to enable rapid solution of the large systems that are required for short wave analysis of acoustics and electromagnetics. Other authors have developed fast iterative solvers, [7], or attempt to sparsify the system matrix using wavelet compression. The approach developed in this work involves the use of a plane wave basis to define the potential. Our approach, which has the benefit of greater simplicity, addresses the problem by reducing the size of the matrix rather than by attempting to introduce sparsity. It may well be that the benefits, presented herein, of the wave boundary elements may be combined with the benefits of the FMM, i.e. the authors believe that the two methods are not mutually exclusive.

Elastodynamic waves are governed by the harmonic form of the wave equation, so that (in the absence of body forces) the displacement field, $\text{u}\text{,}$ (a complex valued vector field) is sought as the solution to

$$\text{μ}\text{∇}{}^2\text{u}\text{+(λ+μ)}\text{∇∇}\text{·}\text{u}\text{+ρω}{}^2\text{u}\text{=0}$$

where λ (not to be confused with the wavelength) and μ are the Lamé material constants, and ρ is the mass density. The situation is more complicated than the acoustic case, not only because of the vector form of the solution, but also since there are different types of wave being propagated, and these have different speeds (and therefore different wave numbers) at the same frequency. We consider the pressure (dilatational) P-wave, having wave number k_p, and the shear (distortional) S-wave, having wave number k_s. These wave numbers arise from the material properties of the medium under consideration, being given by

$$\text{k}{}_{\mathrm{<mtext>p</mtext>}}\text{=ω}\text{ρ}\text{2μ+λ}\text{and}\text{k}{}_{\mathrm{<mtext>s</mtext>}}\text{=ω}\text{ρ}\text{μ}$$

so that typically k_s≈2k_p for a wide range of materials. Rayleigh waves are not considered in the present work.

Conventional mathematical treatment [8] allows the governing differential equation to be cast in the integral form

$$\text{u}\text{(}\text{x}\text{)}\text{2}\text{+}\text{∫}\text{Γ}\text{T}\text{(}\text{x}\text{,}\text{y}\text{)}\text{u}\text{(}\text{y}\text{)}\text{d}\text{Γ}{}_{\mathrm{y}}\text{=}\text{∫}\text{Γ}\text{U}\text{(}\text{x}\text{,}\text{y}\text{)}\text{t}\text{(}\text{y}\text{)}\text{d}\text{Γ}{}_{\mathrm{y}}\text{+}\text{u}{}^{\mathrm{<mtext>i</mtext>}}\text{(}\text{x}\text{)}$$

where $\text{T}$ and $\text{U}$ are the Stokes' traction and displacement tensors that serve as fundamental solutions for this problem, and $\text{t}$ is the surface traction vector. The superscript ‘i’ refers to the incident wave as for the acoustic case, and the left hand side integral is calculated in the Cauchy principal value sense. The solution procedure is thenceforth identical in concept to that described for the Helmholtz problem.