Eliminating the fictitious frequency problem in BEM solutions of the external Helmholtz equation

TL;DR

This paper addresses the fictitious frequency problem in boundary element method solutions of the exterior Helmholtz equation, proposing a modified Green's function and a desingularized BEM to detect and avoid these issues.

Contribution

It introduces a new approach using a modified Green's function and a fully desingularized BEM to eliminate fictitious frequencies in Helmholtz problem solutions.

Findings

Modified Green's function reduces fictitious frequencies.

Desingularized BEM improves numerical accuracy.

Method effectively detects and avoids fictitious solutions.

Abstract

The problem of the fictitious frequency spectrum resulting from numerical implementations of the boundary element method for the exterior Helmholtz problem is revisited. When the ordinary 3D free space Green's function is replaced by a modified Green's function, it is shown that these fictitious frequencies do not necessarily have to correspond to the internal resonance frequency of the object. Together with a recently developed fully desingularized boundary element method that confers superior numerical accuracy, a simple and practical way is proposed for detecting and avoiding these fictitious solutions. The concepts are illustrated with examples of a scattering wave on a rigid sphere.

Tables

Table 1. Table 1: Values of the fictitious frequency that correspond to scattering by a rigid sphere with Neumann boundary condition. The three lowest values that are the solutions to the eigenvalue equation at each m 𝑚 m value given in the right most column are given to 6-7 significant figures.

$m$		$k_{f}R$		Equation: $x\equiv k_{f}R$
	$1^{st}$	$2^{nd}$	$3^{rd}$
0	3.14159	6.28319	9.424778	$\tan x=0$
1	4.49341	7.72525	10.90412	$\tan x=x$
2	5.763459	9.095011	12.32294	$\tan x=\frac{3x}{3-x^2}$
3	6.987932	10.41712	13.69802	$\tan x=\frac{15x-x^3}{15-6x^2}$
4	8.182561	11.70491	15.03966	$\tan x=\frac{105x-10x^3}{105-45x^2+x^4}$
5	9.355812	12.96653	16.35471	$\tan x=\frac{945x-105x^3+x^5}{945-420x^2+15x^4}$

Equations

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

$$c({\mathbf{x}}_0)\varphi ({\mathbf{x}}_0)+{\int }_S\varphi (\mathbf{x})\frac{\partial G(\mathbf{x},{\mathbf{x}}_0|k)}{\partial n}\text{ d}S(\mathbf{x})={\int }_S\frac{\partial \varphi (\mathbf{x})}{\partial n}G(\mathbf{x},{\mathbf{x}}_0|k)\text{ d}S(\mathbf{x}),$$

$$G(\mathbf{x},{\mathbf{x}}_0|k)=\frac{e^{ikr}}{r}$$

$$c({\mathbf{x}}_0)f({\mathbf{x}}_0|k_f)+{\int }_{S}f(\mathbf{x}|k_f)\frac{\partial G(\mathbf{x},{\mathbf{x}}_0|k_f)}{\partial n}\text{ d}S=0.$$

$${\int }_{S}g(\mathbf{x}|k_f)G(\mathbf{x},{\mathbf{x}}_0|k_f)\text{ d}S=0,$$

$$c\Psi ({\mathbf{x}}_0)+{\int }_S\Psi (\mathbf{x})\frac{\partial G(\mathbf{x},{\mathbf{x}}_0|k)}{\partial n}\text{ d}S={\int }_S\frac{\partial \Psi (\mathbf{x})}{\partial n}G(\mathbf{x},{\mathbf{x}}_0|k)\text{ d}S.$$

$$\underset{\mathbf{x}\to {\mathbf{x}}_0}{\lim }\Psi (\mathbf{x})=\varphi ({\mathbf{x}}_0)$$

$$\underset{\mathbf{x}\to {\mathbf{x}}_0}{\lim }\frac{\partial \Psi (\mathbf{x})}{\partial n}=\frac{\partial \varphi ({\mathbf{x}}_0)}{\partial n}$$

$$\displaystyle\Psi(\boldsymbol{x})=\cos\big{(}k\boldsymbol{n}(\boldsymbol{x}_{0})$$

$$\displaystyle+\frac{1}{k}\sin\big{(}k\boldsymbol{n}(\boldsymbol{x}_{0})\cdot[\boldsymbol{x}-\boldsymbol{x}_{0}]\big{)}\frac{\partial\phi(\boldsymbol{x}_{0})}{\partial n}.$$

$$\displaystyle 4\pi\phi(\boldsymbol{x}_{0})+\int_{S}\big{[}\phi(\boldsymbol{x})-\Psi(\boldsymbol{x})\big{]}$$

$$\displaystyle\int_{S}\Big{[}\frac{\partial\phi(\boldsymbol{x})}{\partial n}-\frac{\partial\Psi(\boldsymbol{x})}{\partial n}\Big{]}G(\boldsymbol{x},\boldsymbol{x}_{0}|k)\text{ d}S.$$

$$\frac{\partial \Psi (\mathbf{x})}{\partial n}=\mathbf{n}\cdot \nabla \Psi =$$

$$\displaystyle+\boldsymbol{n}(\boldsymbol{x})\cdot\boldsymbol{n}(\boldsymbol{x}_{0})\cos\big{(}k\boldsymbol{n}(\boldsymbol{x}_{0})\cdot[\boldsymbol{x}-\boldsymbol{x}_{0}]\big{)}\frac{\partial\phi(\boldsymbol{x}_{0})}{\partial n},$$

$$\text{Mean Error}=\frac{\sqrt{{\sum }_{i=1}^{\text{DOF}}{\left(|{\varphi }_{\text{num}}^i|-|{\varphi }_{\text{ana}}^i|\right)}^2}}{\text{DOF}},$$

$$2\pi +{\int }_S\frac{\partial G(\mathbf{x},{\mathbf{x}}_0|k_f)}{\partial n}\text{ d}S(\mathbf{x})=0.$$

$$\int_{S}\frac{\partial G(\boldsymbol{x},\boldsymbol{x}_{0}|k_{f})}{\partial n}\text{ d}S(\boldsymbol{x})=-2\pi\Big{\{}e^{i2k_{f}R}+\frac{1}{ik_{f}R}\big{[}1-e^{i2k_{f}R}\big{]}\Big{\}}$$

$$\sin (k_{f}R)[1-i{k}_{f}R]=0.$$

$$\sin (k_{f}R)=0\text{or}{k}_{f}R=\pi ,2\pi ,3\pi ....$$

$$\tan (k_{f}R)=k_{f}R.$$

$$G_{\text{mod}}(\mathbf{x},{\mathbf{x}}_0|k)$$

$$=G(\mathbf{x},{\mathbf{x}}_0|k)+c_{2}G(\mathbf{x},\mathbf{a}|k)$$

$$\displaystyle c\phi(\boldsymbol{x}_{0})+\int_{S}\phi(\boldsymbol{x})\Big{[}$$

$$\displaystyle\qquad=\int_{S}\frac{\partial\phi(\boldsymbol{x})}{\partial n}\big{[}G(\boldsymbol{x},\boldsymbol{x}_{0}|k)+c_{2}G(\boldsymbol{x},\boldsymbol{a}|k)\big{]}\text{ d}S(\boldsymbol{x}).$$

$$2\pi+\int_{S}\Big{[}\frac{\partial G(\boldsymbol{x},\boldsymbol{x}_{0}|k_{f})}{\partial n}+\frac{\partial G(\boldsymbol{x},\boldsymbol{a}|k_{f})}{\partial n}\Big{]}\text{ d}S(\boldsymbol{x})=0.$$

$$\sin (k_{f}R)+c_2(R/a)\sin (k_{f}a)=0.$$

$$2\pi R+\int_{S}z\Big{[}\frac{\partial G(\boldsymbol{x},\boldsymbol{x}_{0}|k_{f})}{\partial n}+c_{3}\frac{\partial G(\boldsymbol{x},\boldsymbol{a}|k_{f})}{\partial n}\Big{]}\text{ d}S(\boldsymbol{x})=0.$$

$${\int }_{S}z\frac{\partial G(\mathbf{x},\mathbf{a}|k_f)}{\partial n}\text{ d}S(\mathbf{x})$$

$${\int }_{S}z\frac{\partial G(\mathbf{x},\mathbf{a}|k_f)}{\partial n}\text{ d}S(\mathbf{x})=$$

$$\displaystyle\frac{2\pi R}{a}\int_{R-a}^{R+a}\frac{R^{2}+a^{2}-r^{\prime 2}}{2a}\Big{[}-R+\frac{R^{2}+a^{2}-r^{\prime 2}}{2R}\Big{]}\frac{e^{ikr^{\prime}}}{r^{\prime 2}}[ikr^{\prime}-1]\text{ d}r^{\prime}$$

$${\int }_{S}x\frac{\partial G(\mathbf{x},\mathbf{a}|k_f)}{\partial n}\text{ d}S(\mathbf{x})={\int }_{S}y\frac{\partial G(\mathbf{x},\mathbf{a}|k_f)}{\partial n}\text{ d}S(\mathbf{x})=0$$

$$c\varphi (\mathbf{a})+{\int }_S\varphi (\mathbf{x})\frac{\partial G(\mathbf{x},\mathbf{a}|k)}{\partial n}\text{ d}S={\int }_S\frac{\partial \varphi (\mathbf{x})}{\partial n}G(\mathbf{x},\mathbf{a}|k)\text{ d}S,$$

$$4\pi \varphi ({\mathbf{x}}_0)=-{\int }_S\varphi (\mathbf{x})\frac{\partial G(\mathbf{x},{\mathbf{x}}_0|k)}{\partial n}\text{ d}S+{\int }_S\frac{\partial \varphi (\mathbf{x})}{\partial n}G(\mathbf{x},{\mathbf{x}}_0|k)\text{ d}S,$$

Full text

††journal: Engineering Analysis with Boundary Elements

1 Introduction

Recent studies of the boundary integral formulation of problems in time domain acoustic scattering [1], wave propagation in periodic structures [2], dynamic elasticity using the Helmholtz decomposition method [3] and direct field-only formulation of computational electromagnetics [4, 5, 6], all rely on finding accurate and efficient methods of solving the scalar Helmholtz equation. In this regard, it is timely to re-visit the boundary integral method of solving the Helmholtz equation.

It is well-known that the solution of the Helmholtz equation for external problems obtained by the boundary integral method, BIM, (or its numerical counterpart: boundary element method BEM) can become non-unique at certain frequencies (similar problems can occur for multiply connected domains [7]). At these so called fictitious [8] or irregular [9] frequencies, the non-physical solutions that arise are said to correspond to the internal resonance frequencies of the scatterer. Although there are established methods, most notably due to Schenck [10] and to Burton and Miller [11] that have been developed to eliminate such fictitious solutions, these methods require numerical tools beyond the BIM. For instance, the solution of Schenck requires additional numerical algorithms such as least squares minimization and that of Burton and Miller leads to hypersingular integral equations [12, 13, 14, 15]. Here we show that these fictitious solutions, when they do occur, and their corresponding frequencies in the BIM context depend not only on the shape of the object but also on the choice of Green’s function so that these frequencies do not necessarily occur at the corresponding internal resonance frequencies of the object. This observation together with the fact that recently developed desingularized BIM can give sufficiently high precision that the solution is unaffected by such fictitious solutions until the frequency is within about 1 part in $10^{4}$ of a fictitious value. We shall demonstrate how this can be exploited to detect the presence of a fictitious solution. Furthermore, the fictitious frequency spectrum can be changed by using different Green’s functions in the BIM. Taken together, these developments provide a practical way to detect and eliminate the effects of the fictitious solution without additional numerical effort or adjustable parameters beyond the toolkit of the BIM.

The introduction of a modified Green’s function also poses a number of interesting but unanswered questions that can provide stimulus for further theoretical development.

To provide physical context to our discussion on how the fictitious solution arises in the solution of the Helmholtz wave equation using the boundary integral method, we consider the example of the scattering of an incident acoustic wave by an object with boundary $S$ in an infinite medium. In the external domain, assumed to be homogeneous, scattered acoustic oscillations are described by the Helmholtz scalar wave equation in the frequency domain:

[TABLE]

where $k=\omega/c$ is the wave number, $\omega$ the angular frequency and $c$ the speed of sound. The (complex) acoustic potential, $\phi$, is related to the scattered velocity: $\boldsymbol{u}=\nabla\phi$. Since Eq. 1 is elliptic, the Green’s function formalism can be used to express the solution as that of a boundary integral equation [16, 17]

[TABLE]

where

[TABLE]

is the 3D Green’s function with $r=\|\boldsymbol{x}-\boldsymbol{x}_{0}\|$ and $\partial/\partial n\equiv\boldsymbol{n}\cdot\nabla$ is the normal derivative where the normal vector $\boldsymbol{n}$ points out of the domain, and thus into the object. The position vector $\boldsymbol{x}$ in Eq. 1 is located on the boundary $S$. If the observation point $\boldsymbol{x}_{0}$ is located outside the object (i.e. within the solution domain), the solid angle $c=4\pi$, if $\boldsymbol{x}_{0}$ is located inside the object (i.e. outside the solution domain), $c=0$, and if $\boldsymbol{x}_{0}$ is located on the surface, $S$, of the object and that point on $S$ has a continuous tangent plane, then and only then $c=2\pi$, otherwise the value of the solid angle $c$ is determined by the local surface geometry at $\boldsymbol{x}_{0}$.

The advantages of using Eq. 2 over other methods such as using finite difference in the 3D domain are the obvious reduction in the spatial dimension by one and that it is relatively easy to accommodate complicated shapes without deploying multi-scale 3D grids. Also the Sommerfeld radiation condition at infinity [18] is automatically satisfied by Eq. 2.

For the simple example of the scattering of an incoming plane wave specified by $\phi^{\text{inc}}=\Phi_{0}e^{i\boldsymbol{k}\cdot\boldsymbol{x}}$ (with $\Phi_{0}$ a constant and $\|\boldsymbol{k}\|=k$) by a rigid object, the velocity potential, $\phi$ of the scattered wave can be found by solving Eq. 1. The condition of zero normal velocity on the surface is equivalent to the boundary condition on $S$: $\partial\phi/\partial n=-\partial\phi^{\text{inc}}/\partial n$. In this case, the right hand side of Eq. 2 is known so this equation can be solved for the velocity potential, $\phi(\boldsymbol{x}_{0})$, with $\boldsymbol{x}_{0}$ on the surface.

We now demonstrate using this example of a Neumann problem where $\partial\phi/\partial n$ is given on the surface $S$, that there exists certain values of $k=k_{f}$, at which the solution $\phi$ of Eq. 2 is no longer unique. This occurs at those frequencies $k_{f}$ whereby a non-trivial function $f$ can exist to satisfy the following homogeneous equation:

[TABLE]

Consequently Eq. 2 will admit a solution of the form $\phi+bf$ on the surface $S$, where $b$ is an arbitrary constant and $f$, the fictitious solution, also satisfies the integral equation with zero normal derivative on $S$. Thus the fictitious frequency, $k_{f}$ and the corresponding fictitious solution, $f(\boldsymbol{x}|k_{f})$ are the eigenvalue and eigenfunction of Eq. 4, respectively. The existence of fictitious frequencies in boundary integral methods for Helmholtz equations was already identified by Helmholtz in 1860 [19], who said on page 24 (see also page 29 of his book [20]), while discussing the integral equation, Eq. 2:

…aber für eine unendlich grosse Zahl von bestimmten Werthen von $k$ für eine jede gegebene geschlossene Oberfläche Ausnahmen erleidet. Es sind dies nämlich diejenigen Werthe von $k$, die den eigenen Tönen der eingeschlossenen Luftmasse entsprichen.

This text was more or less translated directly by Rayleigh [21] in his book:

For a given space $S$ there is …. a series of determinate values of $k$, corresponding to the periods of the possible modes of simple harmonic vibration which may take place within a closed rigid envelope having the form of $S$. With any of these values of $k$, it is obvious that $\phi$ cannot be determined by its normal variation over $S$, and the fact that it satisfies throughout $S$ the equation $\nabla^{2}\phi+k^{2}\phi=0$.

Note that the internal resonance problem corresponds to a problem with $\phi=0$ on the surface, $S$ and $g\equiv\partial\phi/\partial n\neq 0$ in Eq. 2, is given by

[TABLE]

which is different from Eq. 4. It is not immediately obvious that Eqs. 4 and 5 will produce the same fictitious spectrum and in fact, as we shall see later in Section 4, this is not always the case.

In theory, the fictitious solution only appears if $k$ is exactly equal to $k_{f}$ so that it is not an issue in analytic work nor if computations have infinite numerical precision. With the advent of numerical techniques in the late 1960’s and early 1970’s, the boundary integral equation was transformed into the boundary element method (BEM). The issue of fictitious frequencies now resurfaced once more in the numerical implementations. In the conventional implementation of the BEM [17], the surface $S$ is represented by a mesh of planar area elements and the unknown value of $\phi(\boldsymbol{x})$ on the surface is assumed to be a constant within each planar element and only varies from element to element. The surface integral is thus converted to a linear system in which the values of $\phi$ at different area elements are unknowns to be solved. The practicality of discretization where the representation of the surface $S$ by a finite number of planer elements and round off errors in numerical computation mean that effects of the fictitious solution begin to be important, not only when $k=k_{f}$, but even when the value of $k$ is near $k_{f}$. For instance, in a conventional implementation of the BEM, the apparent location of the fictitious frequency, $k_{f}$ can be in error because of the approximation involved in representing the actual surface by a set of planar elements. Thus the mean relative error can exceed 100% when $k$ is within 1-2% of the actual fictitious frequency (see Fig.1 for examples of a sphere with radius $R$ at $kR\approx\pi$ and $kR\approx 2\pi$). Since the values of $k_{f}$ are not known a priori for general boundary shapes, $S$, the accuracy of any BEM solution of the Helmholtz equation can become problematic.

Two popular methods to deal with this issue that are still in use today are due to Schenck [10] and to Burton and Miller [11]. Schenck introduced the CHIEF method whereby the BEM solution is evaluated at additional internal points inside the scatterer with the requirement that such values must vanish. This results in an over-determined matrix system that requires a least square solution entailing considerable additional computational time, especially for larger systems. However, the CHIEF method does not stipulate how many CHIEF points should be used and where they should be placed. The Burton and Miller [11] method involves taking the normal derivative of Eq. 2, multiplying it by an appropriate complex number and then adding it to the original equation. It is claimed that Eq. 2 and its normal derivative have different resonance spectra and this therefore solves the fictitious frequency problem. Due to the use of the normal derivative of Eq. 2, the Burton and Miller method involves having to deal with strongly singular kernels. This approach therefore has the disadvantage that it requires special quadrature rules for higher order elements [22].

The issue of fictitious solutions is revisited in this article. Clearly, if a numerical implementation of the BEM is not sensitive to the fact that $k$ may be close to a fictitious value $k_{f}$, then the effects of a fictitious solution will be minimized. Furthermore, the spectrum of fictitious frequencies does not only depend on the shape of the object, but also on the choice of the Green’s function. As the classical free space Green’s function or fundamental solution of Eq. 3 is not the only choice that can be used, it can be replaced by other fundamental solutions, as long as they are analytic in the external domain and they satisfy the Sommerfeld radiation condition [23]. Thus using a different Green’s function will shift the spectrum of fictitious frequencies relative to a given $k$ value. Although the theoretical framework of modified Green’s functions has been discussed extensively in the literature [23, 24, 25, 26, 27], only very little attention appears to have been paid to the actual implementation, for example, the cases of Neumann boundary condition and of Dirichlet boundary condition were considered in two nearly identical papers [28, 29]. In this article we address this issue.

The development of our suggestion to eliminate the fictitious frequency problem in BEM solutions of the external Helmholtz equation is organized as follows. In Section 2, we outline how a desingularized implementation of the BEM that is not affected by a fictitious solution unless $k$ is very close to a fictitious value $k_{f}$, can be used to decide if an BEM solution has been adversely affected by the presence of a fictitious component. This framework also enables us to implement higher order elements with ease. In Section 3 the spectrum of fictitious frequencies and corresponding solutions are studied as the solution of an homogeneous integral equation. In Section 4, a modified Green’s function is introduced to show how it can be used to change the spectrum of fictitious frequencies. Thus by employing the desingularized BEM, it is sometimes easy to determine by comparing the solutions obtained from using the conventional Green’s function in Eq. 3, and from a modified Green’s function whether the solutions have been adversely affected by the presence of a solution associated with a fictitious frequency. Some discussion and the conclusion follow in Sections 5 and 6, respectively.

2 Minimize the proximity effects to a fictitious frequency

As noted earlier, discretization and round off errors can cause the spurious solution to become important when the wave number happens to be near a fictitious value. However, since the spectrum of fictitious frequencies is not generally known a priori, the numerical accuracy of a solution obtained by the BEM becomes uncertain. Therefore, to ameliorate the fictitious frequency problem, it is valuable to have an accurate implementation of the BEM that will not produce a fictitious component to the solution unless the frequency $k$ is extremely close to an unknown fictitious frequency. This is provided by a recently developed fully desingularized boundary element formulation [30, 31], a concept that was first introduced for the BEM solution of the Laplace equation by Klaseboer et al. [32]. In this framework, the traditional singularities of the Green’s function and its normal derivative in the BEM integrals are removed analytically from the start.

High accuracy can be achieved in this approach firstly due to the fact that all elements (including the previously singular one) are treated in the same manner with the same Gaussian quadrature scheme. The second reason for the high accuracy lies in the fact that instead of using planar area elements in which the unknown functions are assumed to be constant within such elements, the unknowns are now function values at node points on the surface, and the surface is represented more accurately by quadratic area elements determined by these nodal points. In calculating integrals over the surface elements, variation of the function value within each element is also estimated by quadratic interpolation from the nodal values. The numerical implementation is straightforward, once the linear system is set up, the usual linear solvers can be used. The thus obtained framework is termed Boundary Regularized Integral Equation Formulation (or BRIEF in short [31]).

Here is a brief description of the desingularized boundary element formulation, details of which are given in previous works [30, 31]. Assume we have a known analytic solution, $\Psi(\boldsymbol{x})$, of the Helmholtz equation in Eq. 1 which then also satisfies Eq. 2 as:

[TABLE]

Without loss of generality, we can demand that this solution further satisfies the following two point-wise conditions when $\boldsymbol{x}_{0}$ is on surface $S$:

[TABLE]

[TABLE]

A convenient but not the only possible choice is a combination of two standing waves, one with the node of the wave and the other with the antinode situated at $\boldsymbol{x}_{0}$, both aligned with $\boldsymbol{n}(\boldsymbol{x}_{0})$ [30] as:

[TABLE]

Substituting Eq. 9 in Eq. 6 and subtracting the result from Eq. 2 gives:

[TABLE]

The conditions from Eqs. 7 and 8 guarantee that the terms in $[...]$ on both sides of Eq. 10 cancel out the singularities of the Green’s function and its derivative by noting that

[TABLE]

and the fact that $\boldsymbol{n}(\boldsymbol{x})\cdot\boldsymbol{n}(\boldsymbol{x}_{0})\to 1$, when $\boldsymbol{x}$ approaches $\boldsymbol{x}_{0}$ for any smooth surface. Note that the solid angle in Eq. 10 has been eliminated, but a term with $4\pi\phi(\boldsymbol{x}_{0})$ appears due to the contribution of the integral over a surface at infinity because of the particular choice of Eq. 9. Also note from Eq. 9 that $\Psi$ is a different function for each node on the surface. It is noted that other desingularization methods based on entirely different concepts exist as well in the literature [33].

We now consider the example of solving the scattering problem by a solid sphere with radius $R$ for which the spectrum of resonant frequencies of the internal problem is known. A list of the values of the lower resonant frequencies and the equation that generates them are given in Table 1 where we see that two of the lowest frequencies are at $k_{f}R=\pi$ and $k_{f}R=2\pi$. With the choice of Eq. 3 for the Green’s function, the spectrum fictitious frequencies coincides with the resonant frequencies of the corresponding internal problem. In Fig. 1, we quantify the behaviour of the BEM solution for $kR$ values in the neighborhood these 2 fictitious values in terms of the mean square error defined by

[TABLE]

where $\phi_{\text{num}}^{i}$ and $\phi_{\text{ana}}^{i}$ are, respectively, the numerical (BEM) and analytic solution at node $i$. The number of nodes used in the desingularized BEM, the Degree of Freedom (DOF), is around 2000. We see that the mean squared error even in the small neighborhoods $0.94\pi<kR<1.06\pi$ and $1.94\pi<kR<2.06\pi$ around the 2 fictitious frequencies is extremely localized. In fact, the BEM solutions obtained by the desingularized BEM [30, 31] are unaffected by fictitious solutions until the frequency is within about 1 part in $10^{4}$ of a fictitious value. The results for the conventional boundary integral method (CBIM) are also shown. Note that the fictitious frequency predicted by the CBIM is significantly higher than the known theoretical value in these examples presumably because of the effect of approximating the sphere surface by a finite set of planar elements.

Similar remarks apply for the behavior of the desingularized BEM solution in the neighborhood of the lowest $m=1$ fictitious value $k_{f}R=4.49341$ (see Table 1) shown Fig. 2. Here we show the values of the real and imaginary parts of the solution of nodes at the front and at the back of the sphere. The effect of the fictitious solution can only be discerned in the very narrow window $4.493<kR<4.494$ around $k_{f}R=4.49341$. But outside this window, there is no noticeable effect due to $kR$ being close to the fictitious value, $k_{f}R$. For example, if at the values $kR=3.140$ and $kR=4.490$ as given in Kinsler [16], page 518, the desingularized BEM (BRIEF) was used to to solve the Helmholtz equation, the solution would not register as giving fictitious results. As we shall see below, if a sweep of 10,000 frequencies from $kR=0$ to 10 is performed in steps of 0.001 one would miss many fictitious solutions (since a step size of 0.001 would not be precise enough to detect all of them).

From the above results, we can conclude that the effects of resonance are not observed until one is extremely close to the resonant frequency in our desingularized BEM [30, 31].

3 The genesis of fictitious solutions

In the example of acoustic scattering by a rigid scatterer that was discussed in the previous section, $\partial\phi/\partial n$ on the surface of the scatterer is specified (Neumann boundary conditions), and the variation of $\phi$ on the surface is the unknown to be found. At certain frequencies however, instead of the expected $\phi$, another function say, $\phi+f$ emerges. The frequencies at which this occurs, are often said to correspond to the internal resonance frequency of the same object. The Fredholm integral theory has been used to explain the occurrence of the fictitious frequency and related the fictitious frequency to the corresponding internal resonance frequency of the object [34]. However, by working directly with the integral equation that determines the fictitious solution, $f$, it is easy to demonstrate the origin of the fictitious solution.

First we use the example of scattering on a rigid sphere of radius, $R$, to demonstrate how the fictitious solution and frequency is determined by the Green’s function and the boundary shape. For simplicity, we consider the solution of the Helmholtz equation outside a sphere that has azimuthal symmetry for which the solution on the sphere surface can be expanded in terms of Legendre polynomials of order $m$, $P_{m}(\cos\theta)$ to account for variations in the polar angle, $\theta$. In this case, the fictitious frequencies for different $m$ values are known. We consider in detail the fictitious solution, $f$, and the fictitious frequency, $k_{f}$ for the cases with $m=0$ and $m=1$.

3.1 Case: $f\sim P_{0}(\cos\theta),\text{a constant}$, $m=0$

In this case, the fictitious solution, $f$ is a constant, being proportional to $P_{0}(\cos\theta)$, on the surface of the sphere of radius, $R$ and $c(\boldsymbol{x}_{0})=2\pi$, then Eq. 4, at the fictitious wave number, $k_{f}$, becomes:

[TABLE]

The integral of $\partial G(\boldsymbol{x},\boldsymbol{x}_{0}|k)/\partial n$, can be evaluated (see A) to give

[TABLE]

so that Eq. 13 is equivalent to

[TABLE]

Thus the spectrum of fictitious frequencies corresponding to a constant fictitious function, $f\sim P_{0}(\cos\theta)$, on the surface with $m=0$ is

[TABLE]

see also the first row of Table 1. In the external 3D domain, the fictitious solution $f(\boldsymbol{x})$ that emerges numerically from the BEM solution corresponding to $k_{f}R=\pi$ is: $f(\boldsymbol{x})=c_{3}\;e^{ik_{f}\|\boldsymbol{x}\|}/\|\boldsymbol{x}\|$, where $c_{3}$ is an arbitrary constant and the origin of $\boldsymbol{x}$ taken at the origin of the sphere.

3.2 Case: $f\sim P_{1}(\cos\theta)$, $m=1$

A similar calculation to the one given in Section 3.1, for a fictitious function, $f\sim P_{1}(\cos\theta)$, for $m=1$ leads to (see B)

[TABLE]

The first few solutions to Eq. 17 are given in the $m=1$ row of Table 1. Again, these values are equal to those of the corresponding internal eigenvalue problem, yet they have been derived here purely from a boundary integral equation perspective. Fictitious frequencies for higher order values of $m$ can also be obtained in a similar manner. Table 1 contains all fictitious frequencies below $k_{f}R=10$ for a sphere.

The above derivation that starts from the homogeneous integral equation, Eq. 4 demonstrates the role of the Green’s function and the boundary shape in determining the spectrum of fictitious frequencies and solutions for acoustic scattering by a solid sphere. We can now show how to modify the fictitious frequency spectrum using different Green’s functions.

4 The modified Green’s function

Different forms of the Green’s function can be used to construct the integral equation of the BEM as long as they satisfy the same differential equation in the solution domain and the Sommerfeld radiation condition at infinity as the free space Green’s function. A simple modified Green’s function, $G_{\text{mod}}$, can be taken as

[TABLE]

where the origin is taken to be the center of the sphere and the vector $\boldsymbol{a}$ corresponds to a point inside the sphere $(|\boldsymbol{a}|<R)$ with $c_{2}$ an arbitrary constant. The integral equation that implements the BEM with $G_{\text{mod}}$ becomes:

[TABLE]

The additional term $G(\boldsymbol{x},\boldsymbol{a}|k)$ although singular at the location $\boldsymbol{a}$, does not create any singular behavior on the surface $S$, since $\|\boldsymbol{x}-\boldsymbol{a}\|$ never becomes zero (see also Fig. 3). The modified Green’s function, $G_{\text{mod}}(\boldsymbol{x},\boldsymbol{x}_{0}|k)$, also satisfies the Sommerfeld radiation condition at infinity.

4.1 Case: $f\sim P_{0}(\cos\theta),\text{a constant}$, $m=0$ with modified $G_{\text{mod}}$

Let us now investigate how the modified Green’s function defined in Eq. 18 and 19 can affect the spectrum of fictitious frequencies that is now determined by

[TABLE]

Evaluating the integrals (see C) then gives the equation that determines the spectrum of fictitious frequencies

[TABLE]

Thus the original fictitious frequency spectrum given by $\sin(k_{f}R)=0$ in Eq. 16 due to the use of the unmodified Green’s function in Fig. 4a has been replaced by a different spectrum given by Eq. 21 in Fig. 4b. Furthermore, the precise value of $a$ is not critical. In fact, since $\sin(k_{f}a)/a\rightarrow k_{f}$ as $a\rightarrow 0$, we can put $\boldsymbol{a}=\boldsymbol{0}$, that is, at the center of the sphere. In the results shown in Fig. 4, we have taken $a=0$ and $c_{2}=-1$.

For $m=0$, Eq. 21 will in general assure that the new fictitious frequency spectrum obtained with the modified Green’s function, $G_{\text{mod}}$ will be different from that obtained with the original Green’s function, $G$. However, there are still ways for which this may not be true.

Firstly, it is still possible that both $\sin(k_{f}R)$ and $\sin(k_{f}a)$ vanish, that is, the original spectrum and the modified spectrum contain common values. An example of such a case can be observed when $k_{f}R=2\pi$ and $a=0.5R$ (thus $k_{f}a=\pi$ and $\sin(k_{f}a)=0$). This was tested numerically and indeed for these parameters there is still a spurious solution corresponding to the common fictitious frequency values in the 2 spectra as illustrated in Fig. 5.

• 2.

A second way in which fictitious behaviour can still be observed, is when for particular parameters of $k_{f}$, $R$, $a$ and $c_{2}$, Eq. 21 is still zero. An instance of such fictitious behavior can be observed for the parameters $k_{f}R=0.5$, $a=0.3R$ and $c_{2}=-0.9624563$. The fictitious solution for these parameters is about 100 times the theoretical value in a numerical test. It is interesting to note that a fictitious frequency now appears at $k_{f}R=0.5$, a frequency value that was previously free of fictitious behavior. This is an example of a frequency shift of the lowest fictitious behavior from $k_{f}R=\pi$ to a lower frequency of $k_{f}R=0.5$. However, if $c_{2}=-0.9620000$ is chosen, thus only slightly different from $c_{2}=-0.9624563$, no fictitious behavior is observed at all (see Fig. 6).

• 3.

Finally, the location of the point $\boldsymbol{a}$ should not be chosen too close to the boundary $S$. In order to investigate this, in Fig. 7, the potentials in front and at the back of the sphere are shown, while the location of $\boldsymbol{a}$ of the modified Green’s function is varied from $a=0.0$ to $1.0$. From the figure it can be deduced that $\boldsymbol{a}$ should not be placed closer to the boundary $S$ than roughly the meshsize.

To conclude, for $m=0$, the modified Green’s function approach can indeed remove the fictitious behavior of the solution. In the next section, the $m=1$ case will be investigated.

4.2 The $m=1$ case

In Section 4.1, it was shown that for $f=$constant (or $m=0$), the modified Green’s function can indeed remove the fictitious solutions. A similar proof can now be attempted for $m=1$. In analogy to Eq. 33, it must now be shown that

[TABLE]

Thus the integral

[TABLE]

must be determined. The framework of Eqs. 44, 46 can be adapted immediately, provided that we add $z$ in the equations. With $z=R\cos\alpha=[R^{2}+a^{2}-r^{\prime 2}]/(2a)$:

[TABLE]

This integral can be shown not to be equal to zero. However, a similar calculation for $x$ or $y$ instead of $z$, shows that due to symmetry (provided that $\boldsymbol{x}_{0}$ is still situated on the z-axis):

[TABLE]

From this we can conclude that, unfortunately, the fictitious solutions corresponding to $m=1$ cannot be removed when applying our modified Green’s function in its present form. This is also clear from Fig. 4b, the fictitious behavior corresponding to $m=1$ is still present. A more elaborate Green’s function might still be capable of removing these frequencies as well, but this is beyond the scope of the current article, in which we intend merely to show the proof of concept.

5 Discussion

In both the modified Green’s function and in the CHIEF method, a point in the interior of the domain is chosen on which an integral equation for $G(\boldsymbol{x},\boldsymbol{a})$ is developed. The difference between the modified Green’s function and CHIEF, however, lies in the fact that CHIEF uses the following equation as an extra condition to the system of equations:

[TABLE]

Here, the constant $c=0$, since the point $\boldsymbol{a}$ is situated outside the domain (i.e. inside the object) in the CHIEF method. In the modified Green’s function approach this equation is essentially added to the ’normal’ Green’s function.

A way to check if the solution using our desingularized boundary element code for a general shaped object contains a fictitious component due to $k$ being close to a fictitious value is to repeat the calculation at a very slightly different $k$ value. If the solution differs significantly, the solution is likely to contain a fictitious component.

We further illustrate the concepts with some field values of $\phi$ obtained by post-processing from the following equation

[TABLE]

where $x_{0}$ is not situated on the boundary $S$, but either in the solution domain or inside the sphere (outside the solution domain). If no resonance is present, the solution inside the sphere (and hence outside the solution domain) should be $\phi=0$. In for following examples we use 1442 nodes and 720 quadratic elements in the BEM solution. The first case is the solution for $kR=2\pi$ where in Fig. 8 we plotted the results obtained from both the standard BEM (with fictitious results, Fig. 8a and that obtained using a modified Green’s function Fig. 8b, where the solution inside the sphere is zero.

A second example shows the solution for the resonance frequency $kR=4.49341$ in Fig. 9a. At a frequency nearby at $kR=4.49000$ no resonance behavior is observed in Fig. 9b. This once more demonstrates the extreme accuracy of our desingularized BEM framework.

A third example shows the resonance behavior at $kR=5.76345$ and a nearby value of $kR=5.76000$ in Fig. 10. Again no resonant behavior is observed at the nearby value.

A final example shows the solution at $kR=3\pi$ in Fig. 11, obtained from both with the standard method and with the modified Green’s function.

At present, the modified Green’s function can only remove solutions associated with fictitious frequencies in the “breathing modes” ($m=0$), but these are most likely the first modes to appear with increasing $k$. It would be interesting to find other modified Green’s functions to remove solutions associated with fictitious frequencies in all modes, but we have not as yet been able to develop such an approach.

6 Conclusions

The fictitious frequencies occurring in a BEM implementation of the Helmholtz equation were revisited. From a BEM viewpoint it was highlighted how these fictitious solutions appear and how they can be detected. It was shown that the use of a modified Green’s function can indeed remove certain fictitious frequencies. To the best knowledge of the authors, this is the first time actual numerical results have been obtained with a modified Green’s function. The results presented are a demonstration of the proof of concept. More elaborate modified Green’s functions might be able to remove more fictitious frequencies. If indeed so, then this easy to implement method could be a viable alternative to existing methods.

Fictitious frequencies cannot fully be avoided with the current alternative Green’s function approach, but it is sometimes possible to shift this frequency to another region of the spectrum. Thus the fictitious frequencies do not necessarily coincide anymore with a corresponding internal resonance frequency of the scatterer. The superior accuracy of the desingularized boundary element method further ensures that the fictitious behavior is limited to very narrow bands in the frequency spectrum. The concepts were illustrated with examples of the scattering of a plane wave on a rigid sphere.

Acknowledgements

This work is supported in part by the Australian Research Council grants, DE150100169 & CE140100003, to QS and a Discovery Project Grant to DYCC.

Appendix A Fictitious frequencies for the m=0 case

The normal derivative of the Green’s function, $\partial G(\boldsymbol{x},\boldsymbol{x}_{0}|k)/\partial n$, can be expressed as

[TABLE]

Without loss of generality lets assume that the point $\boldsymbol{x}_{0}$ is located on the z-axis (see also Fig. 12 for the definition of symbols), thus $\boldsymbol{x}_{0}=[0,0,R]$, the vectors $\boldsymbol{x}$ and $\boldsymbol{n}$ can then be presented by $\boldsymbol{x}=R[\cos\theta\sin\alpha,\sin\theta\sin\alpha,\cos\alpha]$ and $\boldsymbol{n}=-\boldsymbol{x}/R$. Then $(\boldsymbol{x}-\boldsymbol{x}_{0})\cdot\boldsymbol{n}=R(-1+\cos\alpha)$. The surface element $\text{d}S=2\pi R\rho\text{ d}\alpha$ can also be expressed as $\text{d}S=2\pi R^{2}\sin\alpha\text{ d}\alpha$:

[TABLE]

With the help of Fig. 12, the term $(-1+\cos\alpha)$ can be rewritten as: $(-1+\cos\alpha)=-2\sin^{2}(\alpha/2)=-r^{2}/(2R^{2})$. From $r=2R\sin(\alpha/2)$, one can deduce $R\text{ d}\alpha=\text{d}r/\cos(\alpha/2)$. With $\sin\alpha=\cos(\alpha/2)r/R$, the singular term $1/r^{3}$ in Eq. 29 will be eliminated and this integral will turn into

[TABLE]

which will finally transform Eq. 13 in:

[TABLE]

Multiplying this equation by $e^{-ikr}$ and rearranging leads to

[TABLE]

Appendix B Fictitious frequencies for the m=1 case

In Section 3.1 and A, it was shown how the fictitious frequencies appear for the simplest case of $f=constant$, corresponding to the lowest order Legendre polynomials with $m=0$. The next least complicated function will be a linear function, corresponding to $m=1$. For simplicity sake, lets take $f=z$ as an example. Taking again $\boldsymbol{x}_{0}$ on the z-axis will give $f(\boldsymbol{x}_{0})=R$ and with $c=2\pi$, Eq. 4 will turn into:

[TABLE]

Eq. 29 is still valid, except that an extra term $z=R\cos\alpha=R[1-r^{2}/(2R^{2})]$ must be included, thus Eq. 30 must be replaced by:

[TABLE]

The first integral in the last expression is the same that appeared in Section A as Eqs. 30, 31 (except for a factor $1/R$), the second integral can be evaluated as:

[TABLE]

Thus Eq. 34 becomes:

[TABLE]

where Eq. 33 was used in the last equality. Multiplying by $e^{-ikR}/(2\pi R)$ and regrouping terms with $e^{-ikR}$ and $e^{ikR}$ leads to:

[TABLE]

Expanding $e^{-ikR}$ and $e^{ikR}$ into $\cos(kR)$ and $\sin(kR)$ terms gives:

[TABLE]

Separating this into real and imaginary parts:

[TABLE]

Both the real and imaginary part lead to the following condition:

[TABLE]

which is the same as the internal resonance condition for $m=1$, with solution $kR=4.49341$ etc. (see Table 1).

Appendix C The m=0 case with a modified Green’s function

The normal derivative of the additional part is:

[TABLE]

As in Section 3 assume that $f=const$ (corresponding to $m=0$) and again assume that the point $\boldsymbol{x}_{0}$ is located on the z-axis, the vectors $\boldsymbol{x}$ and $\boldsymbol{n}$ and $\text{d}S$ are defined the same as in Section 3, while

$\boldsymbol{x}-\boldsymbol{a}=[R\cos\theta\sin\alpha,R\sin\theta\sin\alpha,R\cos\alpha-a]$. Thus $\boldsymbol{n}\cdot[\boldsymbol{x}-\boldsymbol{a}]=-R+a\cos\alpha$. For the length $r^{\prime}$ the following relationship can be found:

[TABLE]

while for $\text{d}r^{\prime}$ one finds:

[TABLE]

Thus, similar to Eq. 29:

[TABLE]

Substituting $\boldsymbol{n}\cdot[\boldsymbol{x}-\boldsymbol{a}]=-R+a\cos\alpha$ and eliminating $\cos\alpha$ with Eq. 42:

[TABLE]

The last equality can be obtained easiest by splitting the integral in two parts and using $\partial e^{ikr^{\prime}}/\partial r^{\prime}=e^{ikr^{\prime}}[ikr^{\prime}-1]/r^{\prime 2}$. Eq. 45 can be simplified to:

[TABLE]

Eq. 31 will now have an additional part as:

[TABLE]

Multiplying by $e^{-ikR}/(4\pi i)$ gives:

[TABLE]

Since the common term $[1+1/ikR]$ can never become zero (k is a real number), this finally simplifies to:

[TABLE]

If this equation is satisfied for a certain wave number, $k$, then $k_{f}=k$.