Closed-form evaluation of potential integrals in the Boundary Element Method

TL;DR

This paper introduces an analytical method for evaluating singular and near-singular integrals in the Boundary Element Method for the Helmholtz equation, improving accuracy and guiding quadrature rule selection.

Contribution

It develops a closed-form analytical approach for potential integrals, enhancing precision and efficiency in boundary element computations for wave problems.

Findings

Achieves accuracy comparable to machine precision.

Provides a criterion for quadrature rule selection.

Demonstrates effectiveness in Helmholtz equation solutions.

Abstract

A method is presented for the analytical evaluation of the singular and near-singular integrals arising in the Boundary Element Method solution of the Helmholtz equation. An error analysis is presented for the numerical evaluation of such integrals on a plane element, and used to develop a criterion for the selection of quadrature rules. The analytical approach is based on an optimized expansion of the Green's function for the problem, selected to limit the error to some required tolerance. Results are presented showing accuracy to tolerances comparable to machine precision.

Equations

$$\varphi (\mathbf{x})$$

$$G(\mathbf{x};{\mathbf{x}}_1)$$

$$R$$

$$I$$

$$R^2$$

$$r$$

$$R^2$$

$$R_{\text{mid}}^2$$

$$\frac{1}{R}$$

$${ϵ}_Q$$

$$P_q(\cos \varphi )\sim {\left(\frac{2}{\pi q\sin \varphi }\right)}^{1/2}\cos \left[(q+1/2)\varphi -\pi /4\right],$$

$${ϵ}_Q$$

$$=\Im \frac{1+j}{R_{\text{mid}}}\frac{{\mathrm{e}}^{\mathrm{j}\varphi /2}}{[\pi (Q+1)\sin \varphi {]}^{1/2}}\frac{{\left(t\cos \varphi {\mathrm{e}}^{\varphi }\right)}^{Q+1}}{1-t\cos \varphi {\mathrm{e}}^{\mathrm{j}\varphi }}.$$

$$E_Q$$

$$\varphi$$

$$s$$

$${\alpha }^2=z^2/S^2,\overline{R}$$

$$I_0={\mathrm{e}}^{\mathrm{j}k|z|}{I}_0^{\prime },I_x={\mathrm{e}}^{\mathrm{j}k|z|}{I}_x^{\prime },I_y={\mathrm{e}}^{\mathrm{j}k|z|}{I}_y^{\prime },$$

$$I_0^{\prime }$$

$$I_x^{\prime }$$

$$I_y^{\prime }$$

$$\text{with}R$$

$$\frac{\partial I_0}{\partial z}$$

$${\mathrm{e}}^{\mathrm{j}k(R-|z|)}$$

$$e_q$$

$$\text{where}\sin x$$

$$I_0^{\prime }$$

$$I_x^{\prime }$$

$$I_y^{\prime }$$

$$K_{q,0}$$

$$K_{q,x}$$

$$K_{q,y}$$

$$k^q{\int }_0^{\overline{r}}{\left(R-|z|\right)}^q\frac{r}{R}\mathrm{d}r$$

$$k^q{\int }_0^{\overline{r}}{\left(R-|z|\right)}^q\frac{r^2}{R}\mathrm{d}r$$

$$=\frac{sS(kS{)}^q}{q+2}{\left(\frac{\Delta }{\cos \theta }-\alpha \right)}^{q+1}\frac{1}{\cos \theta }+\frac{2|z|}{q+2}{J}_q,$$

$$J_q$$

Full text

1 Introduction

A central part of the Boundary Element Method (BEM) is the evaluation of potential integrals, to compute the contribution of an element to the potential field, or to the entries of the solution matrix. It is thus a key factor in the accuracy and efficiency of any implementation, and one which has attracted great interest over many decades. In this paper we develop a method for the evaluation of integrals which arise in the three-dimensional BEM for the wave equation, in particular in acoustics, where the acoustic potential $\phi$, external to a surface $A$, is given by the integral formulation:

[TABLE]

where $\mathbf{x}$ indicates position, subscript $1$ variables of integration on the surface $A$, and $n$ the outward pointing normal to the surface. The Green’s function $G$ is:

[TABLE]

where $k$ is acoustic wavenumber. Given the surface potential $\phi$ and gradient $\partial\phi/\partial n$, the potential, and, after differentiation, its gradient(s), can be evaluated at any point in the field. Also, given a boundary condition for $\phi$ and/or $\partial\phi/\partial n$ on $A$, the integral equation can be solved for $\phi(\mathbf{x})$ and/or $\partial\phi/\partial n(\mathbf{x})$, $\mathbf{x}\in A$.

If the boundary integral equation is solved using a collocation method, the surface $A$ is divided into elements, here taken to be plane triangles, and suitable shape functions are used to interpolate the potential on these elements. The integral equation is transformed to a linear system in the element potentials, with the influence coefficients determined by the potential generated by each element at each node of the surface mesh. This leads to the requirement to evaluate integrals $I$ and $\partial I/\partial n$ where:

[TABLE]

with $A_{e}$ the surface of an element and $(\xi,\eta)$ a coordinate system local to $A_{e}$. The requirement then is to evaluate integrals of $\exp[\mathrm{j}kR]/R$ and its derivatives over a triangular element. This is especially challenging when the field point $\mathbf{x}$ is on, or near, the element, and the $1/R$ singularity must be accommodated in the integration scheme.

There are numerous numerical schemes for the evaluation of the surface integrals, which mainly vary in their approach to dealing with the singularity. There are also a number of analytical schemes for the equivalent integral in the Laplace equation [4, 8, 7, 3, 10, 12, for example], some of which can be used to deal with the singular terms in the acoustic problem and thus ease numerical integration, but there are few analytical methods for the Helmholtz problem. Clearly, given the absence of an exact analytical solution for the retarded potential from a plane element, any closed-form solution is an approximation, but it should be possible to approximate the integral to any required accuracy, in a form amenable to analytical manipulation, so that the result can be used as if it were an analytical formula for the potential. This is especially important for the case of a field point on or near the element, where the ability to handle singularities analytically offers an advantage over purely numerical schemes.

To the author’s knowledge there are two published methods for closed-form or analytical evaluation of the Helmholtz potential from a planar element [13, 9]. These use two different approaches to the problem. In one [13], an expression is derived in the Fourier domain resulting in an expression based on a series of terms defined by integrals of Hankel functions. These integrals can be evaluated analytically in terms of Struve functions, yielding a closed-form solution for the potential from a planar element, but at the expense of using special functions not routinely available in numerical libraries.

The second approach [9], which is similar in spirit to the method of this paper, makes use of results derived for the Laplace problem [7] and approximates $\exp[\mathrm{j}kR]$ as a polynomial over the element. This is justified by noting that the element size is limited by the requirement to avoid aliasing in the representation of the surface potential, so that a relatively low-order approximation containing five or six terms of the Taylor series for $\exp[\mathrm{j}kR]$ is adequate for evaluation of the integrals to the tolerance specified.

The method of this paper uses a similar approach, in that it replaces the exponential with an approximation of controlled error, and uses results from an analysis of the Laplace problem [4] to compute the terms in the resulting expansion. It differs in the form of Laplace solution used, and in the choice of expansion for the exponential, to give a systematic control on quadrature accuracy optimized to require a minimum number of terms. Additionally, a criterion is provided for choosing when, and when not, to use the analytical approach or a purely numerical method, based on an error analysis of integration using a polar coordinate transformation. To the author’s knowledge, this error analysis is novel and may have applications more generally.

2 Integration of the $1/R$ potential

In order to motivate the development of the closed-form expression for the acoustic potential, we begin by analyzing the numerical evaluation of the Laplace potential, which corresponds to the leading-order, singular, part of the Helmholtz potential, which gives rise to the difficulties in numerical integration.

The model problem is shown in Figure 1 and consists of the evaluation of

[TABLE]

over the area of the triangle shown, which lies in the plane $z=0$, with the usual transformation to polar coordinates $(r,\theta)$ for the integration.

The error in the evaluation of this integral, especially at small values of $z$ arises from the singular, or near-singular, term $1/R$. Here we develop an approximate error analysis for the evaluation of this term, which can be used in determining the required order of integration for $r/R$ or when to switch to some other quadrature approach, such as that in the next section. An error analysis for integration using the polar coordinate transformation has been published previously [11] but the analysis presented here appears to be novel and is simple enough for use as an a priori estimator in determining quadrature order in applications.

The analysis depends on an error estimate for the $1/R$ term in a numerical polar integration, such as (8a). If such an integration is performed using Gaussian quadrature, the integrand is being approximated by a polynomial over the interval of integration and the accuracy of the approximation is determined by the number of terms required to approximate the integrand accurately. We perform the analysis by estimating the error in the polynomial expansion of $1/R$ and use it to give an approximation of the order of polynomial required to approximate $1/R$ to a given tolerance. Given this polynomial order, a Gaussian quadrature of sufficiently high order can be selected, or if the order required is too great, the analytical method of the following section can be used.

From Figure 2, we write

[TABLE]

and neglect the case of $\phi=0$ as in this case $r/R\equiv 1$ and the polynomial representation of the integrand raises no difficulties.

Using the generating function for Legendre polynomials [5, 8.921],

[TABLE]

If the expansion is truncated at $q=Q$, the error at any value of $t$ is given by the remainder

[TABLE]

which can be rewritten using the large-order asymptotic form of the Legendre polynomial [5, 8.918],

[TABLE]

so that

[TABLE]

An upper bound for the sum can be found by replacing $q^{1/2}$ with $(Q+1)^{1/2}$ and, upon rearrangement,

[TABLE]

As will be seen, this is an accurate estimate of the remainder in the polynomial expansion of $1/R$ but it is oscillatory as a function of $q$, so we adopt the more convenient measure of the magnitude rather than the imaginary part,

[TABLE]

We note that (5) could be integrated over $t$ to give an estimate of the total error in the integral of $r/R$ but this gives an unwieldy expression with little advantage in implementations. Instead we adopt as error criterion the absolute value given by (6) with the value of $t$ given by the nearest point on the element. In particular, when the projection of the field point lies on the element, i.e. when the triangle in Figure 1 encloses the origin, $t=1$ and the error estimate for $1/R$ is

[TABLE]

From the form of the error estimate, the reason for the difficulty in evaluating near-singular integrals is clear: near the element plane where $\phi\to 0$, approximation of the integrand by a polynomial, implicit in the use of Gaussian quadratures, incurs a very large error, even for quite high order quadratures with large $Q$.

To minimize the computational burden of using the criterion, it is applied in the following manner. Given the transformation into coordinates based on the element plane, the minimum distance $r_{\min}$ from the origin to the triangle can be determined (see Figure 5) with $r_{\min}\equiv 0$ when the triangle encloses the origin. The maximum distance to a vertex $r_{\max}$ is found similarly. Then we set $r_{\text{mid}}=r_{\max}/2$, $t=(r_{\text{mid}}-r_{\min})/\overline{r}$ and other quantities as above. The criterion is then applied by computing $E_{Q}$ for $Q=1,2,\ldots$ until $E_{Q}$ falls below some specified tolerance, and returning the resulting value of $Q$, the order of polynomial required to compute $1/R$ to the specified tolerance over the range of the integral. We note that the error measure here is the maximum error in $1/R$ at any point in the range of integration, which is quite a stringent, though conservative, measure, but it will be found that $E_{Q}$ is a useful assessment of the accuracy of quadrature.

Figure 3 shows the error estimates as a function of $z$ for a test case with a 32nd order polynomial, equivalent to a 16 point Gaussian quadrature. The error estimate $\epsilon_{Q}$ is seen to be very reliable, and the magnitude $E_{Q}$ does indeed match the amplitude of the error quite closely. Figure 4 shows the error as a function of $Q$ for fixed $z$ and again the error behavior is accurately captured by the estimators. Despite the relative simplicity of the error measures, they give reliable indicators of the accuracy of the quadrature or of the order of quadrature required for a given tolerance. We note finally that the quantities used in the error measure are typically computed as part of the geometric transformations required in generating a quadrature on an element, so that there is very little overhead in applying the error estimate.

3 Analysis

The problem to be considered is evaluation of the Helmholtz single- and double-layer potential integrals on a planar triangular element. Integration is performed after transformation of coordinates such that the triangular element is defined by vertices $(x_{i},y_{i},0)$ and the field point lies at $(0,0,z)$. The triangle is then decomposed into up to three triangles each having a vertex at $(0,0,0)$. The process is shown in Figure 5. The approach is similar to that taken in a previous analysis for the Laplace potential [4], though some changes are required to make it suitable for the Helmholtz problem.

Figure 6 shows the basic triangle which is used for the evaluation of the contributions from the subtriangles of Figure 5. It has one vertex at the origin, i.e. at the projection of the field point onto the element plane, and is defined by the lengths of the two sides which meet at the origin, $r_{1}$ and $r_{2}$, and by the angle $\Theta$ between them.

In developing the analysis, we assume that the triangular element conforms to some reasonable standards of quality, in particular that the edge length is no greater than some specified fraction of a wavelength, typically between one sixth and one eighth. This translates into a limit on $k\ell$, where $\ell$ is a typical edge length. Taking into account the need to deal with triangles which are larger than the element proper, such as triangle $023$ in Figure 5, we assume that $k\ell<\pi/2$, which allows us to limit the size of the expansions which will be employed in evaluating the potential integrals, while retaining the required accuracy. If necessary, this limit can be increased, at the expense of extra computational effort and a small increase in stored data.

3.1 Basic integrals

Integration is performed on the reference triangle of Figure 6, using the polar coordinate system $(r,\theta)$. Geometric parameters are defined,

[TABLE]

and auxiliary variables used in performing the integrations are

[TABLE]

The integrals to be evaluated are the zeroth and first order derivatives with respect to $z$ of

[TABLE]

where the basic integrals are

[TABLE]

and correspond to the zero and first order source terms required for linear shape functions on a plane element. The normal derivatives are given by differentiation with respect to $z$, for example,

[TABLE]

where the element is oriented such that the normal lies in the positive $z$ direction, and the upper (lower) signs are taken for positive (negative) $z$.

The integrals are evaluated by expanding the complex exponential in a polynomial approximation, and evaluating term-by-term using analytical formulae defined by recursion relations, as in previous work [4]. The form of the approximation for $\exp[\mathrm{j}kx]$ will be considered in Section 3.3, but for now we write

[TABLE]

Expanding in powers of $k(R-|z|)$ has the advantages of ensuring that the expansion remains valid for large values of $z$ as $(R-|z|)\to 0$ as $z\to\infty$, and providing a natural reduction in the number of terms in (9) for increasing $z$.

Substituting (9) into (8), yields

[TABLE]

where

[TABLE]

The integrals of (11) can be evaluated analytically using a combination of recursions and tabulated integrals. The inner integrals are given by,

[TABLE]

The integral $J_{q}$ can be evaluated using the recursion

[TABLE]

so that all required terms are written in a form suitable for the application of standard formulae for trigonometric integrals [5, 2.58],

[TABLE]

The normal derivatives of the integrals can be evaluated by differentiating terms, yielding,

[TABLE]

All integrals can then be evaluated using the results of Section A. This gives a means of evaluating all required expressions for the integrals on the triangular element, which can then be summed to give the integral over the initial general triangle.

3.2 Hypersingular integral

The results of (3.1) may be used to solve boundary integral problems using a standard Helmholtz equation. It is often desirable to employ a Burton and Miller approach [2] to avoid the well-known problem of fictitious resonances when the wavenumber $k$ in the exterior problem coincides with an eigenvalue of the interior problem. In this approach, the Helmholtz equation is combined with its normal derivative to yield a formulation which is numerically valid for all real wavenumbers, at the expense of requiring the evaluation of hypersingular integrals of the form

[TABLE]

In order to meet continuity requirements, the collocation points in a hypersingular method must lie strictly within elements, though discontinuous elements offer a way around this [6], and so we give a result for the zero-order (constant) element only:

[TABLE]

3.3 Approximation of exponentials

In order to efficiently evaluate the formulae of Section 3.1, we require a means of selecting the polynomial approximation to the exponential, (9). The most obvious choice is to truncate the Taylor series for $\mathrm{e}^{x}$ at some point where the estimated remainder is smaller than a specified tolerance $\epsilon$. For reasons of efficiency, however, we adopt an “economized” series which replaces the truncated Taylor series with a polynomial approximation with minimum deviation and a minimized error over the range where the polynomial is used. Given that the integral terms are evaluated using recursion relations, by reducing the number of terms, we also reduce the chance of numerical error accumulating in moving from term to term.

The economization algorithm is that given by Acton [1, p291–296] and is used to generate a set of polynomial approximations of $\sin x$ and $\cos x$ over a range $0\leq x<\Delta x$, to a tolerance $\epsilon$ where $\epsilon$ is the maximum difference between $\exp[\mathrm{j}x]$ and the polynomial approximation over the range $0\leq x<\Delta x$. For the calculations of this paper, $\Delta x=\pi/16,\pi/8,\pi/4,\pi/2$, and $\epsilon=10^{-n}$, $n=3,6,9,12,15$. In the implementation, a polynomial approximation is chosen which has the required maximum error less than $\epsilon$ and $k\ell<\Delta x$. In the case of $\Delta x=\pi/2$, for example, this gives a reduction in the number of terms required from fifteen for the truncated Taylor series to eight for the economized polynomial when $\epsilon=10^{-9}$.

3.4 Summary of method

The quadrature method of the previous sections can be summarized as follows, for a triangle $(\mathbf{x}_{1},\mathbf{x}_{2},\mathbf{x}_{3})$ which has been rotated into the plane $z=0$ and field point $\mathbf{x}=(0,0,z)$, Figure 5:

determine the closest and furthest points on the triangle boundary and their radial distances $r_{\min}$ (Figure 5) and $r_{\max}$; 2. 2.

compute the required order of quadrature $Q$ for polynomial approximation of $1/R$, Section 2; 3. 3.

if $Q$ falls below the set limit:

1. (a)

evaluate the integrals numerically and terminate;

otherwise

1. (a)

decompose the triangle into up to three sub-triangles centered at the origin; 2. (b)

for each sub-triangle compute the contribution using the formulae of Section 3.1 and accumulate, taking account of sub-triangle orientation.

4 Numerical testing

As a numerical test of the performance of the method, we use the same test case as Pourahmadian and Mogilevskaya [9], Figure 7. Four points are selected in the element plane, as indicated, and we evaluate $I_{0}$ for $k=1$ as a function of $z$, vertical displacement from the element; results for $I_{x}$ and $I_{y}$ are similar. As a reference for error estimation, we follow Pourahmadian and Mogilevskaya and use a polar transformation and a $50\times 50$ point Gaussian quadrature, which is accurate to eight significant figures [9]. The reported error is

[TABLE]

with superscripts ‘a’ and ‘c’ denoting ‘analytical’ and ‘computed’ values respectively. Error is evaluated by specifying the required tolerance $\epsilon^{(a)}$ in the analytical method and computing the resulting $\epsilon$. A second set of error calculations are presented by fixing $\epsilon^{(a)}=10^{-12}$, varying the order of numerical quadrature in the polar transformation, and computing the resulting error $\epsilon^{(c)}$. Figure 8 gives $\epsilon$ as a function of $z$ for varying $\epsilon^{(a)}$, and $\epsilon^{(c)}$ for varying order of Gaussian quadrature, plotted with $Q$ computed for varying values of $E_{Q}$. Figure 8 shows error data for the evaluation of $I_{0}$ and Figure 9 for the normal derivative $\partial I_{0}/\partial z$.

The left-hand column of Figure 8 shows the error estimate for points 1, 2, 3, and 4 in Figure 7 which correspond respectively to field points whose projections lie on a vertex, in the interior, on an edge, and outside the element. The reference integral $I_{0}^{(c)}$ for points 1–3 is computed using the $50\times 50$ Gaussian quadrature after transformation to polar coordinates, and that for point 4 using the 175 point symmetric quadrature of Wandzura and Xiao [14]. Errors are computed with a requested tolerance $\epsilon^{(a)}=10^{-3,-6,-9,-12}$ and the computed errors reflect both the accuracy of the analytical method and the conformity to the requested tolerance. In applications, there is reasonable confidence that the error will be approximately equal to that requested, without excessive computation.

The right-hand column of Figure 8 presents data relevant to the error estimate $E_{Q}$ and the accuracy of Gaussian quadrature in this problem. The darker curves show an error estimate computed as the difference between the analytical method with a requested tolerance of $10^{-12}$ and numerical quadrature of varying order. As expected the low-order methods, e.g. $4\times 4$ points, give a larger error and the high-order approach, $32\times 32$ points, gives accuracy comparable to the analytical technique, except for small values of $z$. In each case, the breakdown of the polynomial approximation for $1/R$ is apparent in the increase of error as $z\to 0$, most clearly for the $16\times 16$ quadrature where the error increases markedly from $z\approx 0.3$. Of interest here is the value of $E_{Q}$ as a criterion for selecting quadrature rules. The lighter curves show the value of $Q$ found from (6) with varying values of $E_{Q}$. The curves do indeed predict quite well the point at which the polynomial approximation to $1/R$ is no longer accurate and the Gaussian quadrature begins to fail, confirming the reliability of the measure as a criterion for the selection of quadrature rules.

Figure 9 gives similar results but for the evaluation of the normal derivative of the layer potential, also required in BEM calculations. The results are similar to those in Figure 8 and the discussion of those data carries over to here, but it is worth noting that though the error behavior of the Gaussian quadratures is different from that in Figure 8 (compare the results for point 3, for example), the curves of $Q$ still function as a reliable criterion for selecting a quadrature method.

5 Conclusions

An analytical method for the evaluation of potential integrals in boundary element codes for the Helmholtz equation has been presented and tested. An error estimator for purely numerical quadrature has been derived and used to establish a criterion for quadrature method selection. The quadrature method has been tested and found to be accurate and reliable; the error criterion is a reliable technique for quadrature selection. We believe that the quadrature method proposed is a suitable plug-in replacement in BEM codes for the wave equation where an a priori error estimate for element integrals and an economical integration are required.

Appendix A Basic integrals

The evaluation of the potential integrals requires a number of elementary integrals which can be computed using results from standard tables combined with recursions. This appendix contains the results required for the evaluation of the trigonometric integrals of the main paper, written in terms of the parameter $\alpha$, $0\leq\alpha<1$, and $\Delta^{2}=1-\alpha^{2}\sin^{2}\theta$. The results are given as the indefinite integral, with a separate result where necessary for the in-plane case $\alpha=0$.

The first basic term is

[TABLE]

where $s=0,1,2,3$. The terms in the summation are pseudo-elliptic integrals which can be evaluated using elementary functions and recursion relations [5, 2.58].

Using the transformation $u=\tan\theta$ and noting that $\Delta/\cos\theta=(1+{\alpha^{\prime}}^{2}\tan^{2}\theta)^{1/2}$

[TABLE]

with $\alpha^{\prime}=(1-\alpha^{2})^{1/2}$.

The integral term can be evaluated using the recursion

[TABLE]

seeding the recursion with

[TABLE]

and using

[TABLE]

For $\alpha=0$,

[TABLE]

A second, similar, integral is

[TABLE]

which can be seeded with [5, 2.584]

[TABLE]

For $\alpha=0$,

[TABLE]

In an implementation of the method of this paper, when the required geometric parameters have been calculated for the reference triangle, and the appropriate expansion for $\exp[\mathrm{j}kx]$ has been selected, the first step is to compute the required elementary integrals (23) and (33) using the initial values and the recursion relations. The computed terms can then be used in the summations of (3.1) to evaluate the potential integrals.

For convenience, we define

[TABLE]

which are readily evaluated using integration by parts. Then,

[TABLE]

and similarly

[TABLE]