Boundary regularised integral equation formulation of the Helmholtz equation in acoustics

TL;DR

This paper introduces a boundary integral method for the Helmholtz equation in acoustics that removes singularities analytically, enabling higher accuracy and efficiency in modeling acoustic radiation problems, even with complex geometries.

Contribution

It develops a singularity-free boundary integral formulation for the Helmholtz equation, improving numerical precision and computational efficiency in acoustic simulations.

Findings

Enhanced numerical accuracy with higher order elements

Effective handling of complex geometries and aspect ratios

Guaranteed solution uniqueness with the CHIEF method

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 due 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.

Tables

Table 1. Table 1: Variation of the condition number of the BRIEF for a pulsating sphere of radius R 𝑅 R for two values of k ​ R 𝑘 𝑅 kR close to resonance at k ​ R = π 𝑘 𝑅 𝜋 kR=\pi , using constant, linear or quadratic surface elements to represent the sphere. The constant element representation predicts an incorrect resonant value very close to k ​ R / π = 1.0038 𝑘 𝑅 𝜋 1.0038 kR/\pi=1.0038 - see also Figure 2.

	$kR/\pi$ = 0.95	$kR/\pi$ = 1.0038
Constant element	6.0	6000
Linear element	6.5	81
Quadratic element	6.5	100

Equations

$$c_{0}p({\mathbf{x}}_0)+{\int }_{S+S_{\infty }}p(\mathbf{x})\frac{\partial G({\mathbf{x}}_0,\mathbf{x})}{\partial n}\text{d}S(\mathbf{x})={\int }_{S+S_{\infty }}\frac{\partial p(\mathbf{x})}{\partial n}G({\mathbf{x}}_0,\mathbf{x})\text{d}S(\mathbf{x}),$$

$$c_{0}q({\mathbf{x}}_0)+{\int }_{S+S_{\infty }}q(\mathbf{x})\frac{\partial G_0({\mathbf{x}}_0,\mathbf{x})}{\partial n}\text{d}S(\mathbf{x})={\int }_{S+S_{\infty }}\frac{\partial q(\mathbf{x})}{\partial n}{G}_0({\mathbf{x}}_0,\mathbf{x})\text{d}S(\mathbf{x}),$$

$${\int }_{S+S_{\infty }}[q(\mathbf{x})-q({\mathbf{x}}_0)]\frac{\partial G_0({\mathbf{x}}_0,\mathbf{x})}{\partial n}\text{d}S(\mathbf{x})={\int }_{S+S_{\infty }}\frac{\partial q(\mathbf{x})}{\partial n}{G}_0({\mathbf{x}}_0,\mathbf{x})\text{d}S(\mathbf{x}).$$

$$\psi (\mathbf{x})\equiv p({\mathbf{x}}_0)g(\mathbf{x})+{\left(\frac{\partial p}{\partial n}\right)}_{0}f(\mathbf{x}).$$

$${\nabla }^{2}g(\mathbf{x})+k^{2}g(\mathbf{x})=0,g({\mathbf{x}}_0)=1,\nabla g({\mathbf{x}}_0)\cdot {\mathbf{n}}_0=0,$$

$${\nabla }^{2}f(\mathbf{x})+k^{2}f(\mathbf{x})=0,f({\mathbf{x}}_0)=0,\nabla f({\mathbf{x}}_0)\cdot {\mathbf{n}}_0=1,$$

$${\int }_{S+S_{\infty }}\left[p(\mathbf{x})-p({\mathbf{x}}_0)g(\mathbf{x})-{\left(\frac{\partial p}{\partial n}\right)}_{0}f(\mathbf{x})\right]\frac{\partial G({\mathbf{x}}_0,\mathbf{x})}{\partial n}\text{d}S(\mathbf{x})$$

$$={\int }_{S+S_{\infty }}\left[\frac{\partial p(\mathbf{x})}{\partial n}-p({\mathbf{x}}_0)\nabla g(\mathbf{x})\cdot \mathbf{n}(\mathbf{x})-{\left(\frac{\partial p}{\partial n}\right)}_0\nabla f(\mathbf{x})\cdot \mathbf{n}(\mathbf{x})\right]G({\mathbf{x}}_0,\mathbf{x})\text{d}S(\mathbf{x}).$$

$$g(\mathbf{x})=\frac{a\cos [k(r_d-a)]}{r_d}+\frac{\sin [k(r_d-a)]}{k{r}_d},f(\mathbf{x})=\frac{a\sin [k(r_d-a)]}{bk{r}_d},$$

$$p({\mathbf{x}}_0){\int }_{S_{\infty }}\left[\nabla g(\mathbf{x})\cdot \mathbf{n}G-g(\mathbf{x})\frac{\partial G}{\partial n}\right]\text{d}S$$

$${\left(\frac{\partial p}{\partial n}\right)}_0{\int }_{S_{\infty }}\left[\nabla f(\mathbf{x})\cdot \mathbf{n}G-f(\mathbf{x})\frac{\partial G}{\partial n}\right]\text{d}S$$

$$4\pi p({\mathbf{x}}_p)=4\pi \psi ({\mathbf{x}}_p)-{\int }_{S+S_{\infty }}$$

$$+{\int }_{S+S_{\infty }}\frac{\partial [p(\mathbf{x})-\psi (\mathbf{x})]}{\partial n}G({\mathbf{x}}_p,\mathbf{x})\text{d}S(\mathbf{x}).$$

$$4$$

$$-{\int }_{S+S_{\infty }}\left[p(\mathbf{x})-p({\mathbf{x}}_0)g(\mathbf{x})-{\left(\frac{\partial p}{\partial n}\right)}_{0}f(\mathbf{x})\right]\left[\frac{\partial G({\mathbf{x}}_p,\mathbf{x})}{\partial n}-\frac{\partial G({\mathbf{x}}_0,\mathbf{x})}{\partial n}\right]\text{d}S(\mathbf{x})$$

$$+{\int }_{S+S_{\infty }}\left[\frac{\partial p(\mathbf{x})}{\partial n}-p({\mathbf{x}}_0)\nabla g(\mathbf{x})\cdot \mathbf{n}(\mathbf{x})-{\left(\frac{\partial p}{\partial n}\right)}_0\nabla f(\mathbf{x})\cdot \mathbf{n}(\mathbf{x})\right]\left[G({\mathbf{x}}_p,\mathbf{x})-G({\mathbf{x}}_0,\mathbf{x})\right]\text{d}S(\mathbf{x}).$$

$$r\left(\frac{\partial p}{\partial r}-ikp\right)\to 0,\text{as}r\to \infty ,$$

$${\int }_{S_{\infty }}\left(\frac{\partial g}{\partial n}G-g\frac{\partial G}{\partial n}\right)\text{ d}S$$

$${\int }_{S_{\infty }}\left(\frac{\partial f}{\partial n}G-f\frac{\partial G}{\partial n}\right)\text{ d}S$$

$$G=\frac{e^{ikr}}{r},$$

$$\underset{r\to \infty }{\lim }\frac{\partial G}{\partial n}=\underset{r\to \infty }{\lim }\frac{e^{ikr}}{r^3}(ikr-1)\mathbf{n}\cdot \mathbf{x}\to ik\frac{e^{ikr}}{r},$$

$$g=\frac{e^{ik(r_d-a)}}{2r_d}\left[a-\frac{i}{k}\right]+\frac{e^{-ik(r_d-a)}}{2r_d}\left[a+\frac{i}{k}\right],$$

$$f=\frac{-ia}{2bk{r}_d}\left[e^{ik(r_d-a)}-e^{-ik(r_d-a)}\right].$$

$$\underset{r\to \infty }{\lim }\frac{\partial g}{\partial n}\to \frac{ik}{2r_d}\left[a-\frac{i}{k}\right]e^{ik(r_d-a)}-\frac{ik}{2r_d}\left[a+\frac{i}{k}\right]e^{-ik(r_d-a)},$$

$$\underset{r\to \infty }{\lim }\frac{\partial f}{\partial n}\to \frac{a}{2b{r}_d}\left[e^{ik(r_d-a)}+e^{-ik(r_d-a)}\right]$$

$${\int }_{S_{\infty }}\left(\frac{\partial g}{\partial n}G-g\frac{\partial G}{\partial n}\right)\text{ d}S=\frac{a}{b}{\int }_{S_{\infty }}\frac{e^{ik(r-r_d+a)}}{r{r}_d}\text{ d}S,$$

$${\int }_{S_{\infty }}\left(\frac{\partial f}{\partial n}G-f\frac{\partial G}{\partial n}\right)\text{ d}S=-ika\left[1+\frac{i}{ka}\right]{\int }_{S_{\infty }}\frac{e^{ik(r-r_d+a)}}{r{r}_d}\text{ d}S.$$

$${\int }_{S_{\infty }}\frac{e^{ik(r-r_d+a)}}{r{r}_d}\text{ d}S$$

$${\int }_{S_{\infty }}\frac{e^{ik(r-r_d+a)}}{r{r}_d}\text{ d}S={\int }_{S_{\infty }}\frac{e^{ika(1+\cos \alpha )}}{r^2}\text{ d}S={\int }_0^{\pi }\frac{e^{ika(1+\cos \alpha )}}{r^2}2\pi r^2\sin \alpha \text{ d}\alpha =-\frac{2\pi }{ika}\left[1-e^{2ika}\right].$$

$${\int }_{S_{\infty }}\left(\frac{\partial g}{\partial n}{G}^p-g\frac{\partial G^p}{\partial n}\right)\text{ d}S,$$

$${\int }_{S_{\infty }}\left(\frac{\partial f}{\partial n}{G}^p-f\frac{\partial G^p}{\partial n}\right)\text{ d}S$$

$$g=\frac{e^{ik[(r_d-a_p)+(a_p-a)]}}{2r_d}\left[a-\frac{i}{k}\right]+\frac{e^{-ik[(r_d-a_p)+(a_p-a)]}}{2r_d}\left[a+\frac{i}{k}\right],$$

$$f=\frac{-ia}{2bk{r}_d}\left[e^{ik[(r_d-a_p)+(a_p-a)]}-e^{-ik[(r_d-a_p)+(a_p-a)]}\right],$$

Full text

{fmtext}

1 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 wave number 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 due 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 and 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: $\nabla^{2}\phi=0$, and the Stokes equation of fluid mechanics whereby all singular terms in the integrals are removed analytically [6]. This regularisation 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 to their characteristic dimensions or where the boundaries possess extreme geometric aspect ratios [8].

In this work, we apply this boundary regularisation 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 regularised 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-singularising 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].

2 Boundary regularised integral equation formulation (BRIEF)

To illustrate the boundary regularised integral equation formulation 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(\mathbf{x})$ outside $S$, obeys the Helmholtz equation: $\nabla^{2}p+k^{2}p=0$ where $k$ is the wave number. The solution can be found by solving the conventional boundary integral equation that follows from using Green’s second identity [10]

[TABLE]

where $G(\mathbf{x}_{0},\mathbf{x})=e^{ikr}/r$, with $r\equiv|\mathbf{x}-\mathbf{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_{\infty}$ at which the Sommerfeld radiation condition is applied (see the Appendix). The points $\mathbf{x}$ and $\mathbf{x}_{0}$ are on the boundaries. For simplicity, we assume that the surfaces $S$ and $S_{\infty}$ have well-defined tangent planes at all points on the surfaces. We refer the readers to our another work [8] for the treatment of surfaces with sharp edges and vertices. The normal derivatives are defined by $\partial p/\partial n\equiv\nabla p\cdot\mathbf{n}(\mathbf{x})$ and $\partial G/\partial n\equiv\nabla G\cdot\mathbf{n}(\mathbf{x})$ where $\mathbf{n}(\mathbf{x})$ is the unit normal vector at $\mathbf{x}$ pointing out of the solution domain. The solid angle, $c_{0}$ at $\mathbf{x}_{0}$ is equal to $2\pi$ if the tangent of the boundary at $\mathbf{x}_{0}$ is defined, otherwise it has to be calculated from the local geometry [11]. Equation (1) provides a relation between the function $p$ and its normal derivative $\partial p/\partial n$ on the surfaces $S$ and $S_{\infty}$. For Dirichlet (Neumann) problems, $p$ ($\partial p/\partial n$) is specified on the surfaces, and (1) can be solved for $\partial p/\partial n$ ($p$).

In spite of the fact that the physical problem may be well-behaved on the boundaries, the mathematical singularities in (1) at $\mathbf{x}=\mathbf{x}_{0}$ due to $G$ and $\partial G/\partial 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 (1) as follows.

But before we tackle (1) that corresponds to the Helmholtz equation for $p(\mathbf{x})$, consider first the simpler case of the Laplace equation: $\nabla^{2}q(\mathbf{x})=0$ for which the corresponding boundary integral equation is

[TABLE]

where $G_{0}(\mathbf{x}_{0},\mathbf{x})=1/|\mathbf{x}-\mathbf{x}_{0}|$. The standard way to remove the singularity arising from $\partial G_{0}/\partial n$ on the left hand side is to note that the constant $q(\mathbf{x}_{0})$ also satisfies the Laplace equation with a vanishing normal derivative on the boundaries, so that equation (2) for $[q(\mathbf{x})-q(\mathbf{x}_{0})]$ becomes

[TABLE]

This is known as the ’constant value subtraction’. However, the integral on the right hand side of (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 $\mathbf{x}\rightarrow\mathbf{x}_{0}$, $[q(\mathbf{x})-q(\mathbf{x}_{0})]$ vanishes as $|\mathbf{x}-\mathbf{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 (1) for the Helmholtz problem. However, since a constant is not a solution of the Helmholtz equation, we need to construct an auxiliary function $\psi(\mathbf{x})$ for a given value of $\mathbf{x}_{0}$ in (1) that satisfies the Helmholtz equation that is analogous to the constant, $q(\mathbf{x}_{0})$, in the case of the Laplace equation. We choose $\psi(\mathbf{x})$ to be linear in the pressure, $p(\mathbf{x}_{0})$ and its normal derivative, $(\partial p/\partial n)_{0}$ at $\mathbf{x}_{0}$

[TABLE]

Now provided the general functions $g(\mathbf{x})$ and $f(\mathbf{x})$ satisfy the Helmholtz equation with the following boundary conditions:

[TABLE]

with $\mathbf{n}_{0}\equiv\mathbf{n}(\mathbf{x}_{0})$ being the outward unit normal at $\mathbf{x}_{0}$, $\psi(\mathbf{x})$ will also satisfy the same boundary integral equation as (1). By taking the difference between the boundary integral equations for $p(\mathbf{x})$ and $\psi(\mathbf{x})$, we obtain

[TABLE]

The conditions imposed on $g(\mathbf{x})$ and $f(\mathbf{x})$ at $\mathbf{x}_{0}$ in (5) will remove both the singularities due to $G$ and $\partial G/\partial n$ in the new boundary integral equation in (2) under the generally valid assumption that as $\mathbf{x}\rightarrow\mathbf{x}_{0}$, $\left[p(\mathbf{x})-\psi(\mathbf{x})\right]$ and $\left[\partial p/\partial n-\partial\psi/\partial n\right]$ vanish as $|\mathbf{x}-\mathbf{x}_{0}|$ or faster. Therefore, if $g(\mathbf{x})$ and $f(\mathbf{x})$ satisfy (5), all singular behaviour in (2) due to $G$ and $\partial G/\partial n$ will be removed. The formal proof is a straightforward generalisation of that given in [6] for the Laplace equation.

Equation (2) is the final key result for the singularity-free [6] boundary regularised integral equation formulation (BRIEF) of the solution of the Helmholtz equation that replaces the conventional boundary integral equation in (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 $\partial p/\partial n$ (Neumann) or a relation between these two quantities on the boundary, $S$, (2) can be readily solved. Also the term involving the solid angle, $c_{0}$, in the conventional boundary integral equation, (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 (2), we see that apart from having to satisfy (5), there is considerable freedom in selecting the functional forms of $g(\mathbf{x})$ and $f(\mathbf{x})$ and hence $\psi(\mathbf{x})$. As a specific example, we can choose the following for $g(\mathbf{x})$ and $f(\mathbf{x})$

[TABLE]

where $\mathbf{x}_{d}$ is any convenient point outside the solution domain in order to ensure $\psi(\mathbf{x})$ is non-singular within the solution domain, with $a\equiv|\mathbf{x}_{0}-\mathbf{x}_{d}|$, $r_{d}\equiv|\mathbf{x}-\mathbf{x}_{d}|$ and $b\equiv(\mathbf{x}_{0}-\mathbf{x}_{d})\cdot\mathbf{n}(\mathbf{x}_{0})/|\mathbf{x}_{0}-\mathbf{x}_{d}|\neq 0$ (see figure 1a). For this particular choice of $g(\mathbf{x})$ and $f(\mathbf{x})$, the integral over the surface at infinity, $S_{\infty}$ in (2) can be found analytically using the Sommerfeld radiation condition [14, 15], (see Appendix A for details)

[TABLE]

In the limit $k=0$, the Helmholtz equation reduces to the Laplace equation with $g(\mathbf{x})\rightarrow 1$ and $f(\mathbf{x})\rightarrow(a/b)[1-a/|\mathbf{x}-\mathbf{x}_{d}|]$, and (2) will become a boundary regularised 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) 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: $\psi(\mathbf{x})$ in the case of the Helmholtz equation considered here. Apart from having to satisfy the governing differential equation and some mild constraints, see (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 $\mathbf{x}_{d}$ in (7). A broad physical interpretation is that we have constructed a pressure field, $\psi(\mathbf{x})$, that cancels the value of $p$ and $\partial p/\partial n$ at $\mathbf{x}_{0}$.

One well-studied uniqueness issue is that the solution of the boundary integral method, (1), with a closed boundary is identical to that obtained from solving directly the differential form of the Helmholtz equation: $\nabla^{2}p+k^{2}p=0$, except at a discrete set of values of $k$. These $k$ values correspond to the resonant wave numbers of the closed boundary [3, 4]. The solution of the BRIEF, (2), 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 Section 5, the higher precision that the BRIEF can attain will alleviate much of the practical numerical difficulty.

3 Accurate evaluation at field points near boundaries

The boundary regularised integral equation formulation (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 due to the near singular behaviour in the evaluation of the requisite integrals that arise in the conventional boundary integral method is often more difficult to deal with than the singularities on the boundaries.

To evaluate the solution $p(\mathbf{x}_{p})$ of the Helmholtz equation at a point $\mathbf{x}_{p}$ in the solution domain, we first use the conventional boundary integral equation to find $[p(\mathbf{x})-\psi(\mathbf{x})]$ at $\mathbf{x}_{p}$, with $\mathbf{x}_{0}$ in (4) taken to be the node on the boundary closest to $\mathbf{x}_{p}$

[TABLE]

The near singular behaviour of the surface integrals due to $G$ and $\partial G/\partial n$ as $\mathbf{x}_{p}$ approaches the boundary, $S$, can be now removed by subtracting the boundary regularised boundary integral equation (2) from (3). Then using (4) for $\psi(\mathbf{x})$ we obtain a numerically robust expression for $p(\mathbf{x}_{p})$ that is free of near singular behaviour irrespective of the distance from $\mathbf{x}_{p}$ to any boundary:

[TABLE]

4 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 conventional boundary integral method (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\omega t)$ [3], giving the boundary condition: $\partial p/\partial n=-i\omega\rho V$, with $k=\omega/c$ where $\rho$ is the density and $c$ the speed of sound in the region outside the sphere. This is a Neumann problem since $p(\mathbf{x})$ is to be found for prescribed values of the normal derivative $\partial p/\partial n$. The analytic solution of this problem is [3]: $p(r)=-[(i\omega\rho VR^{2})/(1-ikR)]e^{ik(r-R)}/r$.

In figure 1b, we compare the maximum relative error of the following six different numerical solutions of the above problem at $kR=\pi/2$ for varying numbers of degrees of freedom or unknowns on the surface of the sphere:

1. 1..

CBIM Constant: The value of the unknown pressure associated with every flat triangular surface element is assumed to be constant within the element. 2. 2..

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 (1) is calculated according to the local geometry around each point $\mathbf{x}_{0}$ [11]. 3. 3..

CBIM Linear $c_{0}=2\pi$: This an often used approach that assumes the incorrect value of $c_{0}=2\pi$ for the solid angle at every node even though the surface tangent plane is clearly not defined at the vertices of the surface elements. 4. 4..

BRIEF Constant: Using BRIEF with the assumption that the unknown pressure associated with every flat triangular surface element is constant within the element. 5. 5..

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

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 to 3 orders of magnitude; or BRIEF Quadratic can achieve the same precision as any CBIM with about $1/10th$ 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\pi$ approach is demonstrated.

In figure 2a, we compare the BRIEF using constant, linear and quadratic elements for a single radiating sphere considered at wave numbers very close to the first resonant value, $kR=\pi$ [3]. We observe that: (i) using the BRIEF with quadratic (BRIEF Quadratic) elements, the maximum relative error is not significant until $|kR-\pi|<10^{-4}$; (ii) when the BRIEF is used with constant elements (BRIEF Constant), the resonant wave number is located incorrectly. This is because the spherical shape is not well represented by constant elements and resulted in a slightly different resonant wave number, 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 to 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\pi$ 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 due to these spheres of radii $R$ and $3R$ at a distance of closest approach, $h/R=0.001$, with $kR=\pi/2$ (see figure 3). The spheres have identical time harmonic radially symmetric normal velocities as in the previous example, but are out of phase by $\pi$. In figure 3a, we show the variation of the magnitude of the pressure, $|p|$, along the longitudes of the 2 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 from 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 the magnitude of the pressure obtained using the boundary regularised integral equation formulation (BRIEF) and the conventional boundary integral method (CBIM) is given in the contour plots for two spheres of radius $R$ and $3R$, at separation $h/R=0.001$ as illustrated in figure 3b. In 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 $\rho\equiv(x^{2}+y^{2})^{1/2}\sim 0$, the BRIEF Linear results remain smooth as expected from the results in figure 3a. In contrast, the numerical errors of the conventional boundary integral method (CBIM) are already quite evident from the roughness of the pressure contours at $\rho\sim 0.7R$, and these errors become unacceptably large for $\rho<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.

5 Conclusions

We have developed the boundary regularised integral equation formulation (BRIEF) of the solution of the Helmholtz equation given by (2) that replaces the conventional boundary integral method (CBIM) given by (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 (3). 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 conventional boundary integral method 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, (2), we have assumed that the boundary $S$ is sufficient 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 wave number 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 since 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 breakdown. 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].

Acknowledgment

This work is supported in part by the Australian Research Council Discovery Project Grant Scheme to DYCC who is a Visiting Scientist at the Institute of High Performance Computing and an Adjunct Professor at the National University of Singapore.

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, $S_{\infty}$, that appear in (8), (9) and (3).

In (8) and (9) the integrals at infinity vanish due to the Sommerfeld radiation condition [14, 15]:

[TABLE]

except for the two terms

[TABLE]

and

[TABLE]

that correspond to Eqs. (8) and (9), respectively. For simplicity, we have applied $g\equiv g(\mathbf{x})$, $f\equiv f(\mathbf{x})$, $\partial{g}/\partial{n}\equiv\nabla g(\mathbf{x})\cdot\mathbf{n}(\mathbf{x})$, $\partial{f}/\partial{n}\equiv\nabla f(\mathbf{x})\cdot\mathbf{n}(\mathbf{x})$, $G\equiv G(\mathbf{x}_{0},\mathbf{x})$, $\partial{G}/\partial{n}\equiv\partial{G(\mathbf{x}_{0},\mathbf{x})}/\partial{n}$ and $\text{ d}S\equiv\text{ d}S(\mathbf{x})$.

Both $g$ and $f$ vanish as $1/r_{d}\sim 1/r$ at the surface at infinity, $S_{\infty}$, and $G$ and $\partial{G}/\partial{n}$ vanish as $1/r$, when $r\rightarrow\infty$, because

[TABLE]

and

[TABLE]

where terms that vanish faster than $1/r$ can be omitted because the surface $S_{\infty}$ only grows as $r^{2}$ at infinity. Consequently, only terms multiplying $\partial{G}/\partial{n}$, that vanish as $1/r$ can give contributions to the integrals.

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

[TABLE]

[TABLE]

Therefore the normal derivatives of $g$ and $f$ as $r\rightarrow\infty$ asymptotes to :

[TABLE]

[TABLE]

where we have again neglected higher order terms.

It can be easily seen that the terms proportional to $e^{ik(r_{d}-a)}$ in the integrals of (13) and (14) remain finite and the individual terms cancel out exactly. However, for the terms proportional to $e^{-ik(r_{d}-a)}$, this is not the case and they contribute to the integrals in (13) and (14) to give:

[TABLE]

and

[TABLE]

Finally we need to determine the integral

[TABLE]

that appears in both of the above equations. The denominator $rr_{d}$ of the integrand can be approximated by $r^{2}$ since $r\approx r_{d}$. For the numerator, we need to invoke the cosine rule and a Taylor expansion to $r-r_{d}+a$, which gives $r-r_{d}+a\approx a(1+\cos{\alpha})$ with $\alpha$ being the angle between the line segments $r_{d}$ and $a$. We then get:

[TABLE]

This leads immediately to (8) and (9).

In (3), we need to evaluate the following two integrals over the surface at infinity, $S_{\infty}$:

[TABLE]

and

[TABLE]

where $G^{p}\equiv G(\mathbf{x}_{p},\mathbf{x})$, $\partial{G^{p}}/\partial{n}\equiv\partial{G(\mathbf{x}_{p},\mathbf{x})}/\partial{n}$. Here we can rewrite $g$ and $f$ in (17) and (18) as

[TABLE]

[TABLE]

in which $a_{p}=|\mathbf{x}_{p}-\mathbf{x}_{d}|$. Following the same procedure we used to evaluate (13) and (14), the key integral that needs to be determined for calculating of (24) and (25) is

[TABLE]

where $r_{p}=|\mathbf{x}-\mathbf{x}_{p}|$. By using (23), we then find

[TABLE]