Determination of scattering frequencies for two-dimensional acoustic problems using boundary element method

Abstract

This paper presents a numerical approach for the determination of scattering frequencies of two-dimensional acoustic exterior resonance problems using boundary element method combined with a contour integral method. Since the solving range is surrounded by a Jordan curve in the complex plane, the complex eigenvalues can be acquired through evaluating the contour integral along the Jordan curve using the contour integral method. Complex scattering frequencies are directly extracted by the proposed approach, and their corresponding eigenmodes are also presented by evaluating the vibration amplitudes of internal points. The behaviors of the fictitious eigenfrequencies for unbounded domain are also studied. Numerical examples demonstrate the effectiveness of the approach for the resonance problems of unbounded domains.

Introduction

When an incident wave is propagating through the interfaces of two different mediums or encounter a rigid obstacle as shown in Figure 1, the maximum amplitude of the scattering wave can be excited at the real frequencies located in the vicinity of the scattering frequencies. The determination of the scattering frequencies and the corresponding scattering modes provides a convenient approach for the study of the dynamic stability of underwater or floating body such as ships, submarine, submerged constructions, etc. The peaks of the response curve of the excited body by waves imply that the maximum energy is transmitted to the body. The approach can also be applied to the sound resonances of the instruments and light scattering of suspended particles in liquid. Scattering frequencies can be obtained by eigenvalue analysis in the complex plane, since they are relevant to the singularity of scattering matrix. ¹

Figure 1.

The scattering wave of an obstacle.

The amplitudes of scattering waves by a cylindrical scatterer have poles which are related to zeros of $H_{\nu }^{(1)}(z)$ and its derivative, where $H_{\nu }^{(1)}(z)$ is the Hankel function of the first kind of order ν. Looking into the literature, the computation of the zeros of Hankel functions and its derivative with respect to the argument z has been investigated^2–4 because of the importance of the poles and their trajectories. These complex poles are called scattering frequencies, whose real and imaginary parts are related to the energy peaks of the steady-state solution and the time-decay of the solution, respectively. In acoustic scattering, the zeros of $H_{\nu }^{(1)}(z)$ are the poles of the scattering wave by a soft cylinder. However, the poles of the scattering wave by a rigid cylinder are determined by the zeros of $H_{\nu }^{\prime (1)}(z)$ , where $()\prime$ denotes $\partial /\partial z$ . For a cylinder with a surface impedance Z, the poles of the scattering wave are given as the zeros of $H_{\nu }^{\prime (1)}(z)+\text{i}Z{H}_{\nu }^{(1)}(z)$ , where i is the imaginary unit.

Among the former investigations of scattering problems, one comes across the previous works presents an exhaustive study of the scattering properties of acoustically soft and hard bodies with geometrically-simple shapes. ⁵ Important analysis methods and numerical techniques frequently used for electromagnetic wave propagating problems can be found in the literature. ⁶ T-matrix^7,8 approach has been widely used in acoustic and electromagnetic scattering. Multiple scattering problems are studied using T-matrix method for acoustic, electromagnetic, and elastic fields. ⁹ The conformal mapping approach ¹⁰ is developed for both small and large deviations from a circular cross-section and is more efficient than T-matrix method under certain conditions. The partial-wave series expansion (PWSE) method ¹¹ in cylindrical coordinates is adopted to analyze acoustic scattering of a cylindrical quasi-Gaussian beam, and the results are compared with the ones obtained by T-matrix. PWSE^11,12 is also applied to the investigation of the radiation forces and torque exerted on a rigid elliptical cylinder. Hybrid smoothed finite element method ¹³ is developed for acoustic scattering unbounded domains. The analytical solution ¹⁴ and semianalytical solution ¹⁵ for the scattering on elliptical cylinders are derived, respectively. Particle pairs with corrugated surfaces are considered in electromagnetic problems ¹⁵ and a direct analogy with acoustic counterpart is discussed. ¹⁶

To simulate an unbounded domain with the finite element method, appropriate boundary conditions have to be imposed on the artificial boundary, which is considered as transparent boundaries.^13,17–19 Otherwise, infinite elements ²⁰ have to be employed. Both of the techniques aim at replacing the radiation condition at infinity. Besides the internal discretization, extra parameters for above techniques have to be chosen carefully to keep good results.

In view of the aforementioned limitation, the boundary element method (BEM) is undertaken to extend the modal analysis to the exterior problems. The numerical model with boundary-only discretization offered by the BEM has been widely applied to a variety of exterior physical problems.^21–23 Furthermore, the adoption of the fundamental solutions makes the BEM satisfy the radiation condition at infinity automatically. For scattering problems, boundary elements can be employed to simulate arbitrary geometrical cross-section with both smooth and corrugated surfaces, and the accuracy can be improved by adopting higher order elements. For eigenfrequency computations, both interior and exterior problems result in nonlinear eigenvalue problems, since the parameter ω (circular frequency) or k (wave number) is involved in the each element of the system matrix from the boundary integral equations. The nonlinearity comes also from the fundamental solution which is in the form of Hankel functions. The direct search method, ²⁴ however, is a low efficient method which searches the zero points of the determinant of the system matrix from BEM. For complex eigenfrequencies originated in scattering problems, it requires a two-dimensional (2D) search on the complex plane so that the computing cost increases significantly. Recently, it has been shown that the BEM stands out in numerical accuracy for the computation of eigenfrequencies.^25–28 In particular, comparing with the finite element method, the BEM yields more accurate eigensolutions corresponding to complex eigenmodes due to the exact interpolation in the internal domain. ²⁹ The domain numerical interpolation may result in an incorrect order of eigenfrequencies, as the mode shape in internal domain is relatively complex comparing with it on the boundary, which means BEM has better results than domain interpolation numerical methods for eigenfrequency analysis. Furthermore, although the BEM takes more computing cost than finite element method or other domain interpolation methods for common static or dynamic problems due to the full system matrix, the combination of BEM and the block Sakurai-Sugiura method (BSSM)^30,31 can lead to a parallel computing program framework. Each thread requires a dynamic analysis process but not an eigenvalue analysis, which makes adoption of the fast algorithms possible. The eigenvalue analysis, actually, is carried out in a reduced eigenspace derived using BSSM, and the dimensions of the reduced eigenspace can be chosen depending on situations. The general efficacy of the BEM and BSSM has posed the BEM as a possible alternative to the conventional numerical techniques for modal analysis.^26,32 Specifically, since the scattering frequencies are complex numbers with relatively large imaginary parts, they locate in the complex plane (each quadrant) but not the real axis. The contour integral method such as the BSSM is suitable for extracting those complex eigenvalues off the real axis through defining a positive-oriented Jordan curve.

Comparing with the conventional methods, the proposed approach is easier to be applied for engineering analysis with more complex analysis of objects such as aircrafts, vehicles, etc. With the development of isogeometric analysis, ³³ and its application to BEM, ³⁴ the modeling of analysis of objects becomes simpler. Therefore, the proposed approach has a potential to be an efficient alternative tool for scattering problems both in acoustic and electromagnetic fields with large and complex models. In this paper, the determination of scattering frequencies for exterior problems is carried out using the BEM, and exterior eigenmodes of cylindrical and square scatterers are investigated. The scattering frequencies of a cylindrical scatterer obtained by the proposed approach show a good agreement with the closed-form solutions. The accuracy and effectivity of the proposed method for solving scattering frequencies of unbounded domains are demonstrated through numerical examples.

Boundary integral equations

Let us consider an exterior domain Ω with boundary S as shown in Figure 2. We suppose that S is piecewise smooth and encloses a simply connected domain B. In this work, acoustic problems of the propagation of the sound are considered, and the sound pressure $p(\mathbf{x})$ is governed by the Helmholtz equation

$${\mathrm{▽}}^{2}p+k^{2}p=0,\mathbf{x}\in \Omega$$

(1)

where $k=\omega /C$ , ω is the circular frequency, C is the sound speed, and p is the total wave that includes the incident wave p_i and scattering wave p_s as

$$p=p_i+p_s$$

(2)

Figure 2.

An exterior domain Ω with a closed smooth boundary S.

For governing equation given in equation (1), usually the following boundary conditions are considered

$$p=\overline{p},\mathbf{x}\in S_p$$

(3)

$$q=\frac{\partial p}{\partial n}=\text{i}\rho \omega \overline{v}=\overline{q},\mathbf{x}\in S_q$$

(4)

where n is the outward unit normal vector at x on S, $\partial p/\partial n$ denotes the outward normal derivative. i denotes the imaginary unit, ρ is the density of the medium, and v is the particle velocity. $\overline{()}$ denotes the known boundary conditions.

Since the scattering wave is radiated from the surface of the domain B to the infinity, the sound pressure p_s must satisfy the Sommerfeld radiation condition

$$\underset{\left|\mathbf{r}\right|\to +\infty }{\lim }|\mathbf{r}{|}^{\frac{\tau -1}{2}}\left(\frac{\partial p(\mathbf{r})}{\partial \left|\mathbf{r}\right|}-\text{i}\mathit{kp}(\mathbf{r})\right)=0$$

(5)

where τ = 2 for 2D case, and $\left|\mathbf{r}\right|$ is the distance between the point r in Ω and a fixed origin.

Green’s function for the 2D case is given by

$$\Phi (\mathbf{x},\mathbf{y})=\frac{1}{4}\text{i}{H}_0^{(1)}(\mathit{kR})$$

(6)

where $R=|\mathbf{x}-\mathbf{y}|$ is the distance between x and y and $H_0^1$ is zero order Hankel function of the first kind.

The boundary integral equation can be formulated with Green’s function (equation (6)) as follows

$$\begin{array}{l}\frac{1}{2}p(\mathbf{x})+{\displaystyle {\int }_S\frac{\partial \Phi (\mathbf{x},\mathbf{y})}{\partial n(\mathbf{y})}}p(\mathbf{y})\text{d}S(\mathbf{y})-{\displaystyle {\int }_S\Phi }(\mathbf{x},\mathbf{y})q(\mathbf{y})\text{d}S(\mathbf{y})=0\hfill \end{array}$$

(7)

where $n(\mathbf{y})$ denotes the outward normal to the boundary S at point y. In order to investigate the behaviors of the complex eigenfrequencies and fictitious eigenfrequencies, Burton–Miller’s method is adopted, and the normal derivative of the boundary integral equation (7) is employed as

$$\begin{array}{l}\frac{1}{2}\frac{\partial p(\mathbf{x})}{\partial n(\mathbf{x})}+{\displaystyle {\int }_S\frac{{\partial }^2\Phi (\mathbf{x},\mathbf{y})}{\partial n(\mathbf{y})\partial n(\mathbf{x})}}p(\mathbf{y})\text{d}S(\mathbf{y})-{\displaystyle {\int }_S\frac{\partial \Phi (\mathbf{x},\mathbf{y})}{\partial n(\mathbf{x})}}q(\mathbf{y})\text{d}S(\mathbf{y})=0\hfill \end{array}$$

(8)

With equations (7) and (8), the boundary integral equation: equation (7)+ α equation (8) for Burton–Miller’s method ³⁵ is formed by a manually specified parameter α. Substituting the homogeneous boundary conditions $\overline{p}=0$ or $\overline{q}=0$ into the system matrices of the BEM, a nonlinear eigenvalue problem with a full matrix A is obtained as

$$\mathbf{A}(\omega )X=0$$

(9)

where the elements in A are related to integral with the kernel function as $\Phi$ or $\partial \Phi /\partial n$ . If the incident wave p_i is considered, the right-hand side of equation (9) is not zero. The determination of the scattering frequencies can be carried out by extracting complex eigenfrequencies of equation (9).

Extraction of complex eigenvalues using the BSSM

Due to nonlinearity of the eigenvalue problem provided by equation (9), in this paper, the BSSM is employed for extracting the scattering frequencies distributed in the complex plane. The procedure of utilizing the BSSM for equation (9) is shown in Figure 3. It should be noted that the BEM yields fictitious eigenfrequencies for modal analysis. In previous studies, ³⁶ interior problems produce both real and complex fictitious eigensolutions. For interior problems with linear mediums without damping, there are only real eigenfrequencies considered as the true solutions. The searching domain is a two-dimensional circular area surrounded by a closed Jordan curve defined in BSSM. Therefore, the complex eigenvalues can also be obtained using BSSM as an eigensolver. However, with the adoption of Burton–Miller’s method, the imaginary part of the complex solutions is shifted by adding an imaginary number, while the real eigenvalues (considered as true solutions) remain unchanged. In the beginning, the Burton–Miller’s method is developed for overcoming the difficulties of non-uniqueness associated with fictitious eigenfrequencies resulted by a single boundary integral equation for exterior problems and the method combines the conventional boundary integral equation (so-called conventional BEM) with normal derivative boundary integral equation to ensure the uniqueness. α which is a parameter for Burton–Miller’s method is investigated in the study by Zheng et al. ³⁷ wherein the fictitious eigenfrequencies move along with the changing of α. It is found that the fictitious eigenfrequencies actually just shift off the real axis. However, in this paper, the exterior problems result in fictitious solutions as real eigenfrequencies and the desired scattering frequencies are usually complex numbers. To identify the fictitious eigenfrequencies, it is necessary to utilize Burton–Miller’s method to shift the fictitious ones. In the literature, ²⁰ resonances of a square with a small opening happen both in the interior area and the exterior surface. Thus, in that case, both the real and complex eigenfrequencies are true solutions. Nevertheless, it is another issue out of the scope of this paper. The identification of fictitious solutions is usually carried out by observing the imaginary parts of the results for a simple connected domain. Because of the numerical approaches adopted in the BEM and BSSM, a manually chosen threshold varies for different models. After the extraction of the eigenfrequencies, the eigenvectors are recovered from the projected eigenspace to the original one, and for exterior problems, the eigenvectors that only contain the quantities on the boundary of the BEM models can be extended to the exterior domain using integral equations.

Figure 3.

A general procedure for utilizing BSSM in modal analysis by BEM. BSSM: block Sakurai-Sugiura method; BEM: boundary element method.

In this study, the conventional circular integral path is first adopted to extract the eigenvalues which contain the fictitious solutions. The fictitious eigenvalues for exterior resonance problems are real numbers resulted by the boundary integral equations. The behavior of the fictitious solutions are also observed as Burton–Miller’s method is applied and the similar phenomena are found like the one in the study by Gao et al. ²⁵ The parameters of the circular path is defined as follows, ρ is the radius of the circle, γ is the geometric center of the integral path (center of the circle for circular path). The moment matrices are given as follows

$${\mathbf{M}}_m=\frac{1}{2\pi \text{i}}{\displaystyle {\int }_{\Gamma }\mathbf{T}}(z)z^m\text{d}z$$

(10)

where Γ is a positively oriented closed Jordan curve in the complex plane for ω, then two Hankel matrices ${\mathbf{H}}_{\mathit{Kl}}^{<}$ and ${\mathbf{H}}_{\mathit{Kl}}$ can be formed by using the moment matrices ${\mathbf{M}}_m$ , where m varies from 0 to $2K-1$ :

$${\mathbf{H}}_{\mathit{Kl}}=\left(\begin{array}{cccc}{\mathbf{M}}_0& {\mathbf{M}}_1& \cdots & {\mathbf{M}}_{K-1}\\ {\mathbf{M}}_1& {\mathbf{M}}_2& \cdots & {\mathbf{M}}_K\\ ⋮& ⋮& \ddots & ⋮\\ {\mathbf{M}}_{K-1}& {\mathbf{M}}_K& \cdots & {\mathbf{M}}_{2K-2}\end{array}\right)$$

(11)

$${\mathbf{H}}_{\mathit{Kl}}^{<}=\left(\begin{array}{cccc}{\mathbf{M}}_1& {\mathbf{M}}_2& \cdots & {\mathbf{M}}_K\\ {\mathbf{M}}_2& {\mathbf{M}}_3& \cdots & {\mathbf{M}}_{K+1}\\ ⋮& ⋮& \ddots & ⋮\\ {\mathbf{M}}_K& {\mathbf{M}}_{K+1}& \cdots & {\mathbf{M}}_{2K-1}\end{array}\right)$$

(12)

By solving the eigenvalues of the generalized eigenvalue problem ${\mathbf{H}}_{\mathit{Kl}}^{<}-\omega {\mathbf{H}}_{\mathit{Kl}}$

$${\mathbf{H}}_{\mathit{Kl}}^{<}\mathbf{w}=\omega {\mathbf{H}}_{\mathit{Kl}}\mathbf{w}$$

(13)

the eigenvalues ${\omega }_1,{\omega }_2,\cdots ,{\omega }_K$ located in the closed curve Γ can be obtained.

It can be proved that ${\omega }_1,{\omega }_2,\cdots ,{\omega }_K$ are the eigenvalues of the original problem equation (9) located inside Γ.

Let $({\omega }_j,{\mathbf{w}}_j)$ be the eigenpairs of the generalized eigenvalue problem (equation (13)) and let

$${\mathbf{S}}_m=\frac{1}{2\pi \text{i}}{\displaystyle {\int }_{\Gamma }{z}^m}\mathbf{A}{(\mathbf{z})}^{-1}\mathbf{V}\text{d}z,m=0,1,\cdots ,K-1$$

(14)

where ${\mathbf{S}}_m$ are Ne × l matrix, with which one can form the matrix $\mathbf{S}=({\mathbf{S}}_0,{\mathbf{S}}_1,\cdots ,{\mathbf{S}}_{M-1})$ , which is a Ne × Kl matrix. Then, the eigenvectors can be solved by using the formula

$${\mathbf{x}}_j=\mathbf{S}{\mathbf{w}}_j$$

(15)

where ${\mathbf{x}}_j$ are eigenvectors of the original nonlinear eigenvalue problem. The proofs of the theorems are given in Asakura et al. ³⁰

The contour integrals in equations (10) and (14) are carried out numerically using N-points trapezoidal rule. A circular integration path is defined as $\Gamma =\gamma +\rho e^{\text{i}\theta }(0\le \theta <2\pi )$ and the collocation points are $p_j=\gamma +\rho e^{2\pi \text{i}(j+1/2)/N},(j=0,1,2,\dots ,N-1)$ . ${\mathbf{M}}_m$ and ${\mathbf{S}}_m$ are calculated numerically

$${\mathbf{M}}_m\approx {\widehat{\mathbf{M}}}_m=\frac{1}{N}{\displaystyle \sum _{j=0}^{N-1}{\left(\frac{p_j-\gamma }{\rho }\right)}^{m+1}}\mathbf{T}(p_j)$$

(16)

$${\mathbf{S}}_m\approx {\widehat{\mathbf{S}}}_m=\frac{1}{N}{\displaystyle \sum _{j=0}^{N-1}{\left(\frac{p_j-\gamma }{\rho }\right)}^{m+1}}\mathbf{A}(p_j)\mathbf{V}$$

(17)

As the number of eigenvalues located inside Γ is not known in advance, the number of the eigenvalues has to be determined. To this end, the singular value decomposition (SVD) of the Hankel matrix is performed. After the Hankel matrices are formed, the SVD of ${\mathbf{H}}_{\mathit{Kl}}$ is carried out as follows

$${\mathbf{H}}_{\mathit{Kl}}=\mathbf{C}\mathbf{\Sigma }{\mathbf{E}}^H$$

(18)

where C is a Kl × Kl complex unitary matrix, $\mathbf{\Sigma }$ is a Kl × Kl diagonal matrix with nonnegative real numbers (singular values of ${\mathbf{H}}_{\mathit{KL}}$ ) on the diagonal, and ${\mathbf{E}}^H$ is the conjugate transpose of E which is a Kl × Kl complex unitary matrix. The $\mathbf{\Sigma }$ matrix can be written as

$$\mathbf{\Sigma }=\left(\begin{array}{ccccccc}{\sigma }_1& \hspace{0.17em}& \hspace{0.17em}& \hspace{0.17em}& \hspace{0.17em}& \hspace{0.17em}& \hspace{0.17em}\\ \hspace{0.17em}& {\sigma }_2& \hspace{0.17em}& \hspace{0.17em}& \hspace{0.17em}& \hspace{0.17em}& \hspace{0.17em}\\ \hspace{0.17em}& \hspace{0.17em}& \ddots & \hspace{0.17em}& \hspace{0.17em}& \hspace{0.17em}& \hspace{0.17em}\\ \hspace{0.17em}& \hspace{0.17em}& \hspace{0.17em}& {\sigma }_{m^{\prime }}& \hspace{0.17em}& \hspace{0.17em}& \hspace{0.17em}\\ \hspace{0.17em}& \hspace{0.17em}& \hspace{0.17em}& \hspace{0.17em}& {\sigma }_{m^{\prime }+1}& \hspace{0.17em}& \hspace{0.17em}\\ \hspace{0.17em}& \hspace{0.17em}& \hspace{0.17em}& \hspace{0.17em}& \hspace{0.17em}& \ddots & \hspace{0.17em}\\ \hspace{0.17em}& \hspace{0.17em}& \hspace{0.17em}& \hspace{0.17em}& \hspace{0.17em}& \hspace{0.17em}& {\sigma }_{\mathit{Kl}}\end{array}\right)$$

(19)

where ${\sigma }_1,{\sigma }_2,\dots ,{\sigma }_{m^{\prime }},{\sigma }_{m^{\prime }+1},\dots ,{\sigma }_{\mathit{Kl}}$ are the singular values of the Hankel matrix and are defined as positive. The singular values are listed in decreasing order ${\sigma }_1>{\sigma }_2,\dots ,>{\sigma }_{\mathit{Kl}}$ . Let δ be a positive threshold value and omit the small singular values ${\sigma }_{m^{\prime }+1}<\delta ·{\sigma }_1$ , then a diagonal $m^{\prime }\times m^{\prime }$ matrix is obtained as follows

$${\mathbf{\Sigma }}_{m^{\prime }}=\left(\begin{array}{cccc}{\sigma }_1& \hspace{0.17em}& \hspace{0.17em}& \hspace{0.17em}\\ \hspace{0.17em}& {\sigma }_2& \hspace{0.17em}& \hspace{0.17em}\\ \hspace{0.17em}& \hspace{0.17em}& \ddots & \hspace{0.17em}\\ \hspace{0.17em}& \hspace{0.17em}& \hspace{0.17em}& {\sigma }_{m^{\prime }}\end{array}\right)$$

(20)

Let $\mathbf{H}={\mathbf{C}}^H{\mathbf{H}}_{\mathit{Kl}}^{<}\mathbf{E}$ , then equation (13) becomes the form as

$$\mathbf{H}\overline{\mathbf{w}}=k\mathbf{\Sigma }\overline{\mathbf{w}}$$

(21)

where $\overline{\mathbf{w}}=\hspace{0.17em}{\mathbf{E}}^H\mathbf{w}$

However, the previous study shows that the threshold δ is changing for different numerical models. In order to improve the stability of the algorithm, the threshold ${\delta }^{\prime }$ is introduced and we consider the following modified singular values

$${\mathbf{\Sigma }}^{\mathbf{\prime }}=\left(\begin{array}{ccccccc}{{\sigma }^{\prime }}_1& \hspace{0.17em}& \hspace{0.17em}& \hspace{0.17em}& \hspace{0.17em}& \hspace{0.17em}& \hspace{0.17em}\\ \hspace{0.17em}& {{\sigma }^{\prime }}_2& \hspace{0.17em}& \hspace{0.17em}& \hspace{0.17em}& \hspace{0.17em}& \hspace{0.17em}\\ \hspace{0.17em}& \hspace{0.17em}& \ddots & \hspace{0.17em}& \hspace{0.17em}& \hspace{0.17em}& \hspace{0.17em}\\ \hspace{0.17em}& \hspace{0.17em}& \hspace{0.17em}& {{\sigma }^{\prime }}_{m^{\prime }}& \hspace{0.17em}& \hspace{0.17em}& \hspace{0.17em}\\ \hspace{0.17em}& \hspace{0.17em}& \hspace{0.17em}& \hspace{0.17em}& {{\sigma }^{\prime }}_{m^{\prime }+1}& \hspace{0.17em}& \hspace{0.17em}\\ \hspace{0.17em}& \hspace{0.17em}& \hspace{0.17em}& \hspace{0.17em}& \hspace{0.17em}& \ddots & \hspace{0.17em}\\ \hspace{0.17em}& \hspace{0.17em}& \hspace{0.17em}& \hspace{0.17em}& \hspace{0.17em}& \hspace{0.17em}& {{\sigma }^{\prime }}_{\mathit{Kl}-1}\end{array}\right)$$

(22)

where

$${{\sigma }^{\prime }}_i=\lg ({\sigma }_i-{\sigma }_{i-1})$$

(23)

if ${{\sigma }^{\prime }}_{m^{\prime }}>{\delta }^{\prime }$ , then we can have equation (20). With the threshold ${\delta }^{\prime }$ , the detection of the gap between the singular values is carried out. The parameter ${\delta }^{\prime }$ , in fact, denotes the span between the desired singular values and small singular values (neglectable singular values). The number of desired singular values is, in fact, the dimension of the reduced eigenspace defined by the closed Jordan curve. For different numerical models such as 2D, 3D or different types of elements such as constant elements, linear elements, etc., ${\delta }^{\prime }$ is not unique even when the parameter setting in BSSM remain unchanged. The behavior of ${\delta }^{\prime }$ is investigated in our previous work^25,32 for 3D and 2D cases, respectively. A separation between the two groups of singular values appears with the increasing of the number of collocation points. The singular values corresponding to the eigenfrequencies locating in the Jordan curve becomes stable and the ones corresponding to those eigenfrequencies out of the Jordan curve tend to decrease. A method of determining the gap between the two groups of singular values is proposed in the study by Zheng et al. ²⁶ wherein a strategy of increasing N is discussed and employed to detect the gap.

Numerical examples

In this section, examples are given to show the accuracy and effectiveness of the method for the determination of the scattering frequencies. In the first example, a circular cylinder is considered as a scatterer with a soft boundary/hard boundary and the scattering frequencies are related to the zeros of Hankel functions in the complex plane. In the second example, a square scatter with a hybrid boundary condition is given. The sound speed of air is assumed to be C = 333.1306 (m/s).

Cylindrical scatterers

Let us see the scatterer depicted in Figure 4 where the radius $\overline{R}=1$ m. The exterior problem is defined in the analyzed domain which is the complementary of the cylindrical scatterer. The boundary of the cylinder is considered to be acoustically soft. The scattering frequencies are determined by the following equations

$$H_{\nu }^{(1)}(z)=0$$

(24)

and the scattering frequencies of a rigid cylindrical scatterer with the hard boundary is determined by

$$H_{\nu }^{(1{)}^{\prime }}(z)=0$$

(25)

where $\nu =0,1,2,3,\dots$ , and $z=\omega \overline{R}/C$ .

Figure 4.

The cylindrical scatterer.

The discretized boundary element model is shown in Figure 5 where the open circular symbols denote the edge points of boundary elements and the filled circular symbols denote the internal points in the unbounded domain. Five hundred internal points are collocated along the radial direction uniformly as five layers. The amplitudes of the eigenmodes at the internal points are computed using integral equations.

Figure 5.

Boundary element discretization for the circular structure.

The closed-form solutions of equation (24) are listed in Table 1. The zeros of $H_{\nu }^{(1)}(z)$ locate in the complex plane which can be defined in four quadrants using $\arg z$

$$B_j\equiv \left\{z:\frac{1}{2}(j-1)\pi <\arg z\le \frac{1}{2}j\pi \right\},j=0,\pm 1,\pm 2,\dots$$

(26)

Table 1.

The closed-form solutions of the soft cylindrical scatterer.

ν	Zeros of $H_{\nu }^{(1)}$ in $B_{-1},B_{-5},B_{-9}$
2	0.429-1.281i, 0.360-1.005i, 0.329-0.894i
3	1.308-1.682i, 1.153-1.392i, 1.081-1.273i
4	2.204-1.978i, 1.991-1.678i, 1.892-1.553i
4	0.433-2.629i, 0.398-2.338i, 0.382-2.211i
5	3.113-2.219i, 2.855-1.908i, 2.735-1.779i
5	1.304-3.135i, 1.217-2.840i, 1.177-2.710i

The scattering frequencies of the cylindrical scatterer located in $[150,1850]$ Hz ( $\rho =850\gamma =(1000,0)$ ) are presented in the form of wavenumber in Table 2  where the desired results are symbolized using * and others are considered to be fictitious ones according to their small imaginary parts. Table 3 shows the true solutions in Table 2 and their corresponding closed-from. It can be found that the numerical results show a good agreement with the closed-form solutions along with their multiplicities. The multiplicity of these eigenvalues is 2 and the convergency of the numerical method is shown in Figure 6. It can be seen that the numerical results converge fast as the number of boundary elements increases. With the number of boundary elements Ne = 100, the numerical solutions are accurate enough and the relative errors E_r become less than 0.1%. Table 4 provides the eigenfrequencies located in $[1500,3500]$ Hz ( $\rho =1000,\gamma =(2500,0)$ ). Without using Burton–Miller’s method, the real eigenfrequencies that are related to fictitious interior resonances locate on the real axis. The identification of complex results is carried out by comparing the imaginary parts with a threshold $C·\eta$ where $\eta =0.05$ is chosen experientially.

Table 2.

Eigenfrequencies in the forms of wave numbers for a cylindrical scatterer with a soft boundary.

No.	$k_i={\omega }_i/C$	No.	$k_i={\omega }_i/C$
1	1.309-1.683i*	7	3.114-2.220i*
2	1.309-1.683i*	8	3.833 + 1.006E-007i
3	2.205-1.979i*	9	3.833 + 9.704E-008i
4	2.205-1.979i*	10	5.138 + 1.653E-007i
5	2.406 + 3.984E-008i	11	5.138 + 1.649E-007i
6	3.114-2.220i*	12	5.522 + 2.023E-007i

Table 3.

The true eigenfrequencies in the forms of wave numbers for a cylindrical scatterer and their closed forms.

No.	$k_i^{*}$	Closed form
1	1.309-1.683i*	1.308-1.682i
2	1.309-1.683i*	1.308-1.682i
3	2.205-1.979i*	2.204-1.978i
4	2.205-1.979i*	2.204-1.978i
6	3.114-2.220i*	3.113-2.219i
7	3.114-2.220i*	3.113-2.219i

Figure 6.

The convergency of the eigensolutions for the soft cylindrical scatterer.

Table 4.

Eigenfrequencies obtained by conventional BEM.

No.	ωi
1	1711.6642-0.1078i
2	1711.6642-0.1078i
3	1745.4070-412.4791i*
4	1745.4077-412.4823i*
5	1839.8151 + 8.9440E-003i
6	2066.7463-435.4751i*
7	2066.7463-435.4751i*
8	2126.4474-0.2599i
9	2126.4474-0.2599i
10	2338.2588-1.2777E-002i
11	2338.2588-1.2777E-002i
12	2389.1749-456.3895i*
13	2389.1768-456.3958i*
14	2529.1057-0.4857i
15	2529.1071-0.4856i
16	2712.5177-475.6729i*
17	2712.5177-475.6729i*
18	2805.4191-9.0456E-002i
19	2805.4198-9.0459E-002i
20	2884.2453 + 2.2190E-002i
21	2923.4140-0.7904i
22	2923.4140-0.7904i
23	3036.6410-493.6287i*
24	3036.6412-493.6298i*
25	3253.2862-0.2279i
26	3253.2862-0.2279i
27	3311.5435-1.1788i
28	3311.5450-1.1788i
29	3361.4471-510.4941i*
30	3361.4471-510.4941i*
31	3390.7659 + 3.4939E-003i
32	3390.7659 + 3.4939E-003i

BEM: boundary element method.

In Table 5, the fictitious real eigenfrequencies are considered as the eigenfrequencies corresponding to the interior problem with Dirichlet boundary conditions ³⁵ and added by components that have relative large imaginary parts. The phenomena are similar with the previous study by Gao et al. ²⁵ in which a complex domain is considered. As a matter of interest, for the exterior resonant problems, however, the behaviors of the complex eigenfrequencies are different from the fictitious ones. It should be noted that the number of eigenfrequencies in Table 5 is same as that in Table 4, and the values of complex eigenfrequencies remain unchanged. The true eigenfrequencies obtained by conventional BEM and Burton–Miller’s method are listed in Table 6 where only true solutions are extracted. The eigenmodes corresponding several scattering frequencies are recovered to the original eigenspace and shown in Figures 7 to 12.

Table 5.

Eigenfrequencies obtained by Burton–Miller’s method.

No.	ωi
1	1721.3362 + 178.6715i
2	1721.3363 + 178.6715i
3	1745.4058-412.2613i*
4	1745.4066-412.2646i*
5	1849.7712 + 181.6902i
6	2066.7570-435.1241i*
7	2066.7570-435.1241i*
8	2133.9703 + 177.6857i
9	2133.9703 + 177.6857i
10	2346.1287 + 181.7958i
11	2346.1287 + 181.7959i
12	2389.2062-455.8707i*
13	2389.2080-455.8772i*
14	2535.2243 + 176.5791i
15	2535.2256 + 176.5793i
16	2712.5803-474.9498i*
17	2712.5803-474.9498i*
18	2811.9042 + 181.2375i
19	2811.9044 + 181.2375i
20	2890.8284 + 182.5324i
21	2928.5466 + 175.4270i
22	2928.5466 + 175.4270i
23	3036.7474-492.6626i*
24	3036.7479-492.6630i*
25	3258.7929 + 180.4273i
26	3258.7930 + 180.4273i
27	3315.9500 + 174.2338i
28	3315.9513 + 174.2339i
29	3361.6124-509.2431i*
30	3361.6124-509.2431i*
31	3396.4163 + 182.4714i
32	3396.4163 + 182.4714i

Table 6.

The true eigenfrequencies obtained by conventional BEM (ω_i) and Burton–Miller’s method ( $\widetilde{{\omega }_i}$ ) (cylindrical scatterer).

No.	ωi	$\widetilde{{\omega }_i}$
3	1745.4070-412.4791i	1745.4058-412.2613i
4	1745.4077-412.4823i	1745.4066-412.2646i
6	2066.7463-435.4751i	2066.7570-435.1241i
7	2066.7463-435.4751i	2066.7570-435.1241i
12	2389.1749-456.3895i	2389.2062-455.8707i
13	2389.1768-456.3958i	2389.2080-455.8772i
16	2712.5177-475.6729i	2712.5803-474.9498i
17	2712.5177-475.6729i	2712.5803-474.9498i
23	3036.6410-493.6287i	3036.7474-492.6626i
24	3036.6412-493.6298i	3036.7479-492.6630i
29	3361.4471-510.4941i	3361.6124-509.2431i
30	3361.4471-510.4941i	3361.6124-509.2431i

Figure 7.

The eigenmode of exterior resonance at ω₃.

Figure 8.

The eigenmode of exterior resonance at ω₆.

Figure 9.

The eigenmode of exterior resonance at ω₁₂.

Figure 10.

The eigenmode of exterior resonance at ω₁₆.

Figure 11.

The eigenmode of exterior resonance at ω₂₃.

Figure 12.

The eigenmode of exterior resonance at ω₂₉.

The true and fictitious solutions of the eigenfrequencies by BEM and Burton–Miller’s method are presented in Figure 13, where the fictitious real eigenfrequencies are shifted by Burton–Miller’s method by adding positive components with large imaginary parts. The complex eigenfrequencies that are considered as the scattering frequencies, however, are not affected by the Burton–Miller’s method. From the study, it can be seen that the scattering frequencies can be obtained using the BEM through the BSSM. Furthermore, the filtering for fictitious eigenfrequencies usually can be carried out in two ways: shifting the fictitious real results off the real axis or adapting a complex contour integral path to filter out the real results. The scatterer is a rigid body in this example, so that all the real eigenfrequencies of the exterior problems are regarded as the fictitious solutions and not required.

Figure 13.

The eigenfrequencies obtained by BEM (circular symbols) and Burton–Miller’s method (cross symbols).

Square scatterer with hybrid boundary conditions

A square scatterer with hybrid boundary conditions is given in Figure 14, where the Neumann boundary conditions are specified on the left and right sides of the square scatterer and The Dirichlet boundary conditions are specified on top and bottom sides. The known quantities on the boundary are considered as p = 0 and q = 0 for modal analysis. The parameters for BSSM are chosen as ρ = 1000, $\gamma =(2500,0)$ and Ne = 200 which also cover a circular region.

Figure 14.

The square scatterer.

Internal points depicted in Figure 15 are utilized to illustrate the amplitude of eigenmodes of the exterior unbounded domain. It should also be noted that the domain is infinite and more internal points can be distributed. The scattering modes are shown in Figures 16 to 19. The results presented in Table 7 includes the fictitious eigenfrequencies (without *) and true solutions (with *). Comparing the results in Tables 7 and 8 in which the eigenfrequencies are obtained through Burton–Miller’s method, the fictitious eigenfrequencies shift from the real axis to the first quadrant. It can be seen that 10 eigenfrequencies are obtained in Table 7; however, the first three eigenfrequencies ω₁, ω_2, and ω₃ are out of the circular region defined by ρ = 1000, $\gamma =(2500,0)$ . It implies that the there is a possibility for users to acquire extra results which locate near to the Jordan curve but out of the solving region. The inaccuracy of the numerical contour integral is affected by the eigenvalues locating on or near the curve. ³⁸ Nevertheless, the accuracy of the eigensolutions almost remain the same. During the numerical computation, the singular values corresponding to those eigenfrequencies near the Jordan curve do not separate rapidly with the increasing of N. These results can be filtered out simply and artificially by using the defined solving region.

Figure 15.

Boundary element discretization for the square scatterer.

Figure 16.

The eigenmode of exterior resonance at ω₂ in Table 7.

Figure 17.

The eigenmode of exterior resonance at ω₃ in Table 7.

Figure 18.

The eigenmode of exterior resonance at ω₆ in Table 7.

Figure 19.

The eigenmode of exterior resonance at ω₇ in Table 7.

Table 7.

Eigenfrequencies obtained by conventional BEM.

No.	ωi
1	1479.9458-0.1462i
2	1682.9335-641.2778i*
3	1754.7483-760.6904i*
4	2339.9160-0.2965i
5	2340.0531-0.3058i
6	2713.1913-737.3175i*
7	2757.6961-835.9571i*
8	2959.6342-0.8587i
9	3309.1556-0.3069i
10	3309.3949-0.3799i

BEM: boundary element method.

Table 8.

Eigenfrequencies obtained by Burton–Miller’s method.

No.	ωi
1	2340.0602 + 353.7223i
2	2340.4073 + 353.4688i
3	2705.9017-738.4883i*
4	2742.9232-844.3138i*
5	2960.2419 + 346.0454i
6	3309.3292 + 358.9504i
7	3309.4382 + 358.6197i

Discussion

Here, some results of scattering frequency analysis based on BEM and BSSM are presented. Since BEM is employing a fundamental solution satisfying the radiation condition and has advantages on unbounded problems, it is convenient to use BEM to simulate unbounded scattering problems. BSSM is employed to project the eigenspace to relatively reduced small eigenspace and extracts the nonlinear eigenvalues derived by BEM. The cylindrical scatter employed in the numerical models provides the scattering frequencies as the zeros of Hankel functions so that we can compare the numerical results obtained by the proposed the method with the closed-form solutions. The relative errors are less than 0.2% with 40 elements. With the effectiveness shown by the accurate results, the solving region is enlarged. Tables 2 and 4 both show the fictitious eigenfrequencies which are real values. The emergence of fictitious eigenvalues is caused by the single boundary integral equation. To identify the fictitious eigenvalues, Burton–Miller’s method is employed, and Table 5 shows that the fictitious eigenvalues are added by positive imaginary parts and move to the first quadrant. To apply this method to scattering frequency problems, users can directly adopt the combined boundary integral equations of Burton–Miller’s method and filter out the fictitious eigenfrequencies by judging their imaginary parts. The real eigenfrequencies, in scattering problems, are all fictitious ones; however, their imaginary parts are not zero due to the numerical errors but they are still very small. Usually, the scattering frequencies of simple shapes are complex values and can be distinguished just by checking the large imaginary parts (comparing with the small nonzero imaginary parts of fictitious eigenfrequencies). However, the imaginary parts of scattering frequencies tend to be zero as the boundary of the scatterer encloses a cavity with small opens (musical instruments such as flutes, holed wind instruments, etc., sometimes result in such cases). Therefore, Burton–Miller’s method is required to identify the true solutions if the true eigenfrequencies have small imaginary parts approximating to zero.

The eigenfrequencies locating near or on the Jordan curve lead to an unstable detection of the reduced eigenspace. From Table 7, it can be seen that the first three eigenfrequencies are out of the Jordan curve but still obtained. However, the first eigenfrequencies $1479.9458-0.1462\text{i}$ move to the first quadrant with Burton–Miller’s method, and its new location is further away from the Jordan curve so that the detection method for the gap of the singular values works again. Although enclosing the eigenfrequencies close to the solving region into the reduced eigenspace does not decrease the accuracy of the solutions, it is a still a significant problem for the stability of the method. More importantly, the detection of gap in singular values is based on a sufficient N which is directly related to the computation cost. Thus, the estimation of the number of eigenvalues in the reduced eigenspace defined by the Jordan curve makes the proposed method more efficient and it is another issue for the future works.

Conclusion

This work for this paper is carried out with the objective of providing a numerical tool for the computation of scattering frequencies for 2D acoustic problems. The scattering frequencies which are complex numbers are given by unbounded domains. The BEM shows a simple modeling procedure for unbounded problems and the accuracy of the numerical solution is relatively high with a small number of elements. Furthermore, the real eigenvalues which are considered as fictitious solutions for unbounded domains are investigated using Burton–Miller’s method. Similar to the bounded interior problems, the fictitious eigenvalues move away from the real axis and the desired scattering frequencies are kept unchanged. The examples show that the numerical approach can carry out the modal analysis for unbounded problems. The investigator believes that the solutions produced in this work are of very high accuracy and there is no missing link in mode sequences in the results.