Non-Singular Burton–Miller Boundary Element Method for Acoustics

Abstract

The problem of non-unique solutions at fictitious frequencies that can appear in the boundary element method for external acoustic phenomena described by the Helmholtz equation is studied. We propose a method to fully desingularise in an analytical way the otherwise hyper-singular Burton–Miller framework, where the original boundary element method and its normal derivative are combined. The method considerably simplifies the use of higher-order elements, for example, quadratic curved surface elements. The concept is validated using the example of scattering on a rigid sphere and a rigid cube, and its robustness and effectiveness for external sound-wave problems are confirmed.

1. Introduction

Understanding acoustic phenomena [1] is of great importance for transitional research and practical applications, ranging from musical instruments [2,3], traffic, or airplane noise prediction [4,5], to sensing and monitoring for structural health [6]. In fluids, acoustic phenomena manifest themselves as the temporal–spatial distribution and evolution of pressure or density waves.

One of the most common acoustic phenomena is sound waves [7] travelling in a fluid at rest, where viscous effects are negligible. Furthermore, assume that the sound wave introduces small perturbations to the fluid density, pressure, and velocity, such that higher-order effects can be ignored. The continuity and momentum equations for linearised acoustic phenomena in a fluid can then be obtained. Meanwhile, under the linearised assumption, the disturbed or fluctuated fluid movement associated with the sound wave is irrotational [8], we can introduce a velocity potential, and the continuity and momentum equations are combined and simplified to become a linear wave equation [9]. The wave equation describing linear acoustics is ${\nabla }^2{\varphi }^{\prime }=1/c^2{\partial }^2{\varphi }^{\prime }/\partial t^2$, with ${\nabla }^2$ the Laplacian operator, c the speed of sound, t the time, and ${\varphi }^{\prime }$ the velocity potential. It is usually convenient and intuitive to analyse periodic wave phenomena in the frequency domain, and we have powerful tools based on the Fourier transform to convert wave problems between the time and frequency domains. In the frequency domain, with ${\varphi }^{\prime }(\mathit{x},t)=\varphi \left(\mathit{x}\right)exp(-\mathrm{i}\omega t)$, $\omega$ being the angular frequency and $\mathrm{i}=\sqrt{-1}$ the imaginary unit, the wave equation turns into the scalar Helmholtz equation,

$${\nabla }^2\varphi +k^2\varphi =0,$$

(1)

where $k=\omega /c$ is the wave number. The pressure, p, can be obtained from the velocity potential as $p=\mathrm{i}\omega {\rho }_0\varphi$ in the frequency domain, in which ${\rho }_0$ is the reference fluid density.

To study sound propagation and scattering in a fluid, we need to solve the Helmholtz equation in Equation (1). Except for a few special cases, such as a spherical or a cylindrical scatterer, it is not feasible to obtain an analytical solution of the Helmholtz equation. Thus, we have to revert to solving the Helmholtz equation numerically and the boundary element method is an effective way to do so [10,11,12,13,14]. The Helmholtz equation (1) is elliptic in nature. Based on Gauss’ theorem and Green’s identities, the boundary element method uses this ellipticity to obtain a relationship between the velocity potential and its normal derivative on the surface of any object or scatterer. Thus, the dimension of the problem can be reduced effectively from a fully three-dimensional (3D) problem to a two-dimensional (2D) surface problem. In addition, the boundary element method is especially advantageous for open space (external) problems, as the Sommerfeld radiation condition [15] at infinity is automatically satisfied and only the surfaces of the objects, S, need to be considered in the meshing. The boundary element method is thus especially suited to solve external problems.

However, a main drawback of external problems when using the boundary element method with large ka values (with a a typical dimension of the problem) is that for certain values of ka, the so-called fictitious wave numbers (or frequencies), the internal resonance solution of the object can appear in the outer solution and, thereby, will produce non-unique (spurious) results. Although it has been shown that it is possible to eliminate certain fictitious frequencies using a modified Green’s function [16], unfortunately, not all fictitious modes can be eliminated in this manner. We have to revert to one of the two most commonly used methods, namely the CHIEF method of Schenck [17] or the method of Burton and Miller [18]. The CHIEF method adds interior points in the internal domain in order to force the internal solution to be zero. The CHIEF method is not entirely rigorous since it is not clearly defined how many and where the internal points should be. In addition, it leads to more equations than variables and a least-square minimisation scheme has to be employed, which can be quite time consuming. The Burton–Miller method, on the other hand, takes the normal derivative of the boundary element equation and adds it to the original boundary element equation. The basic idea is that these two equations always have different fictitious frequencies and it can be shown that, in theory, the combination of these two equations should no longer possess any spurious solutions [19]. The problem with this method is that the integrands now become hypersingular and have to be treated with extreme care [20].

In this paper, the scientific motivation and novelty are to fully desingularise the Burton–Miller formulation, which originally has integrals exhibiting hypersingular integrands. The thus-developed desingularised framework simplifies the numerical implementation considerably since only standard numerical integration schemes are needed. This can greatly increase the opportunity and significantly reduce the threshold to use a technically difficult surface method, the Burton–Miller boundary element method, to solve external acoustic problems robustly, in particular, if one would like to use higher-order surface elements to improve the computational efficiency and accuracy. We demonstrate our method in Section 2 and then validate it in Section 3 using the examples of a rigid sphere and a rigid cube, followed by the conclusion in Section 4.

2. The Burton–Miller Framework

2.1. Overview

In this section, we briefly outline the theory of the following sections. We start with the classical boundary integral equation for the Helmholtz equation in Equation (2). The singular integrals can be rewritten as non-singular integrals by subtracting a function satisfying the Laplace equation and using the corresponding Laplace boundary integral method, as demonstrated in Equation (4) in Section 2.2.

In the Burton–Miller framework, we need the normal derivative of the standard Helmholtz boundary integral equation, as in Equation (5). This equation contains hypersingular integrals and must be treated very carefully. It is, however, possible to desingularise it by realising that subtracting the hypersingular kernels of the Helmholtz and Laplace equations results in a weakly singular integrand (see Equation (11)) and using this result in the non-singular Equation (14). This is described in Section 2.3.

Finally, the Burton–Miller idea is to combine the standard boundary integral equation with the normal derivative variant by multiplying by a complex-valued constant $\beta \mathrm{i}$, as shown in Equation (15). It turns out that we also need to solve an additional boundary integral Equation (16) for a Laplace equation to close the system of equations.

2.2. The Standard Boundary Integral Equation

Given that the Helmholtz equation is elliptic, if the potential and its normal derivative on the surface of the object are given as illustrated in Figure 1, we can obtain all the properties of the acoustic and flow phenomena in the entire fluid domain. However, we normally only know the potential or the normal derivative of the potential (or a linear combination of them) on the surface of the object. To obtain the missing part (see Figure 1 for an example), we can solve the boundary integral equation corresponding to the Helmholtz equation

$$c\left({\mathit{x}}_0\right)\varphi \left({\mathit{x}}_0\right)+{\int }_S\varphi \left(\mathit{x}\right)\frac{\partial G_k}{\partial n}\mathrm{d}S\left(\mathit{x}\right)={\int }_S\frac{\partial \varphi \left(\mathit{x}\right)}{\partial n}{G}_k\mathrm{d}S\left(\mathit{x}\right),$$

(2)

in which $S\left(\mathit{x}\right)$ is the surface of the object under consideration, ${\mathit{x}}_0$ is the observation point (located on the surface S), $c\left({\mathit{x}}_0\right)$ is the solid angle at ${\mathit{x}}_0$, $\mathit{x}$ is the integration point, and $G_k\equiv G_k({\mathit{x}}_0,\mathit{x})=exp\left(\mathrm{i}kr\right)/r$ is the Green’s function for the Helmholtz equation, with $r=|\mathit{x}-{\mathit{x}}_0|$. The symbol $\partial /\partial n=\mathit{n}·\nabla$ represents the normal derivative, where $\mathit{n}$ is the unit normal direction at $\mathit{x}$ pointing out of the domain (thus, into the object).

Suppose we construct a function in the domain of interest as $\psi \left(\mathit{x}\right)=\varphi \left({\mathit{x}}_0\right)+{(\partial \varphi /\partial n)}_0{\mathit{n}}_0·(\mathit{x}-{\mathit{x}}_0)$. Here, $\varphi \left({\mathit{x}}_0\right)$ is the velocity potential and ${(\partial \varphi /\partial n)}_0$ is the normal derivative of that potential at ${\mathit{x}}_0$, which are both constants at a given ${\mathit{x}}_0$, and ${\mathit{n}}_0\equiv \mathit{n}\left({\mathit{x}}_0\right)$ is the unit normal vector at ${\mathit{x}}_0$ on S pointing out of the domain. Clearly, $\psi \left(\mathit{x}\right)$ is a linear function of $\mathit{x}$ that satisfies the Laplace equation as ${\nabla }^2\psi \left(\mathit{x}\right)=0$. The Laplace equation can be taken as a special case of the Helmholtz equation when k = 0, and we can write a similar boundary integral equation for $\psi$ similar to Equation (2) as

$$\begin{array}{cc}\hfill & [c\left({\mathit{x}}_0\right)-4\pi ]\varphi \left({\mathit{x}}_0\right)+{\int }_S\left[\varphi \left({\mathit{x}}_0\right)+{\left(\frac{\partial \varphi }{\partial n}\right)}_0\mathit{n}·(\mathit{x}-{\mathit{x}}_0)\right]\frac{\partial G_0}{\partial n}\mathrm{d}S\left(\mathit{x}\right)\hfill \\ \hfill =& {\int }_S{\left(\frac{\partial \varphi }{\partial n}\right)}_0({\mathit{n}}_0·\mathit{n})G_0\mathrm{d}S\left(\mathit{x}\right),\hfill \end{array}$$

(3)

in which $G_0\equiv G_0({\mathit{x}}_0,\mathit{x})=1/r$ is the fundamental solution of the Laplace equation, the normal derivative $\partial \psi /\partial n={(\partial \varphi /\partial n)}_0({\mathit{n}}_0·\mathit{n})$ on the surface S has been used, and the term with $-4\pi$ is the contribution from the surface at infinity [21,22,23]. Noting that $({\mathit{n}}_0·\mathit{n})=1$ in the limit of $\mathit{x}\to {\mathit{x}}_0$ and subtracting Equation (3) from Equation (2), we obtain [24]:

$$\begin{array}{cc}\hfill & 4\pi \varphi \left({\mathit{x}}_0\right)+{\int }_S\left[\varphi \left(\mathit{x}\right)\frac{\partial G_k}{\partial n}-\varphi \left({\mathit{x}}_0\right)\frac{\partial G_0}{\partial n}\right]\mathrm{d}S\left(\mathit{x}\right)\hfill \\ \hfill =& {\int }_S{\left(\frac{\partial \varphi }{\partial n}\right)}_0{\mathit{n}}_0·(\mathit{x}-{\mathit{x}}_0)\frac{\partial G_0}{\partial n}\mathrm{d}S\left(\mathit{x}\right)+{\int }_S\left[\frac{\partial \varphi \left(\mathit{x}\right)}{\partial n}{G}_k-{\left(\frac{\partial \varphi }{\partial n}\right)}_0({\mathit{n}}_0·\mathit{n})G_0\right]\mathrm{d}S\left(\mathit{x}\right).\hfill \end{array}$$

(4)

The integrands in Equation (4) are now fully regularised.

2.3. The Normal Derivative of the Boundary Integral Equation

We can also perform the normal derivative ${[\partial (·)/\partial n]}_0\equiv \partial (·)/\partial n_0\equiv {\mathit{n}}_0·{\nabla }_{{\mathit{x}}_0}(·)$ operator on Equation (2), which leads to

$$c\left({\mathit{x}}_0\right){\left(\frac{\partial \varphi }{\partial n}\right)}_0+{\int }_S\varphi \left(\mathit{x}\right)\frac{{\partial }^2{G}_k}{\partial n\partial n_0}\mathrm{d}S\left(\mathit{x}\right)={\int }_S\frac{\partial \varphi \left(\mathit{x}\right)}{\partial n}\frac{\partial G_k}{\partial n_0}\mathrm{d}S\left(\mathit{x}\right).$$

(5)

The integral on the left is now hypersingular and must be treated with the utmost care when solved directly. In addition, it is worth noting that to obtain the above equation, we have implicitly used ${[\partial c\left({\mathit{x}}_0\right)/\partial n]}_0=0$. This is because

$$c\left({\mathit{x}}_0\right)+{\int }_S\frac{\partial G_0}{\partial n}\mathrm{d}S\left(\mathit{x}\right)=4\pi .$$

(6)

Performing the ${[\partial (·)/\partial n]}_0$ operator on the above Equation (6), we have

$${\left[\frac{\partial c\left({\mathit{x}}_0\right)}{\partial n}\right]}_0+{\int }_S\frac{{\partial }^2{G}_0}{\partial n\partial n_0}\mathrm{d}S\left(\mathit{x}\right)=0.$$

(7)

Since it has been shown [25,26] that

$$\begin{array}{c}\hfill {\int }_S\frac{{\partial }^2{G}_0}{\partial n\partial n_0}\mathrm{d}S\left(\mathit{x}\right)=0,\end{array}$$

(8)

we conclude that ${[\partial c\left({\mathit{x}}_0\right)/\partial n]}_0=0$.

Either Equation (2) or Equation (5) can be used to solve the Helmholtz equation in Equation (1). As we mentioned before, for external problems, both of them suffer from spurious solutions at fictitious frequencies.

Although the traditional Burton–Miller boundary element method based on Equations (2) and (5) eliminates the spurious solutions for external problems, there remain challenges in its numerical implementation due to the singularities as $\mathit{x}\to {\mathit{x}}_0$ when using linear or high-order surface elements, in particular, the hypersingularity. In addition, if we would like to use higher-order surface elements to represent the object, the calculation of the solid angle $c\left({\mathit{x}}_0\right)$ is tedious. As such, it would be desirable to have a Burton–Miller boundary element method for the Helmholtz equation in which all the singularities and solid angles are fully removed analytically, as demonstrated below.

To handle the main challenge of the integrand on the left-hand side of Equation (5), which is hypersingular as $\mathit{x}\to {\mathit{x}}_0$, let us define a function ${\psi }_1$ that satisfies the Laplace equation instead of the Helmholtz equation and is applied to the same domain and boundaries. Following the same procedure used to obtain Equation (5), we have

$$\begin{array}{c}\hfill c\left({\mathit{x}}_0\right){\left(\frac{\partial {\psi }_1}{\partial n}\right)}_0+{\int }_S{\psi }_1\left(\mathit{x}\right)\frac{{\partial }^2{G}_0}{\partial n\partial n_0}\mathrm{d}S\left(\mathit{x}\right)={\int }_S\frac{\partial {\psi }_1\left(\mathit{x}\right)}{\partial n}\frac{\partial G_0}{\partial n_0}\mathrm{d}S\left(\mathit{x}\right).\end{array}$$

(9)

Assuming now that ${\psi }_1=\varphi$ on the surface S (and only on this surface), subtract Equation (9) from Equation (5), to obtain

$$\begin{array}{cc}\hfill & c\left({\mathit{x}}_0\right){\left(\frac{\partial \varphi }{\partial n}\right)}_0-c\left({\mathit{x}}_0\right){\left(\frac{\partial {\psi }_1}{\partial n}\right)}_0+{\int }_S\varphi \left(\mathit{x}\right)\left[\frac{{\partial }^2{G}_k}{\partial n\partial n_0}-\frac{{\partial }^2{G}_0}{\partial n\partial n_0}\right]\mathrm{d}S\left(\mathit{x}\right)\hfill \\ \hfill =& {\int }_S\frac{\partial \varphi \left(\mathit{x}\right)}{\partial n}\frac{\partial G_k}{\partial n_0}\mathrm{d}S\left(\mathit{x}\right)-{\int }_S\frac{\partial {\psi }_1\left(\mathit{x}\right)}{\partial n}\frac{\partial G_0}{\partial n_0}\mathrm{d}S\left(\mathit{x}\right).\hfill \end{array}$$

(10)

The integrands in Equation (10) are all weakly singular [27] since (see Appendix A)

$$\underset{\mathit{x}\to {\mathit{x}}_0}{lim}\left(\frac{{\partial }^2{G}_k}{\partial n\partial n_0}-\frac{{\partial }^2{G}_0}{\partial n\partial n_0}\right)=\frac{k^2}{2|\mathit{x}-{\mathit{x}}_0|}.$$

(11)

To analytically remove the remaining singularities and the terms associated with the solid angle at ${\mathit{x}}_0$ in Equation (10), we set up the following three functions that all satisfy the Laplace equation:

$$\begin{array}{cc}\hfill {\psi }_2\left(\mathit{x}\right)& =\frac{k^2}{2}\varphi \left({\mathit{x}}_0\right){\mathit{n}}_0·(\mathit{x}-{\mathit{x}}_0)\left(\mathrm{linear}\right),\hfill \\ \hfill {\psi }_3\left(\mathit{x}\right)& ={\left(\frac{\partial \varphi }{\partial n}\right)}_0\left(\mathrm{constant}\right),\hfill \\ \hfill {\psi }_4\left(\mathit{x}\right)& =-{\left(\frac{\partial {\psi }_1}{\partial n}\right)}_0\left(\mathrm{constant}\right).\hfill \end{array}$$

(12)

Note that ${\psi }_2\left(\mathit{x}\right)$ is a linear function and ${\psi }_3\left(\mathit{x}\right)$ and ${\psi }_4\left(\mathit{x}\right)$ are constants in the domain, as $\varphi \left({\mathit{x}}_0\right)$, ${(\partial \varphi /\partial n)}_0$, and ${(\partial {\psi }_1/\partial n)}_0$ are constants for a given ${\mathit{x}}_0$. Following the same procedure used to obtain Equation (3), we can write the boundary integral equations for ${\psi }_2\left(\mathit{x}\right),{\psi }_3\left(\mathit{x}\right),{\psi }_4\left(\mathit{x}\right)$, respectively, as

$$\begin{array}{cc}\hfill & {\int }_S\frac{k^2}{2}\varphi \left({\mathit{x}}_0\right){\mathit{n}}_0·(\mathit{x}-{\mathit{x}}_0)\frac{\partial G_0}{\partial n}\mathrm{d}S\left(\mathit{x}\right)={\int }_S\frac{k^2}{2}\varphi \left({\mathit{x}}_0\right)({\mathit{n}}_0·\mathit{n})G_0\mathrm{d}S\left(\mathit{x}\right),\hfill \\ \hfill & [c\left({\mathit{x}}_0\right)-4\pi ]{\left(\frac{\partial \varphi }{\partial n}\right)}_0+{\int }_S{\left(\frac{\partial \varphi }{\partial n}\right)}_0\frac{\partial G_0}{\partial n}\mathrm{d}S\left(\mathit{x}\right)=0,\hfill \\ \hfill & [4\pi -c\left({\mathit{x}}_0\right)]{\left(\frac{\partial {\psi }_1}{\partial n}\right)}_0-{\int }_S{\left(\frac{\partial {\psi }_1}{\partial n}\right)}_0\frac{\partial G_0}{\partial n}\mathrm{d}S\left(\mathit{x}\right)=0,\hfill \end{array}$$

(13)

in which, $\partial {\psi }_2\left(\mathit{x}\right)/\partial n=k^2\varphi \left({\mathit{x}}_0\right)\left({\mathit{n}}_0·\mathit{n}\right)/2$, $\partial {\psi }_3\left(\mathit{x}\right)/\partial n=0$, and $\partial {\psi }_4\left(\mathit{x}\right)/\partial n=0$ have been used. Subtracting Equation (13) from Equation (10), we obtain

$$\begin{array}{cc}\hfill & {\int }_S\left\{\varphi \left(\mathit{x}\right)\left[\frac{{\partial }^2{G}_k}{\partial n\partial n_0}-\frac{{\partial }^2{G}_0}{\partial n\partial n_0}\right]-\frac{k^2}{2}\varphi \left({\mathit{x}}_0\right)\left({\mathit{n}}_{\mathbf{0}}·\mathit{n}\right)G_0\right\}\mathrm{d}S\left(\mathit{x}\right)\hfill \\ \hfill =& {\int }_S\left[\frac{\partial \varphi \left(\mathit{x}\right)}{\partial n}\frac{\partial G_k}{\partial n_0}+{\left(\frac{\partial \varphi }{\partial n}\right)}_0\frac{\partial G_0}{\partial n}\right]\mathrm{d}S\left(\mathit{x}\right)-{\int }_S\left[\frac{\partial {\psi }_1\left(\mathit{x}\right)}{\partial n}\frac{\partial G_0}{\partial n_0}+{\left(\frac{\partial {\psi }_1}{\partial n}\right)}_0\frac{\partial G_0}{\partial n}\right]\mathrm{d}S\left(\mathit{x}\right)\hfill \\ \hfill & -{\int }_S\frac{k^2}{2}\varphi \left({\mathit{x}}_0\right){\mathit{n}}_0·(\mathit{x}-{\mathit{x}}_0)\frac{\partial G_0}{\partial n}\mathrm{d}S\left(\mathit{x}\right)-4\pi {\left(\frac{\partial \varphi }{\partial n}\right)}_0+4\pi {\left(\frac{\partial {\psi }_1}{\partial n}\right)}_0.\hfill \end{array}$$

(14)

The integrands in Equation (14) are now all regular and this equation is fully desingularised. Note that two plus signs appear in this equation since differentiation of $r=|\mathit{x}-{\mathit{x}}_0|$ with respect to ${\mathit{x}}_0$ will produce a minus sign compared to $\mathit{x}$, that is, ${\nabla }_{x_0}\left(r\right)=-{\nabla }_x\left(r\right)$. Here, it is worth noting that we defined $\partial (·)/\partial n={\nabla }_{\mathit{x}}(·)·\mathit{n}\left(\mathit{x}\right)$ and ${[\partial (·)/\partial n]}_0={\nabla }_{{\mathit{x}}_0}(·)·\mathit{n}\left({\mathit{x}}_0\right)$. For the acoustic flow velocity potential, $\partial \varphi /\partial n={\nabla }_{\mathit{x}}\varphi ·\mathit{n}\left(\mathit{x}\right)$ and ${(\partial \varphi /\partial n)}_0={\nabla }_{{\mathit{x}}_0}\varphi ·\mathit{n}\left({\mathit{x}}_0\right)$ refer to the normal velocity at $\mathit{x}$ and ${\mathit{x}}_0$ on the boundary, respectively. They are just the normal velocity at different locations by definition such that no minus sign difference will appear here. However, $G_0=1/r=1/|\mathit{x}-{\mathit{x}}_0|$ is the kernel that reflects how some physical property at $\mathit{x}$ will affect that at ${\mathit{x}}_0$ or vice versa. Thus, its derivative will depend on the location ($\mathit{x}$ or ${\mathit{x}}_0$) we choose as the observation point and the location we choose as the computation point. In this case, there will be a minus sign difference between $\partial G_0/\partial n$ and $\partial G_0/\partial n_0$ with respect to $\mathit{x}$ and ${\mathit{x}}_0$, respectively, since $r=|\mathit{x}-{\mathit{x}}_0|=\sqrt{{(x-x_0)}^2+{(y-y_0)}^2+{(z-z_0)}^2}$ and ${\nabla }_{x_0}\left(r\right)=-{\nabla }_x\left(r\right)=-(\mathit{x}-{\mathit{x}}_0)/r$.

2.4. The Non-Singular Burton–Miller Formulation

To recover the uniqueness of the solution of the boundary integral method for the Helmholtz equation in the external domain, Burton and Miller proposed to combine Equations (2) and (5) as {Equation (2) + $\mathrm{i}\beta$ Equation (5)} since the fictitious frequencies in Equations (2) and (5) are always different from each other. To balance Equations (2) and (5) in dimension (dimensional homogeneous in length), the parameter $\beta$ should be related to a characterised length of the problem under consideration. The obvious choices are the size of the scattering object or the inverse of the wave number k.

Following the Burton–Miller concept and combining the non-singular boundary integral equations in Equations (4) and (14), we obtain the non-singular Burton–Miller boundary integral method for the Helmholtz equation as

$$\begin{array}{cc}\hfill & 4\pi \varphi \left({\mathit{x}}_0\right)+{\int }_S\left[\varphi \left(\mathit{x}\right)\frac{\partial G_k}{\partial n}-\varphi \left({\mathit{x}}_0\right)\frac{\partial G_0}{\partial n}\right]\mathrm{d}S\left(\mathit{x}\right)\hfill \\ \hfill & +\mathrm{i}\beta {\int }_S\left\{\varphi \left(\mathit{x}\right)\left[\frac{{\partial }^2{G}_k}{\partial n\partial n_0}-\frac{{\partial }^2{G}_0}{\partial n\partial n_0}\right]-\frac{k^2}{2}\varphi \left({\mathit{x}}_0\right)\left({\mathit{n}}_0·\mathit{n}\right)G_0\right\}\mathrm{d}S\left(\mathit{x}\right)\hfill \\ \hfill & ={\int }_S\left[\frac{\partial \varphi \left(\mathit{x}\right)}{\partial n}{G}_k-{\left(\frac{\partial \varphi }{\partial n}\right)}_0\left({\mathit{n}}_0·\mathit{n}\right)G_0\right]\mathrm{d}S\left(\mathit{x}\right)\hfill \\ \hfill & +{\int }_S{\left(\frac{\partial \varphi }{\partial n}\right)}_0{\mathit{n}}_0·(\mathit{x}-{\mathit{x}}_0)\frac{\partial G_0}{\partial n}\mathrm{d}S\left(\mathit{x}\right)\hfill \\ \hfill & +\mathrm{i}\beta {\int }_S\left[\frac{\partial \varphi \left(\mathit{x}\right)}{\partial n}\frac{\partial G_k}{\partial n_0}+{\left(\frac{\partial \varphi }{\partial n}\right)}_0\frac{\partial G_0}{\partial n}\right]\mathrm{d}S\left(\mathit{x}\right)\hfill \\ \hfill & -\mathrm{i}\beta {\int }_S\left[\frac{\partial {\psi }_1\left(\mathit{x}\right)}{\partial n}\frac{\partial G_0}{\partial n_0}+{\left(\frac{\partial {\psi }_1}{\partial n}\right)}_0\frac{\partial G_0}{\partial n}\right]\mathrm{d}S\left(\mathit{x}\right)\hfill \\ \hfill & -\mathrm{i}\beta {\int }_S\frac{k^2}{2}\varphi \left({\mathit{x}}_0\right){\mathit{n}}_0·(\mathit{x}-{\mathit{x}}_0)\frac{\partial G_0}{\partial n}\mathrm{d}S\left(\mathit{x}\right)\hfill \\ \hfill & -\mathrm{i}\beta 4\pi {\left(\frac{\partial \varphi }{\partial n}\right)}_0+\mathrm{i}\beta 4\pi {\left(\frac{\partial {\psi }_1}{\partial n}\right)}_0,\hfill \end{array}$$

(15)

which should be solved simultaneously with the desingularised version of Equation (9):

$$\begin{array}{cc}\hfill & 4\pi \varphi \left({\mathit{x}}_0\right)+{\int }_S\left[\varphi \left(\mathit{x}\right)-\varphi \left({\mathit{x}}_0\right)\right]\frac{\partial G_0}{\partial n}\mathrm{d}S\left(\mathit{x}\right)\hfill \\ \hfill =& {\int }_S\left[\frac{\partial {\psi }_1\left(\mathit{x}\right)}{\partial n}-{\left(\frac{\partial {\psi }_1}{\partial n}\right)}_0\left({\mathit{n}}_0·\mathit{n}\right)\right]G_0\mathrm{d}S\left(\mathit{x}\right)+{\int }_S{\left(\frac{\partial {\psi }_1}{\partial n}\right)}_0{\mathit{n}}_0·(\mathit{x}-{\mathit{x}}_0)\frac{\partial G_0}{\partial n}\mathrm{d}S\left(\mathit{x}\right).\hfill \end{array}$$

(16)

Equation (16) is the non-singular boundary integral formulation similar to Equation (4) to link $\partial {\psi }_1/\partial n$ and ${\psi }_1=\varphi$ on the surface S.

Equations (15) and (16) are the non-singular versions of the Burton–Miller boundary integrals (which turn into the boundary element method once discretised) to solve the Helmholtz equation for external problems in which all the integrands are regular and all the solid angles are removed analytically. In the numerical procedure, only standard numerical integration schemes, such as Gauss quadrature, are needed to calculate all integrals, which significantly simplifies the numerical implementation.

3. Results

3.1. Scattering from a Rigid Sphere

We test our framework of the non-singular Burton–Miller boundary element method given in Equations (15) and (16) via the classical example of the scattered sound wave by a rigid sphere under an incoming plane wave, as illustrated in Figure 2a. The radius of the rigid sphere is a, and its centre coincides with the origin of the coordinate system. The incoming plane wave moves along the positive z direction with wave number k and amplitude ${\varphi }_0$, whose potential profile can be written as ${\varphi }^{\mathrm{inc}}={\varphi }_{0}exp\left(\mathrm{i}kz\right)$. On the surface of the rigid sphere, the boundary condition is that the normal flow velocity is zero, which leads to $\partial {\varphi }^{\mathrm{inc}}/\partial n+\partial {\varphi }^{\mathrm{sc}}/\partial n=0$, where ${\varphi }^{\mathrm{sc}}$ is the potential of the scattered wave.

The scattered sound-wave potential, ${\varphi }^{\mathrm{sc}}$, is calculated by both the non-singular standard boundary element method, as in Equation (4), and the non-singular Burton–Miller boundary element method, as in Equations (15) and (16), with $\beta =1/k$ when ka spans from 1 to 40 with a step of 0.001. The rigid sphere surface is represented by 5762 nodes connected by 2880 higher-order quadratic triangular elements for all ka values.

We choose the scattered velocity potential at point (0, 0, 1.5a) to serve as the test location to compare the results obtained by those two methods. The procedure to obtain the potential at that location is to solve the scattered potential on the surface of the sphere and obtain the potential in the domain through post-processing [28]. In Figure 2b–d, we show the real part, the imaginary part, and the absolute value of the velocity potential of the scattered wave at the location (0, 0, 1.5a). It is clear that the results obtained by the non-singular standard boundary element method given in Equation (4) are not reliable for high ka values. Non-unique or spurious solutions at fictitious frequencies appear more and more as the value of ka increases, as shown by the spiky black curves in Figure 2b–d. On the contrary, with the non-singular Burton–Miller boundary element method given in Equations (15) and (16), no spurious solutions appear as ka is scanned with a tiny step of 0.001 from ka = 1 up to a high value of 40, as shown by the red curves in Figure 2b–d. Clearly, the non-physical solutions at the classical internal resonance frequencies show up in the curve for the non-singular standard boundary element method, for example, at ka = 9.356 and ka = 10.417, corresponding to the theoretical values of ka = 9.3558 and ka = 10.41712, respectively (see also Table 1 in Ref. [16]).

We also compare the results obtained from the non-singular standard boundary element method of Equation (4) and the non-singular Burton–Miller boundary element method in Equations (15) and (16) with the analytical solution [29] for larger ka ranging from 22 to 40. As shown in Figure 3, the results from the non-singular standard boundary element method are heavily polluted by the spurious solutions, whereas those from the non-singular Burton–Miller boundary element method are in good agreement with the analytical solution for large ka, using a moderate number of surface nodes and higher-order quadratic surface elements.

3.2. Scattering from a Rigid Cube

After having performed the verification of a rigid sphere in the previous section, we present a second example of scattering on a rigid cube. As shown in Figure 4a, the cube is aligned with the x-, y-, and z-axes and its centre is situated at (0, 0, 0). The sides of the cube have a length of 2a each and the corners of the cube are rounded in order to facilitate the use of quadratic elements. An incoming wave with wavenumber k travels in the positive z-direction and hits the cube. As in the sphere example, we have placed an observation point at (x, y, z) = (0, 0, 1.5a), thus at a distance of 0.5a behind the cube. In Figure 4b, the normalised real part of the scattered potential ${\varphi }^{rmsc}$ is plotted as a function of ka, where ka ranges from 1 to 40. The graph was generated by increasing ka in steps of 0.001. As in the previous example, the scattered sound-wave potential, ${\varphi }^{\mathrm{sc}}$, is obtained with the non-singular standard boundary element method from Equation (4) and the non-singular Burton–Miller boundary element method from Equations (15) and (16), with $\beta =\mathrm{i}/k$. The cube surface has 11,642 nodes connected by 5820 high-order quadratic triangular elements. As shown in Figure 4b, the red line (Equations (15) and (16)) and the black line (Equation (4)) overlap for up to around ka = 10. For larger values of ka, spurious solutions occur, which are manifested as sharp peaks. A similar trend can be seen in Figure 3c, where the imaginary part of ${\varphi }^{\mathrm{sc}}$ is plotted as a function of ka. It is interesting to note that the peaks that are clearly visible in Figure 4b are not always seen in Figure 4c, and vice versa. Therefore, the absolute value of ${\varphi }^{\mathrm{sc}}$ is also plotted as a function of ka in Figure 4d. The error in the solution for the non-singular standard boundary element method appears to increase with increasing ka, whereas it remains only a few percent for lower ka. However, it can quickly become several times larger, as indicated by the height of the black spikes in Figure 4d.

The observed fictitious frequencies in Figure 4 correspond very closely to the theoretical room-mode formula [29]:

$$f_{room}=\frac{c}{2}\sqrt{\frac{p^2}{L^2}+\frac{q^2}{W^2}+\frac{r^2}{H^2}}$$

(17)

where $f_{room}$ is the frequency where resonance appears in a room with length L, width W, and height H, and p, q, and r are integers 0,1,2 …… Since in our case of a cube, L = 2a, W = 2a, and H = 2a and the relationship between $f_{room}$ and k is $f_{room}=kc/\left(2\pi \right)$, we predict that the non-physical solutions might happen at

$$ka=\frac{\pi }{2}\sqrt{p^2+q^2+r^2}.$$

(18)

For example, the observed numerical non-physical solution peak at ka = 5.878 corresponds to p = 3, q = 2, and r = 1, with ka = 5.87738, and the numerical non-physical solution peak at ka = 9.683 corresponds to p = 5, q = 3, and r = 2, with ka = 9.68304. It appears that, even with a very small interval of 0.001 for ka, we still missed many of the fictitious frequencies. This confirms the findings of [16,28], where it was found that for moderate ka values, the non-singular boundary element method is already quite good at avoiding non-physical solutions at fictitious frequencies.

The validation examples demonstrated here have shown the robustness and effectiveness of our non-singular Burton–Miller boundary element framework to solve acoustic scattering.

4. Conclusions

In this article, it is shown that it is possible to write the Burton–Miller integral equations in a totally non-singular form. Thus, all elements of the structure can be integrated by standard Gaussian quadrature integration formulas, thereby considerably simplifying the numerical implementation. This ensures that higher-order quadratic meshes can be used with ease to achieve better computational efficiency and accuracy. A further advantage is that solid angles are no longer present in the final integral equations. Such advantages can significantly reduce the threshold to use a technically difficult surface method, the Burton–Miller boundary element method, to solve external acoustic problems robustly and effectively. A slight drawback is that a $2N\times 2N$ system must now be solved instead of an $N\times N$ system. On the other hand, the use of quadratic elements does allow us to considerably reduce the number of nodes N to achieve the same accuracy compared to flat or linear elements [28]. The concept was tested with two examples, scattering on a rigid sphere and scattering on a rigid cube. The results from the conventional (desingularised) boundary element method and the desingularised Burton–Miller framework were compared for these two examples.

In this paper, we addressed the fictitious frequency issue for sound-wave phenomena described by a scalar Helmholtz equation. However, such issues also appear in vector wave problems, for example, electromagnetic scattering, which in essence, can be considered as three coupled scalar Helmholtz equations—one for each Cartesian component of the electric field [30,31,32,33,34,35] and possibly even for elastic waves [36]. Our framework demonstrated in this paper can also eliminate the non-unique or non-physical solutions at fictitious frequencies in vector wave problems. Furthermore, using Fourier transforms, we can easily extend the method to pulsed phenomena [37]. As such, our framework can be applied to effectively and robustly solve wave phenomena in multi-disciplinary physics.