## Abstract

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.

## 2. Introduction

Central to acoustic wave theory is solving the Helmholtz equation for the pressure field in the frequency domain. The boundary integral method is commonly used because of the reduction in spatial dimension. However, it is well known that the numerical solution of the boundary integral equation for external problems can become inaccurate when the wavenumber is close to one of the eigenvalues of the internal problem [1,2]. The issue has been addressed by two common methods. The CHIEF method owing to Schenck [3] imposes an additional constraint on the solution of the boundary integral equation by requiring it to vanish at selected positions inside the boundary to suppress the resonant solution that has no physical significance in the external problem. An alternative approach proposed by Burton & Miller [4] involved solving a hypersingular integral equation in which the original boundary integral equation is the real part and the boundary integral equation for the normal derivative is the imaginary part. In both cases, the equations to be solved contain mathematical singularities associated with the conventional formulation of the boundary integral equation and much effort since then has been concerned with the expeditious and efficient treatment of these singularities [5].

Recently, we reformulated the boundary integral solution of the Laplace equation: ∇^{2}*ϕ*=0, and the Stokes equation of fluid mechanics whereby all singular terms in the integrals are removed analytically [6]. This regularization of all the singular behaviour on the boundary means that the surface integrals can be evaluated using any convenient quadrature method. A significant practical consequence of this is that high numerical precision can be obtained with mixed boundary conditions [7] and for problems with multiscale characteristics such as those with boundaries that are very close together compared with their characteristic dimensions or where the boundaries possess extreme geometric aspect ratios [8].

In this work, we apply this boundary regularization to the Helmholtz equation that removes much of the technical effort needed to use linear or quadratic surface elements to represent boundaries and results in significant improvement in the accuracy in the evaluation of surface integrals. Also, it is no longer necessary to calculate the solid angle at each node—a complexity that can discourage the use of higher order surface elements. As a result, far fewer degrees of freedom are needed to achieve the same precision that in turn translates to a significant decrease in computational time—demonstrated here by numerical examples. When the high precision of this boundary regularized formulation is combined with the CHIEF method [3], the possibility of numerical error arising from the resonance solution becomes negligible in practice. The framework given here is applicable for both external and internal problems, for Dirichlet, Neumann or mixed boundary conditions and for radiation as well as scattering problems.

Although examples motivated by problems in acoustics have been used to demonstrate the theoretical formulation and practical numerical advantages, this method of de-singularizing boundary integral equations can be applied in other contexts such as the Laplace equation [8], the equations of hydrodynamics [6,7] and linear elasticity or any equation that belongs to the Moisil–Theodorescu system [9].

## 3. 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]
*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’, **x** and **x**_{0} are on the boundaries. For simplicity, we assume that the surfaces *S* and *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 *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
*G*_{0}(**x**_{0},**x**)=1/|**x**−**x**_{0}|. The standard way to remove the singularity arising from ∂*G*_{0}/∂*n* on the left-hand 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 (3.2) for [*q*(**x**)−*q*(**x**_{0})] becomes
*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 *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 (∂*p*/∂*n*)_{0} at **x**_{0}
*g*(**x**) and *f*(**x**) satisfy the Helmholtz equation with the following boundary conditions:
**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
*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 *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**):
**x**_{d} is any convenient point *outside* the solution domain in order to ensure *ψ*(**x**) is non-singular within the solution domain, with *a*≡|**x**_{0}−**x**_{d}|, *r*_{d}≡|**x**−**x**_{d}| and *b*≡(**x**_{0}−**x**_{d})⋅**n**(**x**_{0})/|**x**_{0}−**x**_{d}|≠0 (figure 1*a*). For this particular choice of *g*(**x**) and *f*(**x**), the integral over the surface at infinity,

In the limit *k*=0, the Helmholtz equation reduces to the Laplace equation with

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.

## 4. 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}
*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:

## 5. 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 *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]: *p*(*r*)=−[(i*ωρV* *R*^{2})/(1−*ikR*)] *e*^{ik(r−R)}/*r*.

In figure 1*b*, 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:

CBIM Constant: The value of the unknown pressure associated with every

*flat*triangular surface element is assumed to be*constant*within the element.CBIM Linear

*c*_{0}: The points on the surface of the sphere on which the pressure is to be found are at the vertices of*planar*triangular area elements with linear shape functions used to represent the boundary. The solid angle*c*_{0}in (3.1) is calculated according to the local geometry around each point**x**_{0}[11].CBIM Linear

*c*_{0}=2*π*: This an often used approach that assumes the*incorrect*value of*c*_{0}=2*π*for the solid angle at every node even though the surface tangent plane is clearly not defined at the vertices of the surface elements.BRIEF Constant: Using BRIEF with the assumption that the unknown pressure associated with every

*flat*triangular surface element is*constant*within the element.BRIEF Linear: Using BRIEF where the points on the surface of the sphere on which the pressure is to be found are at the vertices of

*planar*triangular area elements with linear shape functions used to represent the boundary.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 1*b* 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 2*a*, 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 2*b*. Note that to incur a maximum error exceeding 1%, *r*_{cf}/*R* has to be within 0.004 of the node at

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 3*R* 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 3*a*, 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 3*b*.

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 3*R*, at separation *h*/*R*=0.001 as illustrated in figure 3*b*. In figure 4*a*, 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 4*b*, 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 3*a*. By contrast, the numerical errors of the CBIM are already quite evident from the roughness of the pressure contours at *ρ*∼0.7*R*, 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 4*b* demonstrate the utility of the BRIEF in furnishing accurate results in both the near and far field.

## 6. 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].

## Author contributions

S.Q. carried out the numerical implementation of the theory, contributed to the theoretical development and the draft of the manuscript; E.K. contributed to the theoretical and algorithm development and the draft of the manuscript; B.C.K. commented on the manuscript; D.Y.C.C. contributed to the theoretical development, the design of the project and the draft of the manuscript. All authors gave final approval for publication.

## Funding statement

An Australian Research Council Discovery Project grant to D.Y.C.C.

## Competing interests

We declare we have no competing interests.

## Acknowledgements

D.Y.C.C. is a visiting scientist at the Institute of High Performance Computing and an adjunct professor at the National University of Singapore. He also receives in-kind support from the Particulate Fluids Processing Centre, University of Melbourne.

## Appendix A. Derivation of the integrals over the surface at infinity

In this appendix, we derive analytic expressions for the integrals over the surface at infinity,

In (3.8) and (3.9), the integrals at infinity vanish owing to the Sommerfeld radiation condition [14,15]
*g*≡*g*(**x**), *f*≡*f*(**x**), ∂*g*/∂*n*≡∇*g*(**x**)⋅**n**(**x**), ∂*f*/∂*n*≡∇*f*(**x**)⋅**n**(**x**), *G*≡*G*(**x**_{0},**x**), ∂*G*/∂*n*≡∂*G*(**x**_{0},**x**)/∂*n* and d*S*≡d*S*(**x**).

Both *g* and *f* vanish as 1/*r*_{d}∼1/*r* at the surface at infinity, *G* and ∂*G*/∂*n* vanish as 1/*r*, when *r* can be omitted because the surface *r*^{2} at infinity. Consequently, only terms multiplying ∂*G*/∂*n* that vanish as 1/*r* can give contributions to the integrals.

The sine and cosine terms in *g* and *f*, see (3.7), can be rewritten in complex notation as
*g* and *f* as

It can be easily seen that the terms proportional to *e*^{ik(rd−a)} in the integrals of (A 2) and (A 3) remain finite and the individual terms cancel out exactly. However, for the terms proportional to *e*^{−ik(rd−a)}, this is not the case and they contribute to the integrals in (A 2) and (A 3) to give
*rr*_{d} of the integrand can be approximated by *r*^{2} as *r*≈*r*_{d}. For the numerator, we need to invoke the cosine rule and a Taylor expansion to *r*−*r*_{d}+*a*, which gives *α* being the angle between the line segments *r*_{d} and *a*. We then get

In (4.2), we need to evaluate the following two integrals over the surface at infinity, *G*^{p}≡*G*(**x**_{p},**x**) and ∂*G*^{p}/∂*n*≡∂*G*(**x**_{p},**x**)/∂*n*. Here we can rewrite *g* and *f* in (A 6) and (A 7) as
*a*_{p}=|**x**_{p}−**x**_{d}|. Following the same procedure we used to evaluate (A 2) and (A 3), the key integral that needs to be determined for calculating of (A 13) and (A 14) is
*r*_{p}=|**x**−**x**_{p}|. By using (A 12), we then find

- Received December 14, 2014.
- Accepted December 18, 2014.

© 2015 The Authors. Published by the Royal Society under the terms of the Creative Commons Attribution License http://creativecommons.org/licenses/by/4.0/, which permits unrestricted use, provided the original author and source are credited.