Boundary regularized integral equation formulation of the Helmholtz equation in acoustics

A boundary integral formulation for the solution of the Helmholtz equation is developed in which all traditional singular behaviour in the boundary integrals is removed analytically. The numerical precision of this approach is illustrated with calculation of the pressure field owing to radiating bodies in acoustic wave problems. This method facilitates the use of higher order surface elements to represent boundaries, resulting in a significant reduction in the problem size with improved precision. Problems with extreme geometric aspect ratios can also be handled without diminished precision. When combined with the CHIEF method, uniqueness of the solution of the exterior acoustic problem is assured without the need to solve hypersingular integrals.


Boundary regularized integral equation formulation
To illustrate the boundary regularized integral equation formulation (BRIEF) of the Helmholtz equation, we draw on the problem of calculating the acoustic pressure wave generated by a vibrating boundary specified by a closed surface, or more generally a set of closed surfaces denoted by S. The acoustic pressure, p(x) outside S, obeys the Helmholtz equation: ∇ 2 p + k 2 p = 0, where k is the wavenumber. The solution can be found by solving the conventional boundary integral equation that follows from using Green's second identity [10] where G(x 0 , x) = e ikr /r, with r ≡ |x − x 0 |, is the fundamental solution of the Helmholtz equation. These integrals are taken over the set of closed surfaces, S, specified by the problem and also over the 'surface at infinity', S ∞ , at which the Sommerfeld radiation condition is applied (see appendix A). The points x and x 0 are on the boundaries. For simplicity, we assume that the surfaces S and S ∞ have well-defined tangent planes at all points on the surfaces. We refer the reader to our other work [8] for the treatment of surfaces with sharp edges and vertices. The normal derivatives are defined by ∂p/∂n ≡ ∇p · n(x) and ∂G/∂n ≡ ∇G · n(x), where n(x) is the unit normal vector at x pointing out of the solution domain. The solid angle c 0 at x 0 is equal to 2π if the tangent of the boundary at x 0 is defined, otherwise it has to be calculated from the local geometry [11]. Equation (3.1) provides a relation between the function p and its normal derivative ∂p/∂n on the surfaces S and S ∞ . For Dirichlet (Neumann) problems, p (∂p/∂n) is specified on the surfaces, and (3.1) can be solved for ∂p/∂n (p). Despite the fact that the physical problem may be well behaved on the boundaries, the mathematical singularities in (3.1) at x = x 0 owing to G and ∂G/∂n require careful treatment. To quote Jaswon [12]  and Symm [13]: 'Most integral equations of physical significance involve singular, or weakly singular, kernels, thereby hampering the procedures of both theoretical and numerical analysis'-a statement that is still valid up to this day.
Here we seek to eliminate these singularities analytically by exploiting the linear nature of (3.1) as follows.
But before we tackle (3.1) that corresponds to the Helmholtz equation for p(x), consider first the simpler case of the Laplace equation: ∇ 2 q(x) = 0 for which the corresponding boundary integral equation is where G 0 (x 0 , x) = 1/|x − x 0 |. The standard way to remove the singularity arising from ∂G 0 /∂n on the lefthand side is to note that the constant q(x 0 ) also satisfies the Laplace equation with a vanishing normal derivative on the boundaries, so that equation This is known as the 'constant value subtraction'. However, the integral on the right-hand side of (3.3) still involves an integration over the singularity from G 0 and is normally handled by a change to local polar coordinates [5]. However, this singularity can also be removed analytically [8]. Implicit in this approach is the generally valid assumption that as x → x 0 , [q(x) − q(x 0 )] vanishes as |x − x 0 | or faster.
In this paper, we extend the approach [8] developed for the Laplace equation to remove all singularities in the boundary integral equation (3.1) for the Helmholtz problem. However, as a constant is not a solution of the Helmholtz equation, we need to construct an auxiliary function ψ(x) for a given value of x 0 in (3.1) that satisfies the Helmholtz equation that is analogous to the constant, q(x 0 ), in the case of the Laplace equation. We choose ψ(x) to be linear in the pressure p(x 0 ) and its normal derivative Now provided the general functions g(x) and f (x) satisfy the Helmholtz equation with the following boundary conditions: and with n 0 ≡ n(x 0 ) being the outward unit normal at x 0 , ψ(x) will also satisfy the same boundary integral equation as (3.1). By taking the difference between the boundary integral equations for p(x) and ψ(x), we obtain The conditions imposed on g(x) and f (x) at x 0 in (3.5) will remove both the singularities owing to G and ∂G/∂n in the new boundary integral equation in (3.6) under the generally valid assumption that as x → x 0 , [p(x) − ψ(x)] and [∂p/∂n − ∂ψ/∂n] vanish as |x − x 0 | or faster. Therefore, if g(x) and f (x) satisfy (3.5), all singular behaviour in (3.6) owing to G and ∂G/∂n will be removed. The formal proof is a straightforward generalization of that given in [6] for the Laplace equation. Equation (3.6) is the final key result for the singularity-free [6] BRIEF of the solution of the Helmholtz equation that replaces the conventional boundary integral equation in (3.1). The integrals are now completely free of singularities [6] and can be evaluated by any convenient integration quadrature. Given values for p (Dirichlet) or ∂p/∂n (Neumann) or a relation between these two quantities on the boundary, S, (3.6) can be readily solved. Also the term involving the solid angle, c 0 , in the conventional boundary integral equation, (3.1), has now been eliminated. Thus, higher order area elements, linear, quadratic or splines, can be used to represent the surface to improve numerical accuracy without the need to calculate the solid angle at each node. For our key result in (3.6), we see that apart from having to satisfy (3.5), there is considerable freedom in selecting the functional forms of g(x) and f (x) and hence ψ(x). As a specific example, we can choose the following for g(x) and f (x): where x d is any convenient point outside the solution domain in order to ensure ψ(x) is nonsingular within the solution domain, with a ≡ |x ). For this particular choice of g(x) and f (x), the integral over the surface at infinity, S ∞ in (3.6) can be found analytically using the Sommerfeld radiation condition [14,15] (see appendix A for details) and (3.9) In the limit k = 0, the Helmholtz equation reduces to the Laplace equation with g(x) , and (3.6) will become a boundary regularized integral equation for the solution of the Laplace equation given earlier [8]. This is an example of how the singularity in the integral on the right-hand side of (3.3) may be removed analytically.
We note that the conventional boundary integral formulation of the solution of the Helmholtz equation, or for that matter solutions of the Laplace equation or the equations of hydrodynamics or elasticity, gives rise to singularities in the integrands of the surface integrals. These mathematical singularities originate from the fundamental solution of the equations used in Green's second identity, whereas the physical problem is perfectly well behaved on the boundaries. Therefore, it is not too surprising that these mathematical singularities can be completely removed by subtracting a related auxiliary solution of the governing equation: ψ(x) in the case of the Helmholtz equation considered here. Apart from having to satisfy the governing differential equation and some mild constraints, see (3.5), there is considerable flexibility in choosing the precise functional form of the auxiliary solution or any free parameters that it may contain-such as the value of x d in (3.7). A broad physical interpretation is that we have constructed a pressure field, ψ(x), that cancels the value of p and ∂p/∂n at x 0 .
One well-studied uniqueness issue is that the solution of the boundary integral method, (3.1), with a closed boundary is identical to that obtained from solving directly the differential form of the Helmholtz equation: ∇ 2 p + k 2 p = 0, except at a discrete set of values of k. These k values correspond to the resonant wavenumbers of the closed boundary [3,4]. The solution of the BRIEF, (3.6), also shares this feature. The resonant or homogeneous solutions corresponding to the resonant values (eigenvalues) of the closed boundary will always emerge from solving the integral equation. But as we shall see in §5, the higher precision that the BRIEF can attain will alleviate much of the practical numerical difficulty.

Accurate evaluation at field points near boundaries
The BRIEF of the solution of the Helmholtz equation also offers an accurate and numerically robust method to calculate the solution at field points close to boundaries. Indeed, the loss of precision owing to the near singular behaviour in the evaluation of the requisite integrals that arise in the conventional boundary integral method (CBIM) is often more difficult to deal with than the singularities on the boundaries.
To evaluate the solution p(x p ) of the Helmholtz equation at a point x p in the solution domain, we first use the conventional boundary integral equation to find [p(x) − ψ(x)] at x p , with x 0 in (3.4) taken to be the node on the boundary closest to x p The near singular behaviour of the surface integrals owing to G and ∂G/∂n as x p approaches the boundary, S, can be now removed by subtracting the boundary regularized boundary integral equation (3.6) from (4.1). Then using (3.4) for ψ(x), we obtain a numerically robust expression for p(x p ) that is free of near singular behaviour irrespective of the distance from x p to any boundary:

Examples from acoustics
We use examples drawn from acoustics to illustrate the implementation and advantages of the BRIEF of the Helmholtz equation relative to the CBIM. The first classic problem is the pressure field outside a radiating sphere of radius, R, with a prescribed time harmonic, radially symmetric surface normal velocity v n = V exp(−iωt) [3], giving the boundary condition: ∂p/∂n = −iωρV, with k = ω/c where ρ is the density and c the speed of sound in the region outside the sphere. This is a Neumann problem because p(x) is to be found for prescribed values of the normal derivative ∂p/∂n. The analytic solution of this problem is [3]: In figure 1b, we compare the maximum relative error of the following six different numerical solutions of the above problem at kR = π/2 for varying numbers of degrees of freedom or unknowns on the surface of the sphere:   .7) is at the origin of the sphere. Table 1. Variation of the condition number of the BRIEF for a pulsating sphere of radius R for two values of kR close to resonance at kR = π , using constant, linear or quadratic surface elements to represent the sphere. The constant element representation predicts an incorrect resonant value very close to kR/π = 1.0038, see also figure 2.

BRIEF Quadratic:
Using BRIEF where the points on the surface of the sphere on which the pressure is to be found are at the nodes of quadratic 6-noded triangular area elements with quadratic shape functions used to represent the boundary.
It is clear from figure 1b that results obtained from the BRIEF are at least as good as, and in the case of BRIEF Quadratic far superior to, any results from the CBIM. Depending on the degrees of freedom, the relative error is smaller by 1-3 orders of magnitude; or BRIEF Quadratic can achieve the same precision as any CBIM with about one-tenth the number of degrees of freedom. In all BRIEF approaches, there is no need to deal with singular integrals or to compute the solid angle c 0 at any node. Moreover, the poor performance of the erroneous CBIM Linear c 0 = 2π approach is demonstrated.
In figure 2a, we compare the BRIEF using constant, linear and quadratic elements for a single radiating sphere considered at wavenumbers very close to the first resonant value, kR = π [3]. We observe that: (i) using the BRIEF with quadratic (BRIEF Quadratic) elements, the maximum relative error is not significant until |kR − π | < 10 −4 and (ii) when the BRIEF is used with constant elements (BRIEF Constant), the resonant wavenumber is located incorrectly. This is because the spherical shape is not well represented by constant elements and resulted in a slightly different resonant wavenumber, even though the polyhedra that represent the sphere have the same volume. Even in this case, the resonant solution is only evident when k is within 1% of the (incorrect) resonant value.
In table 1, we compare the condition number at various values of kR. We see that using BRIEF Linear or BRIEF Quadratic, the condition number at kR very close to the resonant value can be reduced by a factor of 60 when compared with BRIEF Constant. This demonstrates the numerical stability that can be achieved using BRIEF Linear or BRIEF Quadratic.
Next, we study the error when BRIEF Quadratic is combined with one CHIEF point located at r cf inside the sphere. With kR = 2π the lowest resonant mode has a node at r = R/2 [3]. The maximum, mean and median errors when the CHIEF point r cf is positioned close to this node at R/2 are compared in figure 2b. Note that to incur a maximum error exceeding 1%, r cf /R has to be within 0.004 of the node at r/R = 1 2 .   Thus in practice, the BRIEF Quadratic approach plus one CHIEF point is very unlikely to encounter problems with the resonant solution.
In a second example, we illustrate the unique ability of the BRIEF in being able to handle boundaries with extreme geometric aspect ratios by considering two radiating spheres that are nearly touching. We calculate the pressure field owing to these spheres of radii R and 3R at a distance of closest approach, h/R = 0.001, with kR = π/2 (figure 3). The spheres have identical time harmonic radially symmetric normal velocities as in the previous example, but are out of phase by π . In figure 3a, we show the variation of the magnitude of the pressure, |p|, along the longitudes of the two spheres that pass through the point of closest approach at z = 0. Near this point, the result from BRIEF Linear remains continuous as expected, whereas the result from CBIM Linear exhibits large errors that have their origin in the adverse influence of the nearly singular behaviour of the kernel of one surface on an adjacent surface that is very close by. The relative difference in the pressure magnitude |p| obtained using BRIEF Linear and CBIM Linear is shown in figure 3b.
A comparison of the accuracy in the magnitude of the pressure obtained using the BRIEF and the CBIM is given in the contour plots for two spheres of radius R and 3R, at separation h/R = 0.001 as  figure 4a, we compare the variation of the magnitude of the pressure in the median plane, z = 0, in the neighbourhood of the osculating point and in figure 4b, pressure variations in the far field. Around ρ ≡ (x 2 + y 2 ) 1/2 ∼ 0, the BRIEF Linear results remain smooth as expected from the results in figure 3a. By contrast, the numerical errors of the CBIM are already quite evident from the roughness of the pressure contours at ρ ∼ 0.7R, and these errors become unacceptably large for ρ < R/2. This is consistent with the comparison shown in figure 3. The pressure variations in the far field for the two spheres in the plane y = 0 shown in figure 4b demonstrate the utility of the BRIEF in furnishing accurate results in both the near and far field.

Conclusion
We have developed the BRIEF of the solution of the Helmholtz equation given by (3.6) that replaces the CBIM given by (3.1). As all integrands in the BRIEF contain no singular behaviour and the term containing the solid angle has been removed, any convenient quadrature method can be used to evaluate the surface integrals thus facilitating the use of more accurate linear, quadratic or spline elements to represent the boundaries.
The BRIEF also has the unique ability to handle problems with extreme geometric aspect ratios such as where boundaries are very close together or where boundaries possess very different characteristic length scales. For such cases, the absence of singular behaviour in the BRIEF means that high precision is always maintained whereas the inherent singularities in the integral of the CBIM will give rise to unavoidable deteriorations in numerical precision.
The absence of singular behaviour in the surface integrals means that the BRIEF provides a simple and numerically robust way to compute the solution at field points located near boundaries using (4.2). This is an important advance in the accurate evaluation of the solution near surfaces or within small gaps between surfaces where the errors in the CBIM become unacceptably large using the same surface mesh. This is an advantage of the BRIEF that is unmatched by the CBIM without very significant additional numerical and analytical effort.
In the derivation of the key result, (3.6), we have assumed that the boundary S is sufficiently smooth whereby the tangent plane is defined at all points on the surface. For points on boundaries at which the tangent plane is not uniquely defined, for example at sharp corners and edges, one can still implement BRIEF using the double node technique to construct the requisite system of linear equation [8].
From numerical case studies drawn from acoustics considered here, the use of linear elements or quadratic elements with the BRIEF can lower the relative error by a factor of 10 to 1000. Conversely, the same precision as the CBIM can be obtained even when the number of degrees of freedom is reduced by a factor of 10 to 100 or more. This represents a very significant speed up in computational time.
Of particular relevance to acoustic problems, the much higher precision results that can be obtained using the BRIEF also means that the uniqueness issue associated with the resonant or normal mode solutions of the Helmholtz equation around closed surfaces is very unlikely to arise. This is simply because the numerical value of the wavenumber must be within a relative deviation of smaller than 10 −3 from the resonant value before the resonant solution can start to contribute to the answer. Indeed, if the CHIEF method of Schenck [3] is also used, the uniqueness issue will not arise in practice as the CHIEF node must now lie within a relative deviation of smaller than 10 −3 of an internal resonance node before the BRIEF-CHIEF combination might break down. The BRIEF is also much simpler than the Burton-Miller [4] approach because no hypersingular integrals are involved.
Given the above advantages of the BRIEF, the case for its wide adoption is very compelling not only for the Helmholtz equation that has applications ranging from wave phenomena in acoustics to electromagnetics but also for boundary integral solutions of the Laplace equation [8] and hydrodynamic problems [7]. r − r d + a, which gives r − r d + a ≈ a(1 + cos α) with α being the angle between the line segments r d and a. We then get This leads immediately to (3.8) and (3.9). In (4.2), we need to evaluate the following two integrals over the surface at infinity, S ∞ : where G p ≡ G(x p , x) and ∂G p /∂n ≡ ∂G(x p , x)/∂n. Here we can rewrite g and f in (