Boundary regularized integral equation formulation of the Helmholtz equation in acoustics

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

2. Introduction

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

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

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

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

3. Boundary regularized integral equation formulation

To illustrate the boundary regularized integral equation formulation (BRIEF) of the Helmholtz equation, we draw on the problem of calculating the acoustic pressure wave generated by a vibrating boundary specified by a closed surface, or more generally a set of closed surfaces denoted by S. The acoustic pressure, p(x) outside S, obeys the Helmholtz equation: ∇²p+k²p=0, where k is the wavenumber. The solution can be found by solving the conventional boundary integral equation that follows from using Green's second identity [10]

$$c_0\hspace{0.17em}p({\mathbf{\text{x}}}_0)+{\int }_{S+S_{\mathrm{\infty }}}p(\mathbf{\text{x}})\frac{\mathrm{\partial }G({\mathbf{\text{x}}}_0,\mathbf{\text{x}})}{\mathrm{\partial }n}\hspace{0.17em}\text{d}S(\mathbf{\text{x}})={\int }_{S+S_{\mathrm{\infty }}}\frac{\mathrm{\partial }p(\mathbf{\text{x}})}{\mathrm{\partial }n}G({\mathbf{\text{x}}}_0,\mathbf{\text{x}})\hspace{0.17em}\text{d}S(\mathbf{\text{x}}),$$3.1

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

Despite the fact that the physical problem may be well behaved on the boundaries, the mathematical singularities in (3.1) at x=x₀ owing to G and ∂G/∂nrequire careful treatment. To quote Jaswon [12] and Symm [13]: ‘Most integral equations of physical significance involve singular, or weakly singular, kernels, thereby hampering the procedures of both theoretical and numerical analysis’—a statement that is still valid up to this day.

Here we seek to eliminate these singularities analytically by exploiting the linear nature of (3.1) as follows.

But before we tackle (3.1) that corresponds to the Helmholtz equation for p(x), consider first the simpler case of the Laplace equation: ∇²q(x)=0 for which the corresponding boundary integral equation is

$$c_0\hspace{0.17em}q({\mathbf{\text{x}}}_0)+{\int }_{S+S_{\mathrm{\infty }}}q(\mathbf{\text{x}})\frac{\mathrm{\partial }{G}_0({\mathbf{\text{x}}}_0,\mathbf{\text{x}})}{\mathrm{\partial }n}\hspace{0.17em}\text{d}S(\mathbf{\text{x}})={\int }_{S+S_{\mathrm{\infty }}}\frac{\mathrm{\partial }q(\mathbf{\text{x}})}{\mathrm{\partial }n}{G}_0({\mathbf{\text{x}}}_0,\mathbf{\text{x}})\hspace{0.17em}\text{d}S(\mathbf{\text{x}}),$$3.2

where G₀(x₀,x)=1/|x−x₀|. The standard way to remove the singularity arising from ∂G₀/∂n on the left-hand side is to note that the constant q(x₀) also satisfies the Laplace equation with a vanishing normal derivative on the boundaries, so that equation (3.2) for [q(x)−q(x₀)] becomes

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

This is known as the ‘constant value subtraction’. However, the integral on the right-hand side of (3.3) still involves an integration over the singularity from G₀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{\text{x}}\to {\mathbf{\text{x}}}_0$, [q(x)−q(x₀)] vanishes as |x−x₀| or faster.

In this paper, we extend the approach [8] developed for the Laplace equation to remove all singularities in the boundary integral equation (3.1) for the Helmholtz problem. However, as a constant is not a solution of the Helmholtz equation, we need to construct an auxiliary function ψ(x) for a given value of x₀ in (3.1) that satisfies the Helmholtz equation that is analogous to the constant, q(x₀), in the case of the Laplace equation. We choose ψ(x) to be linear in the pressure p(x₀) and its normal derivative (∂p/∂n)₀ at x₀

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

Now provided the general functions g(x) and f(x) satisfy the Helmholtz equation with the following boundary conditions:

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

and

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

with n₀≡n(x₀) being the outward unit normal at x₀, ψ(x) will also satisfy the same boundary integral equation as (3.1). By taking the difference between the boundary integral equations for p(x) and ψ(x), we obtain

$$\begin{array}{rl}& {\int }_{S+S_{\mathrm{\infty }}}[p(\mathbf{\text{x}})-p({\mathbf{\text{x}}}_0)g(\mathbf{\text{x}})-{\left(\frac{\mathrm{\partial }p}{\mathrm{\partial }n}\right)}_{0}f(\mathbf{\text{x}})]\frac{\mathrm{\partial }G({\mathbf{\text{x}}}_0,\mathbf{\text{x}})}{\mathrm{\partial }n}\hspace{0.17em}\text{d}S(\mathbf{\text{x}})\\ & ={\int }_{S+S_{\mathrm{\infty }}}[\frac{\mathrm{\partial }p(\mathbf{\text{x}})}{\mathrm{\partial }n}-p({\mathbf{\text{x}}}_0)\mathrm{\nabla }g(\mathbf{\text{x}})\cdot \mathbf{\text{n}}(\mathbf{\text{x}})-{\left(\frac{\mathrm{\partial }p}{\mathrm{\partial }n}\right)}_0\mathrm{\nabla }f(\mathbf{\text{x}})\cdot \mathbf{\text{n}}(\mathbf{\text{x}})]G({\mathbf{\text{x}}}_0,\mathbf{\text{x}})\hspace{0.17em}\text{d}S(\mathbf{\text{x}}).\end{array}$$3.6

The conditions imposed on g(x) and f(x) at x₀ in (3.5) will remove both the singularities owing to G and ∂G/∂n in the new boundary integral equation in (3.6) under the generally valid assumption that as $\mathbf{\text{x}}\to {\mathbf{\text{x}}}_0$, [p(x)−ψ(x)] and [∂p/∂n−∂ψ/∂n] vanish as |x−x₀| or faster. Therefore, if g(x) and f(x) satisfy (3.5), all singular behaviour in (3.6) owing to G and ∂G/∂n will be removed. The formal proof is a straightforward generalization of that given in [6] for the Laplace equation.

Equation (3.6) is the final key result for the singularity-free [6] BRIEF of the solution of the Helmholtz equation that replaces the conventional boundary integral equation in (3.1). The integrals are now completely free of singularities [6] and can be evaluated by any convenient integration quadrature. Given values for p (Dirichlet) or ∂p/∂n (Neumann) or a relation between these two quantities on the boundary, S, (3.6) can be readily solved. Also the term involving the solid angle, c₀, in the conventional boundary integral equation, (3.1), has now been eliminated. Thus, higher order area elements, linear, quadratic or splines, can be used to represent the surface to improve numerical accuracy without the need to calculate the solid angle at each node.

For our key result in (3.6), we see that apart from having to satisfy (3.5), there is considerable freedom in selecting the functional forms of g(x) and f(x) and hence ψ(x). As a specific example, we can choose the following for g(x) and f(x):

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

where x_d is any convenient point outside the solution domain in order to ensure ψ(x) is non-singular within the solution domain, with a≡|x₀−x_d|, r_d≡|x−x_d| and b≡(x₀−x_d)⋅n(x₀)/|x₀−x_d|≠0 (figure 1a). For this particular choice of g(x) and f(x), the integral over the surface at infinity, $S_{\mathrm{\infty }}$ in (3.6) can be found analytically using the Sommerfeld radiation condition [14,15] (see appendix A for details)

$$p({\mathbf{\text{x}}}_0){\int }_{S_{\mathrm{\infty }}}[\mathrm{\nabla }g(\mathbf{\text{x}})\cdot \mathbf{\text{n}}G-g(\mathbf{\text{x}})\frac{\mathrm{\partial }G}{\mathrm{\partial }n}]\text{d}S=2\pi p({\mathbf{\text{x}}}_0)(1+\frac{\mathrm{i}}{ka})[1-{\mathrm{e}}^{2\mathrm{i}ka}]$$3.8

and

$${\left(\frac{\mathrm{\partial }p}{\mathrm{\partial }n}\right)}_0{\int }_{S_{\mathrm{\infty }}}[\mathrm{\nabla }f(\mathbf{\text{x}})\cdot \mathbf{\text{n}}G-f(\mathbf{\text{x}})\frac{\mathrm{\partial }G}{\mathrm{\partial }n}]\text{d}S={\left(\frac{\mathrm{\partial }p}{\mathrm{\partial }n}\right)}_0\frac{2\pi \mathrm{i}}{kb}[1-{\mathrm{e}}^{2\mathrm{i}ka}].$$3.9

In the limit k=0, the Helmholtz equation reduces to the Laplace equation with $g(\mathbf{\text{x}})\to 1$ and $f(\mathbf{\text{x}})\to (a/b)[1-a/|\mathbf{\text{x}}-{\mathbf{\text{x}}}_d|]$, and (3.6) will become a boundary regularized integral equation for the solution of the Laplace equation given earlier [8]. This is an example of how the singularity in the integral on the right-hand side of (3.3) may be removed analytically.

We note that the conventional boundary integral formulation of the solution of the Helmholtz equation, or for that matter solutions of the Laplace equation or the equations of hydrodynamics or elasticity, gives rise to singularities in the integrands of the surface integrals. These mathematical singularities originate from the fundamental solution of the equations used in Green's second identity, whereas the physical problem is perfectly well behaved on the boundaries. Therefore, it is not too surprising that these mathematical singularities can be completely removed by subtracting a related auxiliary solution of the governing equation: ψ(x) in the case of the Helmholtz equation considered here. Apart from having to satisfy the governing differential equation and some mild constraints, see (3.5), there is considerable flexibility in choosing the precise functional form of the auxiliary solution or any free parameters that it may contain—such as the value of x_d in (3.7). A broad physical interpretation is that we have constructed a pressure field, ψ(x), that cancels the value of p and ∂p/∂nat x₀.

One well-studied uniqueness issue is that the solution of the boundary integral method, (3.1), with a closed boundary is identical to that obtained from solving directly the differential form of the Helmholtz equation: ∇²p+k²p=0, except at a discrete set of values of k. These k values correspond to the resonant wavenumbers of the closed boundary [3,4]. The solution of the BRIEF, (3.6), also shares this feature. The resonant or homogeneous solutions corresponding to the resonant values (eigenvalues) of the closed boundary will always emerge from solving the integral equation. But as we shall see in §5, the higher precision that the BRIEF can attain will alleviate much of the practical numerical difficulty.

4. Accurate evaluation at field points near boundaries

The BRIEF of the solution of the Helmholtz equation also offers an accurate and numerically robust method to calculate the solution at field points close to boundaries. Indeed, the loss of precision owing to the near singular behaviour in the evaluation of the requisite integrals that arise in the conventional boundary integral method (CBIM) is often more difficult to deal with than the singularities on the boundaries.

To evaluate the solution p(x_p) of the Helmholtz equation at a point x_p in the solution domain, we first use the conventional boundary integral equation to find [p(x)−ψ(x)] at x_p, with x₀ in (3.4) taken to be the node on the boundary closest to x_p

$$\begin{array}{rl}4\pi p({\mathbf{\text{x}}}_p)& =4\pi \psi ({\mathbf{\text{x}}}_p)-{\int }_{S+S_{\mathrm{\infty }}}[p(\mathbf{\text{x}})-\psi (\mathbf{\text{x}})]\frac{\mathrm{\partial }G({\mathbf{\text{x}}}_p,\mathbf{\text{x}})}{\mathrm{\partial }n}\hspace{0.17em}\text{d}S(\mathbf{\text{x}})+{\int }_{S+S_{\mathrm{\infty }}}\frac{\mathrm{\partial }[p(\mathbf{\text{x}})-\psi (\mathbf{\text{x}})]}{\mathrm{\partial }n}G({\mathbf{\text{x}}}_p,\mathbf{\text{x}})\hspace{0.17em}\text{d}S(\mathbf{\text{x}}).\end{array}$$4.1

The near singular behaviour of the surface integrals owing to G and ∂G/∂n as x_papproaches the boundary, S, can be now removed by subtracting the boundary regularized boundary integral equation (3.6) from (4.1). Then using (3.4) for ψ(x), we obtain a numerically robust expression for p(x_p) that is free of near singular behaviour irrespective of the distance from x_p to any boundary:

$$\begin{array}{rl}4\pi p({\mathbf{\text{x}}}_p)& =4\pi [p({\mathbf{\text{x}}}_0)g({\mathbf{\text{x}}}_p)+{\left(\frac{\mathrm{\partial }p}{\mathrm{\partial }n}\right)}_{0}f({\mathbf{\text{x}}}_p)]\\ & -{\int }_{S+S_{\mathrm{\infty }}}[p(\mathbf{\text{x}})-p({\mathbf{\text{x}}}_0)g(\mathbf{\text{x}})-{\left(\frac{\mathrm{\partial }p}{\mathrm{\partial }n}\right)}_{0}f(\mathbf{\text{x}})][\frac{\mathrm{\partial }G({\mathbf{\text{x}}}_p,\mathbf{\text{x}})}{\mathrm{\partial }n}-\frac{\mathrm{\partial }G({\mathbf{\text{x}}}_0,\mathbf{\text{x}})}{\mathrm{\partial }n}]\text{d}S(\mathbf{\text{x}})\\ & +{\int }_{S+S_{\mathrm{\infty }}}[\frac{\mathrm{\partial }p(\mathbf{\text{x}})}{\mathrm{\partial }n}-p({\mathbf{\text{x}}}_0)\mathrm{\nabla }g(\mathbf{\text{x}})\cdot \mathbf{\text{n}}(\mathbf{\text{x}})-{\left(\frac{\mathrm{\partial }p}{\mathrm{\partial }n}\right)}_0\mathrm{\nabla }f(\mathbf{\text{x}})\cdot \mathbf{\text{n}}(\mathbf{\text{x}})][G({\mathbf{\text{x}}}_p,\mathbf{\text{x}})-G({\mathbf{\text{x}}}_0,\mathbf{\text{x}})]\hspace{0.17em}\text{d}S(\mathbf{\text{x}}).\end{array}$$4.2

5. Examples from acoustics

We use examples drawn from acoustics to illustrate the implementation and advantages of the BRIEF of the Helmholtz equation relative to the CBIM. The first classic problem is the pressure field outside a radiating sphere of radius, R, with a prescribed time harmonic, radially symmetric surface normal velocity $v_n=V\exp (-\mathrm{i}\omega t)$ [3], giving the boundary condition: ∂p/∂n=−iωρV , with k=ω/cwhere ρ is the density and c the speed of sound in the region outside the sphere. This is a Neumann problem because p(x) is to be found for prescribed values of the normal derivative ∂p/∂n. The analytic solution of this problem is [3]: p(r)=−[(iωρV R²)/(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=π/2 for varying numbers of degrees of freedom or unknowns on the surface of the sphere:

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. CBIM Linear c₀: 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₀ in (3.1) is calculated according to the local geometry around each point x₀ [11].

3. CBIM Linear c₀=2π: This an often used approach that assumes the incorrect value of c₀=2π for the solid angle at every node even though the surface tangent plane is clearly not defined at the vertices of the surface elements.

4. BRIEF Constant: Using BRIEF with the assumption that the unknown pressure associated with every flat triangular surface element is constantwithin the element.

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. BRIEF Quadratic: Using BRIEF where the points on the surface of the sphere on which the pressure is to be found are at the nodes of quadratic 6-noded triangular area elements with quadratic shape functions used to represent the boundary.

It is clear from figure 1b that results obtained from the BRIEF are at least as good as, and in the case of BRIEF Quadratic far superior to, any results from the CBIM. Depending on the degrees of freedom, the relative error is smaller by 1–3 orders of magnitude; or BRIEF Quadratic can achieve the same precision as any CBIM with about one-tenth the number of degrees of freedom. In all BRIEF approaches, there is no need to deal with singular integrals or to compute the solid angle c₀ at any node. Moreover, the poor performance of the erroneous CBIM Linear c₀=2π approach is demonstrated.

In figure 2a, we compare the BRIEF using constant, linear and quadratic elements for a single radiating sphere considered at wavenumbers very close to the first resonant value, kR=π [3]. We observe that: (i) using the BRIEF with quadratic (BRIEF Quadratic) elements, the maximum relative error is not significant until |kR−π|<10⁻⁴ and (ii) when the BRIEF is used with constant elements (BRIEF Constant), the resonant wavenumber is located incorrectly. This is because the spherical shape is not well represented by constant elements and resulted in a slightly different resonant wavenumber, even though the polyhedra that represent the sphere have the same volume. Even in this case, the resonant solution is only evident when k is within 1% of the (incorrect) resonant value.

In table 1, we compare the condition number at various values of kR. We see that using BRIEF Linear or BRIEF Quadratic, the condition number at kR very close to the resonant value can be reduced by a factor of 60 when compared with BRIEF Constant. This demonstrates the numerical stability that can be achieved using BRIEF Linear or BRIEF Quadratic.

Next, we study the error when BRIEF Quadratic is combined with one CHIEF point located at r_cf inside the sphere. With kR=2π the lowest resonant mode has a node at r=R/2 [3]. The maximum, mean and median errors when the CHIEF point r_cf is positioned close to this node at R/2 are compared in figure 2b. Note that to incur a maximum error exceeding 1%, r_cf/R has to be within 0.004 of the node at $r/R=\frac{1}{2}$. Thus in practice, the BRIEF Quadratic approach plus one CHIEF point is very unlikely to encounter problems with the resonant solution.