The Boundary Element Method in Acoustics: A Survey

The boundary element method (BEM) in the context of acoustics or Helmholtz problems is reviewed in this paper. The basis of the BEM is initially developed for Laplace’s equation. The boundary integral equation formulations for the standard interior and exterior acoustic problems are stated and the boundary element methods are derived through collocation. It is shown how interior modal analysis can be carried out via the boundary element method. Further extensions in the BEM in acoustics are also reviewed, including half-space problems and modelling the acoustic field surrounding thin screens. Current research in linking the boundary element method to other methods in order to solve coupled vibro-acoustic and aero-acoustic problems and methods for solving inverse problems via the BEM are surveyed. Applications of the BEM in each area of acoustics are referenced. The computational complexity of the problem is considered and methods for improving its general efficiency are reviewed. The significant maintenance issues of the standard exterior acoustic solution are considered, in particular the weighting parameter in combined formulations such as Burton and Miller’s equation. The commonality of the integral operators across formulations and hence the potential for development of a software library approach is emphasised.


Introduction
Acoustics is the science of sound and has particular applications in sound reproduction, noise and sensing [1,2]. Acoustics may be interpreted as an area of applied mathematics [1,3]. In special cases, the mathematical equations governing acoustic problems can be solved analytically, for example if the equations are linear and the geometry is separable. However, for realistic acoustic problems, numerical methods provide a much more flexible means of solution [4,5]. Numerical methods are only useful in practice when they are implemented on computer. It is over the last sixty years or so that numerical methods have become increasingly developed and computers have become faster, with increasing data storage and more widespread. This brings us to the wide-ranging area of computational acoustics; the solution of acoustic problems on computer [6][7][8][9][10][11], which has significantly developed in this timescale.
In the context of this work, the acoustic field is assumed to be present within a fluid domain. If there is an obstacle within an existing acoustic field, then the disturbance it causes is termed scattering. If an object is vibrating and hence exciting an acoustic response, or contributing to an existing acoustic field, then this is termed radiation. If the fluid influences the vibration of the object or structure, and energy is exchanged in both directions between them, then this is termed coupled fluid-structure interaction [12] or, in the acoustic context, vibro-acoustics [13,14]. Alternatively, if the acoustic field is an outcome of a background flow, for example the generation of noise by turbulence, then this is termed aero-acoustics [11,15]. The determination of the properties of a vibrating or scattering object (for example, its shape, the surface impedance or the sound intensity) from acoustic measurements in the field is termed inverse acoustics [16].
Vibration analysis is not the focus of this work, but it clearly cannot be divorced from the study of acoustics. In vibration analysis, the domain of the structure oscillates, perhaps under variable loading or forcing. In the case of vibro-acoustics, the vibration can excite or be excited by the acoustic field and, in any analysis, both need to be modelled [17][18][19], and the equations coupled. The computational modelling of vibration is normally carried out by the finite element method (FEM) [20,21], and this method is significantly established in this application area. Vibrational modelling may be carried out in the time domain. However, with the nature of structural vibration in that it is normally dominated by modes and their corresponding resonant frequencies, and hence vibration is often close to periodic. Similarly, the excitation forces on a structure, such as the driver on a loudspeaker of the explosions in an engine can similarly be close to periodic in short time scales. It is, therefore, found that vibratory problems are routinely analysed in the frequency domain, both in the practical work and computational modelling.
An acoustic field is most straightforwardly interpreted as a sound pressure field, with the sound pressure varying over the extent of the domain, and with time. The transient acoustic field can be computationally modelled by the finite element method [22,23], the finite difference-time domain method (a particular finite difference method that was originally developed for electromagnetic simulation [24][25][26]) can also be applied in acoustics [27][28][29][30]) and the boundary element method [31][32][33][34], but acoustics and vibration are more often analysed and modelled in the frequency domain. In acoustic problems, the most likely fluid domains are air or water and, in many cases in these fluids, the linear wave equation is an acceptable model. By observing one frequency at a time, the wave equation can be simplified as a sequence of Helmholtz equations. Again, there are a variety of methods for the numerical solution of Helmholtz problems; the finite difference method, the finite element method [23,[35][36][37][38][39][40][41][42][43] and, the subject of this paper, the boundary element method.
The boundary element method is one of several established numerical or computational methods for solving boundary-value problems that arise in mathematics, the physical sciences and engineering [44]. Typically, a boundary-value problem consists of a domain within or outside a boundary in which a variable of interest or physical property is governed by an identified partial differential equation (PDE). A computational method for solving a boundary-value problem is tasked with finding an approximation to the unknown variable within the domain. Mostly, this is carried out by domain methods, such as the finite element method [45], in which a mesh is applied to the domain. However, the boundary element method works in an entirely different way in that in the first stage further information is found on the boundary only; the solution at the domain points is found in the second stage by integrating over the known boundary data. The boundary element method requires a mesh of the boundary only, and hence is generally easier to apply than domain methods. The number of elements in the BEM is therefore expected to be much less than in the corresponding finite element method (for the same level of required accuracy or element size), and there is therefore often a potential for significant efficiency savings. The BEM is not as widely applicable as the domain methods; when problems are non-linear, for example, the application of the development of a suitable BEM requires significant further adaption [46][47][48][49]. Typically, the boundary element method has found application in sound reproduction modelling, such as loudspeakers [50,51], sonar transducers [52] and in modelling noise from vehicles [53][54][55][56] and, more recently, aircraft [57][58][59].
For the acoustic boundary element method to be accessible, and hence widely used, it has to be implemented in software. This was precisely the rationale behind the development of the software [60] and the monograph [61], the latter also serving as a manual. Several texts on the same theme were published at about the same time [62,63], following on from the earlier collection of works [64]. A chapter of a recent book contains a modern introduction to the acoustic boundary element method [65].
Implementing the acoustic BEM as software continues to be challenging in terms of scoping, the choice of sub-methods and efficiency. The method requires matrices to be formed with each component being the result of an integration, with some of the integrals being singular or hypersingular. The standard BEM requires the solution of a linear system of equations that is formed from the matrices, or, if the BEM is used for modal analysis, a non-linear eigenvalue problem. There have been significant reliability issues with the BEM for exterior problems. These challenges have been met, but much of the research is focused on future-proofing the method, with the increasing expectations in scaling-up and maintaining reasonable processing time. There is a continual desire to progress to higher resolutions of elements, particularly as this is necessary for modelling high frequency problems.

Acoustic Model
In this Section, the underlying acoustic model of the wave equation that governs the sound pressure in the domain is stated. It is shown how the model can be revised into a sequence of Helmholtz problems for periodic signals. The classes of domains that can be solved by boundary element methods are summarised and a generalised boundary condition is adopted. The other acoustic properties, such as sound intensity, radiation efficiency and sound power that that are mostly used in exterior 'noise' problems, are defined. Traditionally, acoustics properties are presented in the decibel scale, and it is shown how they are converted.

The Wave Equation and the Helmholtz Equation
The acoustic field is assumed to be present in the domain of a homogeneous isotropic fluid and it is modelled by the linear wave equation, where Ψ (p, t) is the scalar time-dependent velocity potential related to the time-dependent particle velocity V(p, t) by V(p, t) = ∇Ψ (p, t) and c is the propagation velocity (p and t are the spatial and time variables). The time-dependent sound pressure Q(p, t) is given in terms of the velocity potential by Q(p, t) = −ρ ∂ ∂t Ψ (p, t) where ρ is the density of the acoustic medium. The time-dependent velocity potential Ψ (p, t) can be reduced to a sum of components each of the form Ψ (p, t) = Re ϕ(p) e −iωt , where ω is the angular frequency (ω = 2πν, where ν is the frequency in hertz) and ϕ(p) is the (time-independent) velocity potential. The substitution of the above expression into the wave equation reduces it to the Helmholtz (reduced wave) equation: where k 2 = ω 2 c 2 and k is the wavenumber. The complex-valued function ϕ relates the magnitude and phase of the potential field.
Similarly, the components of the particle velocity have the form V(p, t) = Re ∇ϕ(p) e −iωt . Often the boundary normal velocity v(p) is given as a condition or is required and this is defined as follows, v(p) = ∇ϕ(p)·n p = ∂ϕ(p) ∂n p , where n p is the unit normal to the boundary at p.

Acoustic Properties
To carry out a complete solution, the wave equation is written as a series of Helmholtz problems, through expressing the boundary condition as a Fourier series with components of the form in Equation (2). For each wavenumber and its associated boundary and other conditions, the Helmholtz equation is then solved. The time-dependent velocity potential Ψ (p, t) can then be constituted from the separate solutions, but it is more usual that the results are considered in the frequency domain.
The sound pressure p(p) at the point p in the acoustic domain is one of the most useful acoustic properties, and it is related to the velocity potential by the formula p(p) = iωρϕ(p). In practice, the magnitude of the sound pressure is measured on the decibel scale in which it is evaluated as the sound pressure level as 20 log 10 p(p) √ 2p * , where p * is the reference pressure of 2 × 10 −5 Pa. Particularly for 'noise' problems, the sound power, the time-averaged sound intensity and radiation efficiency are often considered to be useful properties. The normal sound intensity I(p) at points p on a boundary is defined by the formula where p represents the complex conjugate of p. The sound power W is an aggregation of the sound intensity to one value by direct integration, where H is a boundary. The sound power is also often expressed in decibels as the sound power level as 10 log 10 W W * , where W * is the reference sound power of 10 −12 watts. The radiation ratio is defined as

The Scope of the Boundary Element Method in the Solution of Acoustic/Helmholtz Problems
In applying the boundary element method to acoustic or Helmholtz problems, the user is in effect adopting a model for the boundary/ies and domain(s), the nature of which determines the integral equation(s) that is/are employed. Traditionally, the BEM has been developed to solve the acoustic or Helmholtz problem interior or exterior to a closed boundary in the standard physical domains are two-dimensional, three-dimensional and axisymmetric three-dimensional space. The basis of the method is the integral equations that arise through applying Green's theorems (direct BEM) or by using layer potentials (indirect BEM). The exterior problem has received much more attention, because of its value in solving over an infinite domain from a surface mesh and because of the difficulty in resolving its reliability problem. Domain methods, such as the finite element method, may also be applied; more elements are required, but the matrices are sparse and structured and the FEM is more established than the BEM in general. If the FEM is applied to exterior problems then techniques such as infinite elements [66] can be used to complete the outer mesh or the perfectly matched layer may be used to absorb outgoing waves [67][68][69]. The solution of the interior and exterior acoustic problem by the BEM is considered in Sections 4.1 and 4.3.
An outline of the equations that arise in the finite element and related methods is given in Section 5.3.1. In domain methods like the FEM, whether it is applied to a structure or an enclosed fluid, a linear eigenvalue problem is a natural consequence. With the FEM, it is routine to extract the modes and resonant frequencies and use the modal basis to determine the response under excitation. From that perspective, it is natural to include the eigenvalue problem within the boundary element library. However, the application of the BEM leads to a non-linear eigenvalue problem, which requires special solution techniques and hence the construction of solutions through the modal basis is not a natural pathway in the BEM. Acoustic modal analysis via the BEM is considered in Section 4.2.
Although significantly restrictive, one of the simplest acoustic models is the Rayleigh integral [70]. The model is that of a vibrating plate set in an infinite rigid baffle and radiating into a half-space. The Rayleigh integral relates the velocity potential or sound pressure at any point to the velocity map on the plate. As the Rayleigh integral shares its one operator with the integral equation formulations in the boundary element method, the Rayleigh integral method (RIM) [71] can be adopted into the boundary element method fold. The Rayleigh integral is a special case of half-space problems and these are considered in Section 5.1 An acoustic model that uses the same operators-and hence fits in the BEM context-is that of a shell discontinuity. For example, this can be used to model the acoustic field around a thin screen or shield. The integral equation formulation for the shell are derived from those used in the traditional BEM, but taking the limit as the boundary becomes thinner [72,73], leaving a model that relates a discontinuity in the field across the shell. Shell elements in acoustics are considered in Section 5.2.
The development and application of the BEM from the models described are reviewed in this paper. Although this is a fairly exhaustive list of basic models as things stand, hybrid models can also be developed by using superposition or applying continuity and they will also be considered. The BEM may be applied in vibro-acoustics, aero-acoustics and inverse acoustics and these areas are discussed in Sections 5.3-5.5. Hybrid boundary element methods that include half-space formulation are covered in Section 5.1.2 and shells in Section 5.2.2.
There have been recent developments in the development of the BEM for periodic structures, with noise barriers being the usual application [74][75][76][77][78][79][80][81]. This special case will not be developed further in this paper.

Boundary Conditions
To maintain reasonable generality in the software, the author has generally worked with the boundary condition of the following general (Robin) form with the condition that α(p) and β(p) cannot both be zero at any value of p. This model includes the Dirichlet boundary condition by setting β(p) = 0, the Neumann boundary condition by setting α(p) = 0 and an impedance condition by setting f (p) = 0. In the boundary condition (7), p is any point on the boundary and α, β and f are complex-valued functions that may vary across the boundary. Although the boundary condition for the shell is an adaption of this, this generalised boundary condition model is apparently achievable and seems to cover most current expectations. An explanation of typical boundary conditions that occur in acoustics can be found in Marburg and Anderssohn [82]. For exterior problems, it is also necessary to introduce a condition at infinity, in order to ensure that all scattered and radiated waves are outgoing. This is termed the Sommerfeld radiation condition. In two-dimensional space the condition is and in three dimensions. A thorough discussion on the Sommerfeld radiation condition can be found in Ihlenburg [43] (pp. 6-8).

The Boundary Element Method and the Laplace Equation Stem
Laplace's equation is the special case of the Helmholtz Equation (3) with k = 0, Although Laplace's equation models many phenomena, it is not directly useful in acoustic modelling. However, many of the issues that have to be tackled in solving the Helmholtz equation by the boundary element method are also found with Laplace's equation. The methods applied in developing the BEM for Laplace's equation also can be developed to be used in the Helmholtz problems. As such, applying the BEM to Laplace's equation is a useful entry point in initiating work on the Helmholtz equation. For example, much of the author's work in solving Helmholtz problems has been underpinned by corresponding work on Laplace's equation [83][84][85][86][87]. Thus, in this section, the BEM is communicated in its simplest, but still realistic and practical context as a stepping-stone on the journey to the full scope of the BEM in acoustics, the subject of this paper. The theoretical basis is set out in Kellogg [88] and development of integral equation methods in Jaswon and Symm [89]. In this Section, the boundary element method is derived for Laplace's equation, and this forms the foundation for its application to acoustic problems. In terms of comunication, it is helpful to sustain a notational conformance, that continues through this paper, and is helpful in relating mathematical expressions and their discrete equivalent as software components. Integration methods for finding the matrix components are considered. Interesting and useful properties of some of the operators and matrices are revealed.

Elementary Integral Equation Formulation for the Interior Problem
The boundary element method is not based on the direct solution of the PDE, but rather its reformulation as an integral equation. Historically, the integral equation reformulation of the PDE has followed two distinct routes, termed the direct method and the indirect method. The direct method is based on Green's second identity where ϕ and ψ are twice-differentiable scalar function in a domain D that is bounded by the closed surface S. If ϕ is a solution of Laplace's equation, Green's third identity can be derived from the second identity by choosing ψ(q) = G(p, q), where G is a Green's function. A Green's function is a fundamental solution of the partial differential equation, in this case Laplace's equation, that is the effect or influence a unit source at the point p has at the point q, and is defined by ∇ 2 G(p, q) = −δ(p − q). For the two-dimensional Laplace equation G(p, q) = − 1 2π ln(r) and for three-dimensional problems G(p, q) = 1 4πr where r = p − q . The substitution of the Green's function into Equation (12) gives Hence, as a result of the properties of the Dirac delta function, where c(p) is the angle (in 2D, divided by 2π) or solid angle (in 3D, divided by 4π) subtended by the interior domain at the boundary point p. (Similarly, for exterior problems, c(p) similarly relates to the exterior angle). If the boundary is smooth at p then c(p) = 1 2 , and, for simplicity, this is the value that will be used for the remainder of this paper.
In the indirect method, ϕ is presumed to be related by a layer potential σ on the boundary Differentiating Equation (15) with respect to a unit outward normal to a point on the boundary that that passes through p, gives the following equation: As p approaches the boundary, the integral operator has a 'jump' similar to the direct method (14) resulting in the following equation: The BEM can be derived from the equations in this Section. The equation defined on the boundary (the last one in Equation (14) for the direct method, Equations (15) and (17) for the indirect method (p ∈ S)), the boundary integral equations, are used to find the unknown functions from the known data on the boundary. The corresponding integrals defined in the domain (the first one in Equation (14) for the direct method and Equation (15) for (p ∈ D) in the indirect method) return the solution at the chosen domain points.

Operator Notation and Further Integral Equations for the Laplace Problem
In this Subsection, operator notation is introduced and further integral equations are introduced in order to illuminate some useful properties. Operator notation provides a shorthand and is an aid to communication. The Laplace integral operators are defined as follows: where Γ is a boundary (not necessarily closed), n q is the unique unit normal vector to Γ at q, v p is a unit directional vector passing through p and µ(q) is a function defined for q ∈ Γ. With this notation, the Equations (14), which form the basis of the elementary direct method for the interior problem, can be written as Similarly, the equations for the indirect method ((15) and (17)) in operator notation for the interior problem are as follows: Further integral equations for the direct formulation can be obtained by differentiating (21) (as in the derivation of Equations (16), (17)) and for the indirect formulation by introducing a double layer potential. Although they are generally unnecessary in solving Laplace problems, the equations resulting from differentiating the elementary integral equations or using double layer potentials are often used in the exterior acoustic problem, considered in the next Section. These equations are also useful in general and they will illuminate some important aspects. For example, differentiating Equation (21) with respect to the outward normal to the boundary, which can be written in operator notation, in which a new operator, N, has been introduced, Further indirect integral equations may be obtained through presuming the field is the result of a double layer potential and, taking into account the jump discontinuity, The integral equations for the exterior Laplace problem may be derived in the same way. In general, the equations are the same as for the interior problem, except for a change of sign. The equations that make up the direct formulation are The equivalent of equations for the indirect method ( (15) and (17)) in operator notation for the exterior problem are similar, with just a couple of sign changes to indicate that the jump discontinuity in M and M t in the limit from the exterior, rather than the interior:

Derivation of the Boundary Element Method
There are several integral equation methods that can be applied to transform integral equations into boundary element methods, but the most popular and straightforward method is that of collocation [90]. The development of the boundary element method from the selected integral equation requires that the boundary is represented by a set of panels or a mesh. An integral equation method is applied to solve the boundary integral equation. The solution in the domain can then be achieved by effecting the appropriate integration over the boundary.

Boundary Element Approximation
In the spirit of the author's previous work [61,91], in order to maintain generality, the discrete forms of the Laplace operators (18)- (20), (26) are sought, in order to effectively become a software component. Let the boundary Γ in Equation (18) be represented by the approximationΓ, a set of n panels: The boundary function µ is replaced by its equivalentμ on the approximate boundaryΓ: In general, the function on the boundary is replaced by a sum of a set of basis functions. The simplest approximation is that of approximating the boundary functions by a constant on each panel: whereẽ = 1. For example, for the simple boundary integral Equation (22) with p ∈ S, The most common method of solving boundary integral equations is collocation, in which a linear system is formed through setting p to take the value of each collocation point in turn:

Solution by Collocation
The linear system of approximations (41) may be written in matrix-vector form where ϕ where ϕ Hence approximations to the solution within the domain ϕ D may be found by the matrix-vector multiplication

The Galerkin Method
Although most of the implementations of the boundary element method are derived through collocation, other techniques can be used, most typically the Galerkin method. The Galerkin method and collocation can both be derived from a more generalised approach that are termed weighted residual methods. For example, for the approximation (40), the residual is the difference between the approximation and the exact solution; In weighted residual methods, R σ S ; p is integrated with test basis functions χ i (p) (p·S) and the methods arise by setting the result to zero; S R σ S ; p i χ i p i dS = 0 at points p i on the boundary, for i = 1, 2, . . . n S . If χ i (p) = δ p − p i , the Dirac delta function, the collocation method is derived; which leads to the methods outlined above. In the Galerkin method, the test functions are the same as the original basis functions. For example, for constant elements, the basis and the test functions are defined as Substituting this and the definition of the residual equation The reason for the relative unpopularity of the Galerkin approach in boundary elements is that the matrix is now the result of a double integration, rather than the single integration in the collocation method. However, the matrix is understood to be symmetric.

Properties and Further Details
In the boundary element method, the boundary is represented by a mesh of panels. The boundary functions are approximated by a linear combination of basis functions on each panel. The boundary element is the combination the panel and the functional representation. In the previous Subsection, constant elements were mentioned as an example. In the finite element method, the degree of the approximation has to be at least half the order of the PDE. Fortunately, this is not the case for the application of the boundary element method, where constant elements are widely used. The panels that make up the boundaries are most simply represented by straight line panels in 2D, triangles in 3D and conic sections for axisymmetric 3D problems.
In this Subsection, important properties and further details of Laplace's equation, the related boundary integral equations and operators are discussed. An overview of methods for carrying out the integrations is included. The issue of non-uniquesness, that is an important feature of the boundary integral equation formulations for the exterior acoustic problem, is first addressed with Laplace's equation and useful outcomes from this are outlined.

Integration
In the previous Subsection, it was shown that up to four operators are involved in the boundary element method for Laplace problems, and these four operators extend to acoustic/Helmholtz problems. For the Laplace problem, it may be possible to develop analytic expressions for the integrals [92], but in general-and particularly for acoustic/Helmholtz problems-numerical integration is necessary. In the main, the integrals are regular, and these are most efficiently evaluated by Gauss-Legendre quadrature [93]. This is straightforward to apply to straight-line panels and in the generator and (typically in composite form) azimuthally. There are also published points and weights for Gaussian quadrature on a triangle [94].
However, the integrals corresponding to the diagonal components of the L SS matrix are weakly singular and the diagonal components of the N SS matrix are hypersingular, that is when the point lies on the panel. Moreover, if the point is close to the panel, for example at the centre of a neighbouring panel then the integrals are said to be nearly singular and may require a more accurate numerical integration rule, or special treatment [95]. Special numerical integration methods can be applied in order to evaluate the singular integrals and expressions for the hypersingular integrals may be found through limiting process. The weakly singular integrals, corresponding to the diagonal components of the L SS matrix have a O(ln r) singularity in 2D and an O 1 r singularity in 3D, where r is the distance from the central collocation point. For the axisymmetric 3D case, the azimuthal integration resolves the O 1 r singularity to an O(ln r) singularity on the generator. For the simple elements discussed, analytic expressions for the diagonal components of the L SS and N SS matrices are available for the 2D and 3D (non-axisymmetic) cases for Laplace's equation. Further work on singular integration can be found in the following references [96][97][98][99][100][101].
Although for simplicity, the four operators are often lumped together as integral operators, N is not an integral operator: N is termed a pseudo-differential operator, and it therefore has distinct properties. The full expressions for the straight line and triangular panels are given in the author's book [61] (pages [49][50]. For illustration the expressions are given for the straight-line panel of length h (2D) and for an equilateral triangular panel (3D) with each side of length h: for the equilateral triangular panel with side of length h (3D) (46) There is one simple but noteworthy remark to be made about the expressions (45) and (46). Normally, if the domain of integration is reduced in size, the integral is similarly reduced, at least in the limit as the domain size converges to zero. However, with these hypersingular integrals, the opposite is found to be the case! A significant consequence of this property is considered in Section 6.2.

Non-Uniqueness and Its Useful Outcomes
The non-uniqueness of some boundary integral equations in the exterior acoustic problem at a set of frequencies is well-known, and the issue and methods of resolution will be considered in the next Section. Because of the importance of the non-uniqueness of solution within the context of this paper, it is helpful to visit this early, with a simpler equation. Although the term non-uniqueness often prefixes the word problem, it is found that some very useful outcomes arise from this analysis.
The non-uniqueness is found with Laplace's equation itself; if ϕ is a solution of the interior Neumann problem then ϕ + c is also solution, where c is any constant. The non-uniqueness in the interior Laplace problem with a Neumann boundary condition must be reflected in the boundary integral equations. It therefore follows from Equation (21) that the operator M + 1 2 I is degenerate, and similarly for M t + 1 2 I from Equation (23) and N from Equations (25) and (28). As we will see in the next Section on acoustic/Helmholtz problems, the operators from the interior problem equations are shared with the exterior problem, and this is shown for Laplace's equation in Section 3.2. For the exterior problem, the derivative direct formulation (31) is unsuitable for solving exterior Laplace problems as both the N and M t + 1 2 I operators are degenerate. Similarly, the indirect formulations, formed from double layer potentials, (35) and (36) are unsuitable. This correlation between the interior Neumann problem and the exterior derivative (direct) and double-layer potential (indirect) formulations carries through to the Helmholtz equation.
Most simply, if v(p) = 0 on the boundary then ϕ is any constant (e.g., ϕ = 1) throughout the domain is a solution. Substituting these values into the discrete form of the boundary integral Equation (21), the resulting matrix-vector equation is M SS + 1 2 I 1 ≈ 0, where 0 is a vector of zeros and 1 is a vector of ones; every row of the M SS matrix must approximately sum to − 1 2 . Similarly, from Equation (25), N SS 1 ≈ 0; every row of the N SS matrix must approximately sum to zero. A potentially useful outcome of this is that the hypersingular diagonal components of the N SS matrix can be computed from the others, which are all 'regular' (at least for simple elements). This can be taken a step further and, the panel in question may be linked to a fictitious boundary made up of panels and the value of the hypersingular integral determined by summing the other integrals. The values of the singular and hypersingular integrals in the BEM solution of Laplace's equation may be stored and be used to subtract out the singular and hypersingular components of the same integrals for the Helmholtz equation. The method of inventing a fictitious surface to evaluate the hypersingular integrals is precisely the method used in the author's axisymmetric programs [60], as there were analytic expressions for the hypersingular integrals for the panels used in 2D and general 3D, as discussed in the previous Sub-subection, but no other similar way forward for the axisymmetric panels.
The results in the previous paragraph also provide useful methods of validation. One row of the M SS and the N SS matrices need to be evaluated that that is when the point p is on the boundary. The result of summing the rows of M SS can verify that a 'closed' boundary is actually closed, if their values are − 1 2 and zero, and, if they are not then it indicated that the boundary may be open or there are errors in the mesh. Finally, substituting these values into Equation (30) gives the following, where e is the unit function; useful in validating that the solution points are within the domain.
These simple validation methods, arising from potential theory, are applicable to all BEMs and beyond; any computer simulation involving surface meshes may benefit from these simple techniques. The author includes these methods of validation routinely in his boundary element codes.

The Standard Boundary Element Method in Acoustics
In the author's book and software [60,61], three classes of acoustic problem were considered; the determination of the acoustic field within a closed boundary, outside of a closed boundary and the interior modal analysis problem. In this Section, the three methods are revisited and recent developments and applications are included. Recently, the software has been re-written in Python [102].

The Interior Acoustic Problem
In this Subsection, the BEM is developed for the solution of acoustic problems in a domain that is interior to a closed boundary [61,103]. The method has been applied to room acoustics [104][105][106][107][108], the interior of a vehicle [109][110][111], modelling sound in the human lung [112,113] and in biological cells [114].

Integral Equation Formulation
The direct integral equations for the interior acoustic problem, reformulating the Helmholtz Equation (10), but follows the same format as the formulation for Laplace's Equation (21); Equation (48) introduces two Helmholtz integral operators that are analogous to the L and M operators for Laplace's equation, and are defined similarly; where the Green's function is defined as follows: , for three-dimensional problems.
As with the interior Laplace Equations (22) and (23), the indirect integral equation is derived through presuming that the field is generated by a layer potential, defined on the boundary; where, similarly with Equation (20), the operator M t k is defined as follows:

The Boundary Element Method for the Generalised Boundary Condition
Substituting the expressions (53) and (54) into the boundary condition (7) gives Following the derivation of the BEM through collocation, this resolves to a linear system of the form where D α and D β are diagonal matrices, with the values of α and β aligned along the diagonal. Equation (56) is solved in the primary stage of the boundary element method, yielding an approximation σ S to σ S . In the secondary stage, the discrete equivalent of Equation (53) returns the solution at the interior points: ϕ For the direct formulation, the discrete equivalent of Equation (48) which is solved with the discrete equivalent of the boundary condition (7), Comparing the linear systems for the indirect method (56) with that of the direct method (59) and (60) illustrates an apparent significant advantage in the indirect approach, the number of unknowns in the system corresponding to the direct method is twice that of the indirect method. However, through exchanging columns, the system for the direct method can be reduced to an n × n system, and this method has been automated [115].
There have been several papers discussing the accuracy of the interior acoustic boundary element method near to the resonance frequencies. It is considered that real eigenvaulues are shifted into the complex plane following discretisation. This phenomenon is termed numerical damping [116][117][118][119].

Equations of the First and Second Kind
Integral equations with a fixed region of integration are termed Fredholm integral equations. Fredholm integral equations are categorised as first kind or second kind. Equation (53) is a typical first kind equation, in which we are solving over integral operator(s) alone, in this case the L k operator, to find σ from ϕ. With second kind equations, we are solving not just over integral operators, but also the identity (or diagonal) operator. For example, Equation (54) is a second kind equation, solving over the M t k + 1 2 I operator, in order to find σ from v. In general, first kind equations have poor numerical properties and are avoided. Although pure first kind equations only occur in particular Dirichlet cases, they can be avoided, through using the derivatives of the integral equation formulations (direct) or double layer potentials (indirect). However, from experience, first kind equations with singular kernels, as we find with the L or L k operator, do not have poor convergence properties. It follows, therefore, that the formulations and methods developed thus far in this Subsection are generally suitable in solving the interior acoustic problem.

Derivative and Double-Layer Potential Integral Equations for the Interior Helmholtz Problem
As with the equations listed in Section 4.1.1, the derivative equations are unnecessary in solving the interior Helmholtz problem. However, the derivative equations provide an alternative formulation and help with the general understanding. The form of the equations is analogous with the equations for the interior Laplace equation, developed in Section 3.2.
For the direct method, the derivative boundary integral equation is analogous to Equation (25): where For the indirect method, the boundary integral equations follow the form of Equations (28) and (29):

Interior Acoustic Modal Analysis: The Helmholtz Eigenvalue Problem
An enclosed acoustic domain has resonance frequencies and associated mode shapes. Mathematically, these are the solutions k * (the eigenvalues-that relate to the resonance frequencies) and ϕ * (the eigenfunctions-that relate to the mode shapes) of the Helmholtz equation with the homogeneous form of the boundary condition (7) (i.e., with f (p) = 0). This problem is more typically solved by the finite element method and a finite element model of an acoustic or structural problem is developed in Section 5.3.1. In this Subsection the modal analysis of an enclosed fluid via the boundary element method is outlined.

The Generalised Non-Linear Eigenvalue Problem from the Boundary Element Method
In the indirect boundary element method, the eigenvalues k * and the eigenfunctions σ * are found through solving that follows from Equation (56), with the true eigenfunctions ϕ * can then be found with Equation (53). Through applying collocation, this is equivalent to solving the non-linear algebraic eigenvalue problem which follows from Equation (57).
The eigen-solution of Equation (66) returns the approximationsk * to the eigenvalues and the approximation σ S to the layer potential eigenfunction. The approximations to the physical eigenfunctions at the chosen domain points can then be found with Equation (58). The method described can also be applied with the direct integral equation formulation but requiring the row-exchanging method mentioned in Section 4.1.2.
Although in this work, the focus is on the generalised boundary condition, most of the examples in the literature consider the Dirichlet and Neumann eigenvalues, and it is revealing to focus on those. The Dirichlet/Neumann interior eigenvalues and eigenfunctions are the solutions of the equations of this Subsection with ϕ(p) = 0 / v(p) = 0 (p ∈ S). Hence the Dirichlet eigenvalues are the eigenvalues of the L k operator from Equations (48) and (53), of the M k − 1 2 I operator from Equation (63) and the M t k − 1 2 I operator from Equation (61). Similarly, the Neumann eigenvalues are the eigenvalues of the M k + 1 2 I operator from Equation (48), M t k + 1 2 I from Equation (54) and the N k operator from Equations (61) and (64).

Solving the Non-linear Eigenvalue Problem
The standard algebraic eigenvalue problem has the form where A is a square matrix. Because of the analogy between the Helmholtz Equation (3) and the standard algebraic eigenvalue problem (67), the Helmholtz eigenvalues are sometimes termed the eigenvalues of the Laplacian [120][121][122][123][124].
The generalised algebraic eigenvalue problem has the form where A and B are square matrices. Standard methods exist for solving these problems, termed the QR and QZ algorithm. As we have noted, the application of the boundary element method in Equation (66) results in a non-linear eigenvalue problem, and this has the form Although this is a significantly complicated problem, at least-in the case of the boundary element matrices-the individual components of the components matrix are continuous. Earlier methods tended to find the eigenvalues by finding the zeros of |A k |. Further developments and applications can be found in the following research papers [61,.

The Exterior Acoustic Problem
The boundary element solution of the exterior acoustic problem is the most popular area of research in the context of this paper. A method that can solve over a theoretically infinite domain from data on a surface mesh has significant value in the context of acoustics. However, it was found more than half a century ago that the boundary integral equations, derived in the same way as the ones for the interior problem, resulted in unreliable boundary element methods. Much has been achieved on the road to developing a more successful outcome. However, to obtain a reliable exterior acoustic boundary element method remains a tantalising goal.

The Integral Equation Formulations of the Exterior Helmholtz Equation and Their Properties
The direct boundary element reformulation of the Helmholtz equation follows the format for Laplace's Equation (30) and is as follows: The indirect formulation is However, the L k , M k − 1 2 I and M t k − 1 2 I operators are degenerate at the eigenvalues of the interior Dirichlet problem (see Section 4.2.1). The issues in using the boundary integral equations in (70)-(72) as general purpose exterior acoustic/Helmholtz equation solvers has been known for over 50 years. The wavenumbers or frequencies in which these boundary integral equations are unsuitable are often termed the characteristic or irregular wavenumbers or frequencies in the literature. These characteristic wavenumbers are physical in the interior problem, but they are unphysical in the exterior problem in which they are not a property of the Helmholtz equation model but are manifest in the boundary integral equations. In Section 3.4.2 the effects of the non-uniqueness of the solution of the interior Laplace problem with a Neumann boundary condition were discussed. The boundary integral Equations (70)-(72) have a similar property at the characteristic wavenumbers, and this is also often termed the non-uniqueness problem.
Although it is 'unlikely' in practice that the value of wavenumber k is equal to a characteristic wavenumber k * , it is shown in Amini and Kirkup [230] that the numerical error as a consequence of being 'close' to a characteristic wavenumber is O( 1 |k−k * | ). Hence, one technique of resolving the problem is to use finer meshes in the neighbourhood of the characteristic wavenumber in order to offset this error [231]. However, this strategy is likely to be prohibitive, requiring the overhead of creating a range of meshes and increased computational cost. The values of the characteristic wavenumbers are also generally unknown, although they can be found (as discussed in Section 4.2), but at a substantial computational cost. It is also found that the characteristic wavenumbers cluster more and more as the frequency increases. In conclusion, therefore, the boundary integral equations are-in practice-only useful for frequencies reasonably below a conservatively estimated first characteristic wavenumber, severely restricting the methods to the low frequency range in practice.

The Derivative and Double Layer Potential Integral Equations
For the direct method, the derivative boundary integral equation is analogous to Equation (31): For the indirect method, the boundary integral equations follow the form of Equations (35) and (36):

The Schenck Method
The Schenck method [149] is often termed the CHIEF method in the literature and it is a development of the standard direct method, based on Equation (70). Given the potential non-uniqueness, discussed in Section 4.3.1, the system of equations that form the discrete equivalent of the boundary integral equations are regarded as potentially underdetermined and they are augmented with equations related to a set of points in the interior D: By adding more equations, the expectation is to eliminate the non-uniqueness and determine a solution.
The equations can be solved by the least-squares method for Dirichlet and Neumann problems, the column exchanging algorithm [115] could be used in the case of the general boundary condition (7). There are several reported implementations and applications of the Schenck method [216,[232][233][234].
Obviously, there are immediate questions about the number and position of the interior CHIEF points. Juhl [235] develops an iterative method for selecting points and halting when the results are unchanged. Equation (47) can verify that CHIEF points are interior, so this could be usefully included in the method. There has been a significant number of implementations and testing of the Schenck method. In general, more and more points are required to offset the non-uniqueness, as k increases. This increases the computational overhead with respect to the wavenumber, additional to the potential need for more elements at higher wavenumbers. Several improvements in the method have been put forward and these are summarised in Marburg and Wu [236]. There have been several evaluations of the CHIEF and combined methods and their variants [157,237,238].

The Combined Integral Equation Method
In this method, a boundary integral equation is formed through a linear combination of the initial boundary integral equation and its corresponding derivative equation (for the direct method and double-layer potential for the indirect metod). This concept was initially derived for the indirect method and is attributed to Brakhage and Werner [239], Leis [240], Panich [241] and Kussmaul [242]. The corresponding direct integral equations are attributed to Burton and Miller [243].
The Burton and Miller method is based on a boundary integral equation that is a linear combination of the initial one (70) with the derivative (73), where µ is a complex number, a coupling or weighting parameter. Similarly, for the indirect method, the equation is based on writing ϕ as a linear combination of a single and double layer potential This returns the following boundary integral equations The issue of the determination of the values for the coupling parameter will be revisited in Section 6.2.

Scattering
The boundary integral equations for the exterior acoustic problem are readily applicable to radiation problems. The equations can also be used for the scattering problem, but with extra terms involved that model the incident field. These same techniques can be applied in the interior problem. For example, for the indirect method, the exterior acoustic field (77) is presumed to be made up of the layer potentials, superposed with the (known) incident field: This similarly adjusts Equations (78) and (79):

Extending the Boundary Element Method in Acoustics
In this Section, extensions in the standard boundary element methods of the previous section are considered. These include the Rayleigh integral method for computing the acoustic field surrounding a vibrating plate set in an infinite reflecting baffle, as a case of the more general half-space methods. This Section also includes shell elements in which a revision in the boundary integral equation for the exterior problem returns a model for the acoustic field surrounding a thin screen. Through principles of continuity and superposition, hybrids of these models and the standard models also significantly extend the range of acoustic problems that come under the boundary element fold. Vibro-acoustic, aero-acoustic and inverse acoustic problems are also considered in this Section.

Half-Space Methods
In Section 3 it was stated that the boundary element solution of Laplace's equation was a useful entry to the BEM in acoustics. The Rayleigh integral method (RIM) is also a good a starting point, in that it required only one of the four Helmholtz operators, and, for Neumann problems, it is an integral, rather than an integral equation. In this Subsection, the Rayleigh integral method is defined and further developments are reviewed.
The Rayleigh integral method can be viewed as a solution to a half-space problem. If further boundaries are placed in the half-space, then the integral equation formulation, with a simple modification of the Green's function, forms the model with ϕ = 0 or ∂ϕ ∂n = 0 on the plane. Further development of the Green's function have been researched in order that an impedance condition is modelled on the plane, a useful model in outdoor sound propagation.
On the other hand, if there is a cavity in the plane then the interior boundary element method can model the acoustic field within the cavity and this is coupled to the Rayleigh integral. The model is based on based on the interior formulation to model the cavity, applying a false flat boundary across the opening and coupling the interior formulation with that of the half-space. This method is a hybrid of the boundary element method and the Rayleigh integral method and is termed BERIM.

The Rayleigh Integral Method
In the operator notation of this paper, the Rayleigh integral is as follows: where Π is the flat plate, radiating into the half-space E + . The solution of the Neumann problem, finding ϕ from v, is simply the evaluation of an integral. For the general boundary condition of the form (7), the technique is again to solve the boundary integral Equation ((80) with p ∈ Π) to obtain v and then evaluating the integral to obtain ϕ at any point in E + . In the Rayleigh integral method [71], the plate is divided into elements, as discussed and applying collocation or, the equivalent for the integral, product integration. Substituting (80) into the boundary condition (7) gives Using the discrete notation of this paper, the equation following the application of collocation is as follows: There have been several reported implementations and developments of the Rayleigh integral method [244][245][246][247][248][249][250][251][252][253][254], including vibro-acoustics [197,255,256]. Applications of the method include sandwich panels [257][258][259][260][261][262][263], engine or machine noise [203], electrostatic speaker [264] and transducers [52,265,266].

Developments on Half-Space Problems
Integral equations for scattering or radiating boundaries above an infinite plane can be developed through altering the Green's function [247] in order that it also satisfies the condition on the plane. For the perfectly reflecting plane the Green's function must satisfy the Neumann condition on the plane; ∂G * ∂n = 0 and hence G * (p, q) = G(p, q) + G(p, q * ), where q * is the point that corresponds to q, when reflected through the plane. Similarly, if there is a homogeneous Dirichlet condition on the plane then the revised Green's function is G * (p, q) = G(p, q) − G(p, q * ). More generally, the modified Green's function is G * (p, q) = G(p, q) + R G(p, q * ), with R representing the reflection of the plane (−1 ≤ R ≤ 1). An impedance boundary condition is generally required in modelling environmental noise problems [220,[267][268][269][270].
However, the methods in the previous paragraph are applicable when the acoustic scatterers or radiators are on or above the plane. Another set of methods apply if there is a cavity in the plane. Early examples of this type of problem have arisen in harbour modelling, in which the Helmholtz equation has been used to model the waves in a harbour (the cavity), which is open to the sea bounded by the coastline (the plane) [271,272]. The Boundary Element-Rayleigh Integral Method (BERIM) [169], is applicable to open cavity problems in acoustics. The interior boundary element method (Section 4.1) models the sealed interior and the Rayleigh integral method models the field exterior. The advantage in this model over the exterior model could be substantial; the mesh covers the interior of the cavity and the opening only, rather than the inside and outside. There is an important issue with the model, as it presumes that the cavity opens out on to an infinite reflecting baffle. However, this could be a small price to pay, and some problems-such as the loudspeaker problems considered by the author [169]-the mouth opens onto a front face of the cabinet. Motivated again by environmental noise problems, several methods have been developed for a cavity opening on to an impedance plane [273].

Shell Elements
The derivation of integral equations and methods for modelling thin shells (that is an open boundary modelling a discontinuity in the field) in the boundary element context takes us back to the works such as Krishnasamy [274], Gray [275], Terai [276] and Martin [277]. In these, and various other references, the shell is also termed a crack. For the Helmholtz equation, the derivation of the integral equations is set out in Warham [73] and Wu [72].
Shell problems may be solved using the standard boundary element methods already described in this paper. One method is to mesh the shell as a closed boundary, with a finite thickness. However, if the thickness of the shell is a fraction of the element size, then the boundary integral equations representing collocation points on either side of the shell are similar, and the equation approaches degeneracy. Alternatively, many elements may be required, to ensure that the elements on the shell are not disproportionate in comparison with those along the edges.
A more productive method of using existing methods is to artificially extend the shell to form a complete boundary, with the interior and exterior BEM applied to the inner and outer domains, and continuity applied over the artificial boundary. However, this approach, in many cases, would be prohibitive, as the number of elements would be substantially greater. However, analytic test problems for shell problems are difficult to develop and using the more-established interior and exterior BEM in this way can provide comparative solutions. As with the Rayleigh integral method of the previous Subsection, the same Helmholtz operators re-occur with the shell model. Hence, the inclusion of shells in boundary element software significantly extends the functionality of the library, at little extra cost.

Derivation of the BEM for Shells
In this model, ϕ is the solution of the exterior Helmholtz Equation (3) in the exterior domain, surrounding an open boundary Ω. The boundary Ω is presumed to be an infinitesimally thin discontinuity and so, at the points on the boundary, ϕ and its normal derivative v take two values, one at each side of the discontinuity. The two sides of the shell are labelled "+" and "-": it doesn't matter which way round this is, as long as it remains consistent. Hence, on the shell, four functions are modelled, ϕ + (p), ϕ − (p), v + (p) and v − (p) (p ∈ Ω), where the normal to the boundary, which orientates v + and v − , is taken to point outward from the '+' side of the shell. However, rather than working with these functions, it is more straightforward to work with the difference and average functions; for (p ∈ Ω) and where, for simplicity, it is presumed that the boundary is smooth. The boundary condition is defined in a similar way as in Equation (7), but as there are double the number of unknown function, two boundary conditions are required (p ∈ Ω): The integral equations that govern the field around the shell discontinuities derived from the exterior direct formulation (70) [73], by taking the limit as the boundary becomes thinner, and they are as follows: There are few tests of methods based on these equations in the literature. The only known issue is at the edge, as it is likely that the solution is singular there. Therefore, mesh grading, using smaller and smaller elements close to the edge may improve efficiency.

Mixing Opem with Closed Boundaries
Again, using the superposition principle, discussed in Section 4.3.5, shell boundaries can be mixed with the traditional boundaries. For example, for the interior problem made up of a domain D with boundary S and with shell discontinuities Γ within the domain, the superposition of the direct Equations (48) and (61) with Equations (85)-(87) returns the following equations: Methods based on this analysis and equations were developed and demonstrated by the author for the interior Laplace equation [85], the exterior acoustic/Helmholtz problem [278], and for the interior acoustic/Helmholtz problem [135,279].

Vibro-Acoustics
Problems that involve structural vibration, as well as an acoustic response, are termed vibro-acoustics. In this Subsection, the modelling by a domain method, such as the finite element method is outlined. The FEM can be applied to the structure and/or the acoustic/fluid domain, but, in the context of this Subsection, the FEM is used as the structural model. When the structure and the fluid significantly influence each other's dynamic response then they are modelled as coupled fluid-structure interaction.

Discrete Structural or Acoustic (Finite Element) Model
The properties of the interior acoustic problem parallel the expected response an excited elastic structure, presuming no damping, or, more simply, simple harmonic motion. In this discussion, let us consider this further in order to bring context. With a mass M and a stiffness K, the equation of (unforced) simple harmonic motion is M ..
where q is the displacement and ..
q is the acceleration. The same equation results from a system of masses or from the finite element method solution with M and K termed the mass and stiffness matrices and q and .. q are vectors of displacement and acceleration of the individual masses or nodes in the FEM.
The phasor solution q = Qe jωt returns the following equation, Hence, applying the finite element method to the structural or interior acoustic modal analysis problem returns a generalised algebraic eigenvalue problem (of the form of Equation (68)). The matrices are sparse and hence are amenable to more efficient methods of solution. Although the matrices are larger and the domain needs to be meshed, the BEM with its nonlinear eigenvalue problem struggles to compete with the FEM in acoustic modal analysis. Let ω 2 j and Q j for j = 1, 2, . . . be the eigenvalues (natural frequencies) and corresponding eigenvectors (mode shapes) of Equation (88). It follows that Let us also now generalise the model in order to include and external excitation force g: Following the phasor substitution g = Ge jωt Equation (89) is modified as follows: For the homogeneous equations (G = 0), the above may be cast as a generalised or standard algebraic eigenvalue problems, as discussed in Section 4.2.2. Let us write the response and the excitation in terms of the modal basis Q = j γ j Q j and M −1 G = j a j Q j .
Considering each eigen-solution in turn relates the coefficients, so that In practice, a structure or an enclosed fluid experience damping, the simplest and usual model is to relate damping to the velocity: M ..
where C is termed the damping matrix. Following the phasor substitutions, Equation (89) is modified as follows: The response of the system to an applied boundary condition across a frequency sweep is to have a smoothed peak at the resonance frequency and a more gradual phase change. In general, the response of the system can be modelled as a sum of weighted modes (90) with the coefficients that are relatable as follows:

Coupled Fluid-Structure Interaction
The finite element model in the previous Sub-subsection is directly applicable to a structure when there is no significant coupling between the structure and a fluid. Similarly, the acoustic analysis methods of Sections 4 and 5.1 and Section 5.2 are directly applicable when there is no significant coupling between the fluid and a structure. However, for many dynamical systems, it is appropriate to couple the structural model, of the FEM form outlined in the previous Sub-subsection, with the acoustic model of the fluid with which it is in contact. The finite element method can be applied in both domains, with appropriate properties. In the context of this work, the boundary element method provides the computational acoustic model. The models are coupled together, through continuity of the particle velocity at the interface, and the resultant forcing on the structure is affected by the sound pressure. The discrete coupling is applied at the elements that coincide with the boundary.
In Section 4.2, the acoustic modal analysis of an enclosed fluid was discussed. Similarly, in the previous Section, structural modal analysis via the finite element method was outlined. Of course, when coupling occurs, the eigenvalue analysis need to be applied to the coupled system. For example, a structure typically exhibits different resonant frequencies in vacuo than it does when immersed in a fluid. In the literature, these are often termed the dry and wet modes (at least for structures placed outside of and in water). Modal analysis via the finite element method returns a standard and sparse eigenvalue problem (88), modal analysis by the boundary element method yields a non-linear eigenvalue problem (69) and hence the hybrid coupled FEM-BEM system is also non-linear. There are several reported implementations of coupled fluid-structure modal analysis using the boundary element method [19,[255][256][257]265,284,287,[289][290][291][292][293][294][295][296][297][298][299][300][301]. The coupled matrices are generally much larger than for the boundary element method alone and this can be resolved by determining the response in the dry structural modal basis and coupling that with the boundary element model [299,301].

Aero-Acoustics
While the boundary element method has been applied to coupled fluid-structure interaction problems for almost as long as the method has been around, the BEM in aero-acoustics is a much more recent development. The standard exterior BEM in a vibro-acoustic settting has been applied in aircraft noise [210,211]. The wave Equation (1) does not include convective flow and is only an adequate model in aero-acoustics in the special case of insignificant flow and the formulations and methods outlined in this work are no longer directly applicable. The Navier-Stokes equations model fluid flow and their solution by numerical methods is termed computational fluid dynamics (CFD). There are reports of applying the standard acoustic BEM and CFD to aircraft [47,48,[302][303][304] and, similarly, to underwater vehicles [305,306].
To model the noise from aircraft, the domains are typically large and significant resolution is required to capture the higher frequencies, and hence domain methods come with a high computational cost. Computational aero-acoustics [11] has arisen for developing and applying numerical methods in this particular area. The attraction of the BEM in computational aero-acoustics is the same as it is in standard acoustics, a significant reduction in meshing and hence the potential for faster computation.
Work on the adaption of the BEM to a wider scope of problems has been developed, for example by the dual and multiple reciprocity method, which has also been applied to variants of the Helmholtz equation [46][47][48]103,131,148]. A generalisation of the boundary element method in acoustics that includes convection, is applicable to problems with uniform flow, but this can also form a useful approximation method with the mean flow rate substituting the value for the uniform flow rate [307]. Similar to the approach in half-pace problems, discussed in Section 5.1.2, aero-acoustic problems are adapted for the BEM by revising the Green's function [308][309][310]. Recently the BEM in acoustics has been adapted to model viscous and thermal losses [311][312][313].
Methods in aero-acoustics that use the BEM generally involve setting a fictitious surface, enclosing the significant effects such as noise generation and turbulence, within which typical methods of CFD are used to model the Navier-Stokes equations; the sound generation and sound propagation are modelled separately. The boundary element method models the outer domain, but requiring the fictitious surface to be meshed, modelling the uniform flow and the radiation condition from the boundary and into the far-field. There are several reported implementations and aero-acoustic and related applications of these methods [307,[314][315][316][317][318][319][320][321][322][323][324][325][326].

Inverse Problems
An inverse problem in acoustics, involves measuring properties of the acoustic field and processing that information in order to determine something about its origin. Over recent decades, the main text that guides inverse acoustic (and electromagnetic) (scattering) problems is that of Colton and Kress [16], now in its third edition. Colton and Kress provide a mathematical analysis of inverse problems, and, in acoustics, their focus is on determining the shape of the boundary of a scatterer from the (disturbed) sound pressure at points in the far-field. An example of an application of this is identifying bodies on the sea floor [327].
Inverse problems are ill-posed. This means a range of solutions can give rise to the same measurements, in contrast to the forward problems, studied so far, which are usually well-posed, with a unique solution. In practice, the linear system that is returned by the boundary element method (or any standard method) is significantly ill-conditioned. If a conventional method is used to solve the linear system then the solution will not be acceptable. In the author's previous work on the inverse diffusion problem (classically, the backward heat conduction problem) [328], also included in Visser [237]) it was discussed that there is effectively insufficient information in the data to determine the solution of an inverse problem. The notion of an 'acceptable' solution, the bias of the observer, provides the final constraint that enables the ill-posed problem to be substituted by a nearby well-posed problem. This returns an acceptable solution, even though it is a less accurate solution of the discrete equations. In the literature, this technique is termed regularisation.
Acoustic holography-determining the sound field near the source from acoustic properties at a distance from the source-is one of the main application areas of the inverse boundary element method. In general, an array of pressure or velocity transducers provides the data and the inversion returns the acoustic properties on the surface. Acoustic holography can determine the sources of noise from a radiating structure, which can help guide a re-design. For example, starting with the discrete equivalent of Equation (71), the field data ϕ E (effectively the sound pressure, see Section 2.1) is related to the layer potential σ S . If we could find a discrete approximation to σ S , then the approximation to the surface potential (sound pressure) and velocity could be found by the matrix-vector multiplication of the discrete equivalent of Equations (53) and (54) and the surface intensity could then be found by the Equation (5). However, as discussed, the solution of the linear system will not yield acceptable results. Even if the number of data points massively exceeds the number of elements, the system is stll effectively underdetermined.
Perhaps the most popular method of resolving the ill-posedness is to use Tikonov regularisation. This involves minimising the residual in system, along with a penalty function that constrains the variability of the solution in some sense. For example, Equation (93) is replaced by where η is a parameter that can be 'tuned' to achieve and acceptable solution and the norm is the 2-norm. A related regularisation method is termed truncated singular value decomposition (TSVD). For example, the singular value decomposition (SVD) of the L ES,k matrix factorises it as follows: where U is a n E × n S matrix, V is n S × n S and Σ is a diagonal matrix containing n S singular values s i , non-negative values in non-decreasing order. In Equation (95) where n * < min (n S , n E ).
There have been several reported implementation of acoustic holography using the methods discussed [237,[329][330][331][332][333][334]. Similar methods have also been developed for finding the impedance of the surfaces in rooms from measured sound pressure data [335,336].

Meshless Methods
One of the main advantages of the boundary element method over domain methods, such as the finite element method is the reduced burden of meshing. With meshless methods, a mesh is not required at all. The simplest concept of a meshless method is the equivalent sources method (ESM) in which it is presumed that the acoustic field is equivalent to a field generated by a finite set of point sources: where the q j are the positions of the sources that are outside of the acoustic domain and the γ j are the unknown source strengths and G k is the appropriate Green's function. For the Dirichlet boundary condition, by matching the Dirichlet data on the boundary ϕ p j for p j ∈ S, gives a linear system of equations, provided there are at least as many source points as there are items of boundary data. Similarly, by differentiating the above equation with respect to the normal to the boundary then approximate solution can be found by matching the values with Neumann data, and using a combination. Usually, more internal points than items of boundary data are chosen and a least-squares solution is sought. The meshless methods avoid the problem of singular integration. However, the ESM is not based on a boundary integral equation formulation and hence it is an alternative to the boundary element method and is therefore beyond the scope of this paper. Lee [337] provides a recent review of the ESM in acoustics. An alternative method for developing meshless methods has more in common with the standard boundary element method. The method can be applied to the standard interior or exterior problem. For the exterior problem, the methods relates to the Schenck or CHIEF method, of Section 4.3.3, but only the integral equations for the internal points are applied, and the integrals are approximated by using only the midpoint value M DS,kφ S = L DS,kv S .
With the number of interior points exceeding the number of boundary points the solution can be found in the least-squares sense. For the exterior problem, the Burton and Miller form is expected to have superior numerical properties: These methods still avoid integration, particularly singular integration and can also be used on the modal problem. There are several reports of implementations of these methods [144,145,[338][339][340][341].
As methods for solving acoustic problems, the meshless methods, are well behind the standard boundary element method in the sense of developing robust software. As with the Schenck or CHIEF method of Section 4.3.3, there is the added issue of determining the placement of the equivalent sources. The methods outlined in Section 3.4.2 could be used in determining whether source points are interior or exterior.

Areas of Discussion
Much of the development work on the boundary element method in acoustic has been on relatively simple shapes and relatively low frequencies. For practical problems, the existing methods must be applicable to more complicated domains and to high frequencies. In this Section, two significantly challenging areas in the acoustic BEM are surveyed. The first is that of efficiency, the relationship between the accuracy achieved and the computational effort. Error analysis of the acoustic/Helmholtz BEM has received significant attention [230,342,343]. Three categories of error arise in the BEM; the discretisation error due to the approximation of the boundary and boundary functions, the quadrature induced error resulting from the numerical integration method and the error in the solution of the linear system of equations. In general, efficiency is maximised by balancing the errors. Much of the focus in the development of the BEM in acoustics is on improving its efficiency so that a full frequency sweep, particularly including the resolution required to model high frequencies, can be achieved with reasonable computer time and memory requirements. In Section 6.1 the efficiency of the acoustic BEM is discussed and methods for improving efficiency are reviewed.
The most valuable acoustic boundary element method-the standard exterior problem, outlined in Section 5.3-has also been found the most difficult to maintain. The combined integral equations of Section 4.3.4 have held the most promise. The earlier work on the acoustic BEM suggested that the coupling parameter was somewhat arbitrary. In terms of scaling up the method, however, the coupling seems to have become a significant issue and, in Section 6.2, the choice of coupling parameter is revisited.

Efficiency
It is often stated that the boundary element method had a significant efficiency advantage over domain methods, such as the finite element method. However, efficiency concerns remain critical, not least for the acoustic boundary element method, in which it has been a focus of research for many decades. Efficiency, relates the accuracy of the output to the computational resources required to achieve it. The computational resources include the computer processing time and/or the memory requirements. One approach for improving the execution time of numerical methods is to use parallel processing, so that instructions are executed simultaneously, rather than sequentially, and such techniques have been applied in the boundary element method [304,344].
In this Subsection, the efficiency of the acoustic BEM is analysed and techniques for improvement are reviewed. Although a range of classes of boundary element methods have been outlined in this paper, the analysis is fairly generic, and can be applied to the boundary-value problems that arise in acoustics. The modal analysis or eigenvalue problem is not considered, but much of the analysis carries over. For the extended problems considered in Section 5, the analysis in this Subsection is only relevant insofar as the BEM is implemented with the wider context.

The Computational Cost of the Acoustic Boundary Element Method at One Frequency
In this Sub-subsection, the boundary element method in acoustics is developed in its most typical way for one frequency. The efficiency of the method is analysed and discussed. More disruptive methods for improving the efficiency are considered subsequently.
As discussed, the boundary element method is composed of two stages, the first stage involving determining the boundary functions and the domain solution in the second stage. In the first stage one or more n·n matrices are formed, where n is the number of boundary solution points (e.g., collocation points, or elements for simple constant elements). In the second stage, usually one n P ·n matrix is formed, where n P is the number of domain points. Hence the storage requirement is O n 2 + n P n complex numbers.
The matrices are normally computed by numerical integration. In general, the number of quadrature points required on each panel would be varied with the size of the panel and the distance between the panel and the point [345]. For each quadrature point, the Green's function would be computed, together with the other geometric information, although the former is likely to incur the greater part of the computational cost. The discretisation could be carried out separately for each operator; however, for efficiency reasons, the computed values should be shared between the operators, as, for example, followed in the author's previous work [61,91]. The diagonal components of the square matrix (or matrices) in the initial stage may be the result of evaluating singular and hypersingular integrals, as discussed, and, although they require special treatment, it is usually more efficient if the quadrature points are similarly the same for all required discrete operators. Lumping together the matrices, for the n × n and n P × n components, let the average number of quadrature points be N Q and let the average cost of the evaluation of the Green's function be C G and the average cost of the geometrical properties be C g , then the total cost (or computer time) of computing the matrices is N Q n 2 + n P n C G + C g .
Clearly C g will be larger for the combined operators, used in the primary stage for exterior problems. Alternatively, for the Schenck method, the matrices are augmented by the discrete operators for the interior points.
Once, the matrix vector system is formed, the next step is to solve it. The most straightforward method of solution of a square system is the Gaussian elimination or LU factorisation method (re-writing a matrix as the product of a lower and upper triangular marices). The overall computer time for the acoustic boundary element method at one frequency can be summarised with the equation where C f represents the cost of a floating-point operation, such as multiplication or division. The non-square system that arises in the Schenck method can also easily be resolved as a square system through pre-multiplying both sides by the transpose of the matrix over which the solution is sought and this is equivalent to the least-squares solution.
In its earlier development, the cost of computing the matrices generally far outweighed the cost of solving the system. However, in the first stage of the boundary element method, the computational cost of setting up the linear system is O n 2 , whereas the cost of solving it (using the methods stated) is O n 3 . For this reason, it has also been understood for a long time that LU factorisation or related methods for solving the linear system was unsustainable, as progress is made towards the finer meshes that are particularly required for high frequency problems. However, it is important also to state the particular advantages that LU factorisation has. LU factorisation is a robust method and, in the case when there are a range of inputs to a problem with a fixed boundary and boundary condition, once the O n 3 factorisation as been stored, it can be used repeatedly again with O n 2 cost.
The fast multiple method (FMM) is a popular method of speeding up the computation of the matrices in the boundary element method and it has been applied to acoustic/Helmholtz problems [220,227,254,284,287,290,320,346]. In this method the Green's function is approximated by local polar expansions that can be translated through the domain, rather than re-computed. The FMM is often used with iterative methods.
In this Subsection alternatives to standard LU factorisation are reviewed. Faster solution methods include the use of hierarchical matrices [190] or panel clustering [347] and iterative methods are also considered. If the solution of the linear system is potentially dominant, in terms of computational cost, then the attention also re-focuses on the integral equation method. For it may be advantageous to choose an integral equation method that results in a linear system that has faster convergence properties, rather than one that is the easiest to apply or the most efficient (such as the collocation method that is highlighted in this work).

The Frequency Factor, Multi-Frequency Methods and Wave Boundary Elements
Typically, in vibration and acoustics, the time-dependent signals are resolved into frequency components. For example, for air-acoustics, the audible range for human being is up to 20 kHz and typically his could be resolved into components with a 10Hz resolution; multiplying the computational cost, considered in the previous Sub-subsection, by around 2000. The prospect of solving acoustic problems over the full frequency range has led to another set of techniques focus on reducing the computational burden by solving a range of frequencies together and these are usually termed multi-frequency methods [297,[348][349][350][351][352][353][354][355]. However, with the nature of an acoustic field, originating typically from structural vibration, is such that structural and acoustic resonances, and the driving profile, dominate the total acoustic response. From a computational point of view, effort should therefore be concentrated on these 'peak' frequencies. Hence, this adaption returns us to the standard mono-frequency problem, or applying the multifrequency method, centred on each peak.
Moreover, in general the frequency of observation is matched in the acoustic and vibratory properties. In practice, as the frequency increases, a higher elemental resolution is required and hence n = O(k) for two-dimensional and axisymmetric problems and n = O k 2 in three dimensions. Although this seems to imply a revision of the mesh at every frequency, it is more likely that one mesh that is suitable for the highest frequency is used throughout the range, or separate meshes are applied in significant sub-ranges of the frequency range, in order that the burden of meshing is proportionate within the overall project.
A review of the actual number of elements required to capture the solution is provided by Marburg [356][357][358]. In general, for the simple elements that are often used, it is considered that 6-10 elements per wavelength is a reasonable guideline. Obviously, higher-order boundary element methods [100,[359][360][361][362][363] would often return the same accuracy with fewer elements. There is significant interest in isogeometric elements, in which the boundary and the boundary functions are modelled with the same basis functions (typically splines), so that the BEM can be more easily integrated with computer-aided design [364][365][366][367][368][369][370][371]. For air acoustics, at high frequencies reaching 20 kHz, the wavelength is less than 2 cm. For example, even for a 10cm cube, the number of elements required to model the high frequencies is in the tens of thousands and for a 1m cube, over a million elements would be required. A variation on the boundary element method has arisen for developing sinusoidal basis functions or wave boundary elements, in order to more accurately model the acoustic functions and the term PUBEM or partition of unity BEM is often used to identify such methods [267,[371][372][373][374], which has similarity with the application of the Treffetz method [292][293][294].

Iterative Methods and Preconditioning
The scalability of the standard BEM in acoustics has been in question for several decades. The computational estimate relates two areas of particular concern, the O(n 3 ) nature of LU factorisation and the O(n 2 ) cost of computing the matrices. Hence, much of the current research on the BEM in acoustics is focused on reducing its computational burden, so that the methods can be casually applied to more significant problems and high frequencies.
The conjugate gradient method was identified as a useful iterative solution method for the acoustic BEM [155,375,376]. The conjugate gradient method and related methods are generally termed Krylov subspace methods and these variants have been significantly tested on the linear systems arising from the boundary element method in acoustics [106,[377][378][379]. In general, if the underlying operator is compact, corresponding to a clustering of the corresponding matrix eigenvalues, then iterative methods, such as the Krylov subspace methods, work well [380]. However, one of the operators, N k , is not compact, the eigenvalues of its matrix equivalent are spread out, and hence the raw iterative methods are only applicable in particular cases, in which the hypersingular operator is not in play. Hence, in the spirit of generality that is sought in this work, the iterative methods are of limited value in solving the untreated equations. However, in the acoustic boundary element method, these iterative methods are usually applied following an intervention with the original linear system, termed preconditioning [381].
In the acoustic BEM, preconditioning has been advocated since its early days. The original concept can be derived from the integral equation formulations. By various substitutions with the boundary integral equations, operator identities can be derived. One of the most useful is the substitution of (74) and (75) into (70) (p ∈ S) (or substituting (63) or (64) into (48) (for p S)) giving The equation is preconditioned, as the N k operator has been eliminated and all the operators are compact. However, the motivation for this technique can be with the elimination of the hypersingular integration, with incidental preconditioning. Methods based on (96) along with L k N k = M t 2 k − 1 4 , that is obtained through the substitution of (53) and (54) into (61) or the substitution of (71) and (72) into (73) (p ∈ S), are often termed the Calderon equations [382,383]. Methods that make use of this substitution have been developed [155,157,241]. Alternative approaches in developing well-conditioned integral formulations for iterative solution have also been the subject of research [384][385][386][387][388]. Langou [389] develops methods for the iterative solution of linear systems with multiple right-hand sides, echoing the stated advantage of the LU factorisation method, discussed in Section 6.1.1.
Matrix preconditioners are often based on constructing an approximate inverse and following a fixed point/ defect correction/ contraction method. Introducing a preconditioning matrix also introduces an O(n 2 ) matrix multiplication at each step, but this can be reduced to O(n) if the approximate inverse matrix is sparse or banded [378,380]. Methods based on an incomplete LU factorisation [390] have been tested in the acoustic BEM [391]. A similar approach is the construction of a DtN (or DN) map (Dirichlet to Neumann) or NtD (or ND) map, also termed an on surface radiation condition (OSRC), as the approximate inverse [190,[392][393][394].

The Coupling Parameter
The direct and indirect integral formulations of the exterior problem that that were free of the characteristic frequencies or non-uniqueness, were introduced in Section 4.3.4. The equations were defined with a coupling parameter µ, which weighs the contribution of derivative equation with the original equation in the direct method, or between the double layer potential and the single layer potential equations in the indirect method. Mathematically, the coupling parameter has a non-zero imaginary part that ensures that the equations are free of characteristic wavenumbers, and, in general, µ is presumed to be an imaginary number. In accepting that µ is imaginary, the starting point is −∞ < Im(µ) < ∞. If |µ| is very small or very large then one or the other of the original equations would dominate, and the issues with the dominant equation would be evident with the hybrid equation, and this narrows its range; 0 |µ| ∞. The argument quickly shifted from introducing µ in order to avoid the potential catastrophic errors in the original equations of Section 4.3.1 to considering an optimal value. For several decades, based originally on the works of Kress [395][396][397][398] and the further works of Amini [399], there has been a strong recommendation, with supporting research, that µ = i k is reasonably close to the 'optimal' choice, as least for the simple boundaries, like spheres. Obviously i k is unsuitable for low wavenumbers, as the parameter would be large, violating the condition |µ| ∞, and the second equation would dominate and a cap on its value is recommended, for example, The term 'optimal' in these paragraphs, refers to the condition of the combined operators over which the equation is being solved. The stronger the condition of the operator, the less able it is to magnify the solution, or magnify its error. The original equations, on their own, have 'spikes' of ill-conditioning around their respective eigenvalues. The rationale is that in combining the equations, these spikes are significantly reduced, and the condition of the combined operator steers an even keel through the frequencies. It has been shown, for example for a sphere that i k provides a generally well-conditioned combined operator. A recent paper by Zheng [400] showed that each eigenvalue of the combined operator follows a loop-shaped locus as µ varies, joining the real axis when |µ| = 0, ∞, and supporting µ = i k , as this reasonably approximated the point on the locus that was furthest from the real axis. Zheng terms the inclusion of the derivative equation as 'adding damping', relating the analogy with, for example, Equation (91), wherein the damping term similarly shifts the eigenvalue off the real axis; the combined integral equations are not free from irregular frequencies, rather they are simply moved from the real axis into the complex plane.
There is a simple case for the parameter with an inverse proportionality with the wavenumber. Considering the three-dimensional case, the Green's function (52) is the kernel of the L k operator, and its magnitude does not change with k e ikr = 1 . Hovever its derivative, a factor in the kernel of the M k and M t k operators, is O(k). Hence, combining the operators of L k with M k or M t k , as in Section 4.3.4, a parameter that is inversely proportional to k equalises the contribution from the two operators. Similarly, differentiating again, the kernel of the N k is O(k 2 ), and the same parameter has a similar function, when M k or M t k are combined with it. A recent paper by Marburg [401] considered principally the sign of µ. Marburg noted that many papers had inadvertently used µ = − i k as the coupling parameter. Given that its value was not thought to be critical, it might be thought that this would work as well as µ = i k . However, Marburg demonstrated that i k usually returns significantly more accurate results, and much faster convergence with the iterative methods for solving the systems of equations. The parameter µ = − i k was not used out of choice in the papers reviewed therein, but rather through the confusion of the signs that underlie the basic definitions in the mathematical model. These results have been confirmed by Galkowski [116].
If L k is coupled with its derivative, as it is in the Burton and Miller equation, then the following equation is obtained, for one side of the equation (L k + µ(M t k + With µ = i k there is significant cancellation in r + µikr ∂r ∂n p when ∂r ∂n p ≈ 1, possibly with a diagonalising effect on the operator, and this is also presumed for the other combinations of operators. Another approach is to consider the condition of the matrices that arise in the BEM [402][403][404]. Earlier, it was discussed that the N k operator was a pseudo-differential operator, rather than an integral operator. As a result of this, the N SS,k increases with the number of elements n s , whereas for the other integral operators, the norm of the matrix is independent of the number of elements. This is also supported by Equations (45) and (46); the diagonal components of the N SS,k matrix are O 1 h = O(n S ) for two-dimensional and axisymmetric problems and O √ n S for three-dimensional problems. It follows that N SS,k = O(n S ) or O √ n S . The scene is therefore set for two conflicting interests in choosing the coupling parameter µ = i k is near-optimal in conditioning the combined operator, as a result of mathematical analysis whereas µ = i n S or µ = i √ n S (put simply), balancing the matrix norms, as a result of numerical analysis. Without the latter correction, the N SS,k matrix will apparently provide an increasing dominant potential for (quadrature-induced) error as the number of elements increases.
However, in Section 6.1, the modern form of the combined method was outlined, involving pre-conditioning and iterative methods for solving the linear systems. For example, with the preconditioner in (96), the N k operator is avoided and, in such cases, the analysis of the previous two paragraphs is not applicable.

Concluding Discussion
The main purpose of this paper is to encapsulate the modern scope of the boundary element method in acoustics; how it can be adopted for general 2D, 3D and axisymmetric dimensions, interior, exterior and half-space problems, thin shells and modal analysis. The standard BEM can be directly applied to an enclosed domain with cavities or to multi-boundary problems in an exterior domain. Boundary element methods are founded on a boundary integral equations formulations and these can be fused to form additional methods: thin shells can be used to model discontinuities in an interior or exterior domain defined by closed boundaries, half-space problems can be modelled by the Rayleigh integral or modifying the Green's function, domain methods can be linked to the BEM so that vibro-acoustic and aero-acoustic problems can be tackled and-through regularization of the BEM equations-inverse problems can be addressed.
The boundary element method is a numerical method that only becomes a potentially useful tool for solving real-world problems in acoustic engineering when it has been implemented in software. Clearly, it is important that the executing code fits within the memory of the computer and completes the computation in reasonable time, and efficiency issues have been considered in Section 6.1. There are issues of generality and, in that regard, this work has focused on the generalised Robin boundary condition and general topologies. Reliability and maintenance are also very important issues in software development, and hence much of the focus has been on the standard exterior acoustic problem, which has had issues and solution approaches documented for more than half a century. The robustness of the BEM, particularly considering the validity of a defined elemental boundary, has had little attention in the field, but methods for assisting with this are summarised in Section 3.4.2.
There is much commonality across the various topologies to which the boundary element method can be applied in acoustic/Helmholtz problems; the equations for different dimensional problems are literally the same, with a change in the definition of the Green's function and line integrals become surface integrals when we move from 2D to 3D. All the core equations, whether they are interior, exterior, shell problems or the straightforward half-space problems, have the same essential operators within a chosen dimensional space. Hence there is significant scope for component-wise development of software, adopting a 'library' approach and unifying the method.
Significant issues have always surrounded the N k operator. This operator was initially included in the combined integral equation formulations for the exterior problem, outlined in Section 4.3.4. Firstly, for surface (collocation) points its integral definition is hypersingular, which is difficult to interpret and evaluate, and perhaps for that reason alone, the alternative Schenck/CHIEF method is the preference of many. Once this significant issue is surpassed, the N k is not compact, it is a pseudo-differential operator and this causes further issues for solving over by iterative methods and hence preconditioning has been introduced to circumvent this. Preconditioning can eliminate the hypersingular operator N k for the standard exterior problem (96). However, N k is still required for shell elements and hence it has to be included in a general library.
Much of the current research theme in the area, put simply, is to shunt the BEM in acoustics from its traditional comfort-zone of problems with ∼10 3 elements to modern problems, to include high frequencies, with ∼10 6 elements and from straightforward problems to complicated applications [190,238,405]. The computational bottle-neck in the traditional BEM is in the use of LU factorisation or similar methods for solving the linear system. Hence the LU factorisation is the first necessary casualty of this move and iterative methods are favoured, although this has been presaged for several decades. The shift to ∼ 10 6 elements could also have a significant demand on computer memory and hence the expectation of storing matrices must also be relaxed. With iterative methods, the effect of the approximation of the matrix components is more controllable and methods such as the fast multipole method and panel clustering are more able to broker the issues of storage and computation time to this end.
Probably the majority of research on the BEM in acoustics is on the seemingly intractable exterior problem. Our expectation that the solution in its infinite domain can be thread through the boundary is ultimately a questionable one, as the method is scaled up. The combined methods, and most prominently, the Burton and Miller method, have held centre-stage in this endeavour, seen by most as numerically superior to the Schenck/CHIEF method [238]. However, the two erstwhile competing methods may have to be married in our best efforts to achieve ∼10 6 elements and high frequencies; using the equations from interior points to provide more stability in the Burton and Miller equation, as they did originally with the elementary equations to form the CHIEF method. The system could be augmented further with directional derivates of Equation (70) for points in the interior D.
The combined methods throw up two issues, one is the inclusion of the N k operator, as discussed, and the other is the coupling parameter. The systems of equations arising in the BEM using combined methods require that the N k operator is preconditioned in order to achieve convergence with iterative methods, as discussed in Section 6.1. The choice of coupling parameter is discussed in Section 6.2. However, the approach to choosing the coupling parameter is already to potentially optimise the system and hence it is itself a pre-conditioner, as pointed out in Betcke et al. [190]. In Harris and Amini [406], the coupling parameter is generalised from a singular value to a function over the surface. Hence, another marriage is proposed. As the coupling parameter and the preconditioner have a similar purpose, then the Burton and Miller Equation (76), for example, could take the following form where U k is the preconditioning operator, fusing the coupling parameter and the preconditioner into one entity. The objective then is to determine U k , or an analogous preconditioning matrix, U SS,k . For example, following the method in Equation (96), U k = i k L k ; however, in general, the goal is to set U k or U SS,k to best-prepare the system for iterative solution.