Wave boundary elements: a theoretical overview presenting applications in scattering of short waves
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 equationwhere φ is a potential function whose complex nature provides information on both the magnitude and phase of oscillation. The wave number, k, is given bywhere λ 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 equationwhere 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 bywhere r is the distance from point x to point y, i is the imaginary square root of −1 and H0 is a Hankel function of the first kind of order zero. The fundamental solution is simpler in 3D, being given byUse 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.where 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, φi, 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, (a complex valued vector field) is sought as the solution towhere λ (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 kp, and the shear (distortional) S-wave, having wave number ks. These wave numbers arise from the material properties of the medium under consideration, being given byso that typically ks≈2kp 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 formwhere and are the Stokes' traction and displacement tensors that serve as fundamental solutions for this problem, and 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.
Section snippets
Wave boundary elements
The approach described herein is to include the wave behaviour into the shape function of the element. This has been considered for Helmholtz problems by several authors, although not normally in a boundary element context but instead using infinite and finite elements. Astley [9] and Chadwick and Bettess [10] developed wave envelope methods, in which it is the shape of the envelope and not that of the potential that becomes the primary variable. Later, Melenk and Babuska [11] introduced the
Some numerical considerations
The approach described in the previous section has been implemented in computer codes for 2D analysis of both Helmholtz and elastodynamic wave propagation. A further code has been written to consider the 3D acoustic problem [25]. Before presenting engineering applications and a comparison with conventional BEM analysis, it is important to mention some of the numerical features of such an implementation.
We recall, firstly, that the use of the plane wave basis causes the governing matrix to
Applications and error analysis
The wave boundary elements have been used to study a variety of classical problems, with known analytical solutions, and other problems involving geometries that are selected to test the numerical stability of the formulation.
Conclusions
This paper has presented the basic theoretical foundations for wave boundary elements in both Helmholtz and elastodynamic wave problems, along with some considerations relating to the BEM implementation and numerical stability of the solutions. Error analysis for simple scattering problems shows the new elements to perform extremely well in comparison with conventional, direct collocation, quadratic boundary elements, exhibiting accuracy improvements of up to 8 orders of magnitude and extending
Acknowledgements
This work was carried out under grant GR/N09879 from the Engineering and Physical Sciences Research Council (EPSRC), UK, whose support is gratefully acknowledged. Peter Bettess is also grateful to EPSRC for support under grant SF/000169.
References (31)
Rapid solution of integral equations of classical potential theory
J Comput Phys
(1985)The fast multipole method: numerical implementation
J Comput Phys
(2000)- et al.
Mapped wave envelope elements for acoustical radiation and scattering
J Sound Vib
(1994) - et al.
The partition of unity finite element method. Basic theory and applications
Comput Meth Appl Mech Engng
(1996) - et al.
Computational aspects of the ultra weak variational formulation
J Comput Phys
(2002) - et al.
Plane wave approximation in direct collocation boundary element method for radiation and wave scattering by smooth objects: numerical aspects and applications
J Sound Vib
(2003) - et al.
Boundary element methods in acoustics
(1991) - et al.
The application of integral equation methods to the numerical solution of some exterior boundary value problems
Proc R London Ser A
(1971) Improved integral formulation for acoustic radiation problems
J Acoust Soc Am
(1968)- et al.
Fast solution methods in electromagnetics
IEEE Trans Antennas Propag
(1997)