A boundary element eigensolver for acoustic resonances in cavities with impedance boundary conditions Open Access

Analysis of acoustic resonances in cavities can provide information for cavity noise control (Zhang et al., 2019). In general, the finite element method (FEM) is a powerful numerical tool for the resonance analysis of interior cavities. However, FEM is troublesome when dealing with open cavities or cavities pasted with sound absorbing materials (Hein et al., 2007; Moheit and Marburg, 2018; Panneton and Atalla, 1997). The aim of this letter is to develop a numerical scheme for the analysis of acoustic resonances in cavities with impedance boundary conditions (IBCs) based on the boundary element method (BEM). Since the coefficient matrices of the boundary element (BE) system are frequency-dependent, the resultant eigenproblem is nonlinear. To solve the problem effectively, the nonlinear eigenproblem (NEP) is transformed into a normal and small linear one which is tractable with common eigensolvers by using a contour integral approach (Asakura et al., 2009).

Similar to the BE solution for exterior acoustic problems, spurious eigenfrequencies also arise in the BE solution for interior acoustic problems. The emergence of spurious eigenfrequencies could be explained by the relevance of the boundary integral equations for interior and exterior problems. Analytical and numerical discussions have been given by Chen et al. (2003) and Zheng et al. (2019) for multiply-connected and simply-connected domains, respectively. Unlike the spurious eigenfrequencies of BEM for exterior acoustic problems, the spurious ones of BEM for interior acoustic problems are often complex since they are consistent with the eigenfrequencies of the corresponding exterior problems. Due to damping, the true eigenfrequencies of the cavities with absorbing boundary conditions are also complex numbers for which imaginary parts have the same sign as the spurious ones. As a consequence, it is hard to identify the spurious eigenfrequencies directly. In this letter, a variant-parameter scheme, based on the combined boundary integral formulation (Burton and Miller, 1971), is suggested to identify and filter out the spurious eigenfrequencies.

Time-harmonic acoustic problems are considered in this letter. The partial differential equation (PDE) governing the steady-state linear acoustics is the Helmholtz equation

where $∇2$ is the Laplacian operator, $p(x)$ is the complex sound pressure at point x, $k=ω/C=2πf/C$ is the wave number, C is the ambient speed of sound, ω is the angular frequency, f denotes the frequency, Ω represents the acoustic domain bounded by the closed boundary Γ, and d denotes the space dimension. The harmonic time dependence $e−iωt$ is assumed, where $i=−1$ and t is time. As the IBC on Γ is considered, we have

where $n→(x)$ is the normal vector, $ϱm$ is the medium density, $Z(x)$ is the acoustic impedance at point x on Γ, and $β(x)=iϱmω/Z(x)$⁠. It is known that the IBC generalizes the sound-soft and sound-hard boundary conditions to address a big amount of cases between these two extremes.

Equation (1) can also be treated as an eigenvalue PDE to solve for eigenmodes and eigenfrequencies. However, to solve the problem with the BEM, it is recast into the Kirchhoff-Helmholtz boundary integral equation (HBIE) as

where $c(x)$ is a constant depending on the solid angle at x, $c(x)=1$ if $x∈Ω, c(x)=1/2$ if $x∈Γ$ and Γ is smooth at x, $G(x,y)$ denotes the Green's function and is given by $G(x,y)=iH0(1)(kr)/4$ for two-dimensional (2D) problems. Herein, $H0(1)$ is the zeroth order Hankel function of the first kind, and $r=|x−y|$ denotes the Euclidean distance.

Discretizing Γ into a set of BEs and collecting the discretized HBIEs for all boundary nodes, we can obtain the BEM system equation given by

where $B$ is a diagonal matrix given by $B=diag(β1, β2, ⋯, βNd)$⁠, N_d is the number of degrees of freedom (DOF) of the problem, $H$ and $G$ are the system matrices whose entries are expressed as

where δ_ij denotes the Kronecker delta, $δij=1$ when i = j, $δij=0$ when $i≠j$⁠, and $ϕj$ represents the interpolation function for approximation.

Looking at the kernels in the entries of H and G, we can find both H and G are matrix-valued functions with respect to the wave number or frequency implicitly. As a result, it leads to a nonlinear eigensystem

where $A=H−GB$ and $ψj$ is the eigenvector associated with the eigenvalue λ_j. In general, solving such a NEP is not straightforward and easy. Therefore, a number of approaches have been given in the last several decades (Chen et al., 2001, 2002; Güttel and Tisseur, 2017). In Sec. 2.2, an approach based on the contour integral is employed to transform the NEP described by Eq. (6) into a normal and small linear one which is tractable with common eigensolvers such as the lapack routines.

The contour integral approach proposed by Asakura et al. (2009) is applied to convert the NEP described by Eq. (6) into a standard linear one. In this approach, the conversion of the NEP is achieved through two Hankel matrices consisting of the moment matrix

where $C$ denotes a closed Jordan curve in complex frequency domain, $T(z)=PA−1(z)Q, A∈ℂNd×Nd$ is the frequency-dependent system matrix in Eq. (6), $P∈ℂL×Nd$ and $Q∈ℂNd×L$ are nonzero matrices whose entities are set randomly, and L is an integer which should be bigger than the maximum algebraic multiplicity of the eigenvalues located inside $C$⁠. Usually, $Q$ can be specified as the conjugate transpose of P.

Then, the Hankel matrices $H1$ and $H2$ can be defined by $Mℓ$ as

where K is a positive integer requested to ensure $KL⩾Ne$⁠, and N_e is the total amount of eigenvalues lying inside $C$⁠. Usually, special schemes (Di Napoli et al., 2016; Maeda et al., 2011) could be used to estimate N_e in advance. The multiplicities of the eigenvalues sometimes could also be estimated in advance, for example, according to the problem symmetry.

The eigenpairs of the matrix pencil $H2−λH1$ have been proved to be related to the ones of the original NEP (Asakura et al., 2009). To compute them, the singular value decomposition (SVD) is applied on $H1$ to obtain

where $U,V∈ℂKL×KL$ are unitary matrices, $VH$ is the conjugate transpose of matrix V, and $Σ$ is a diagonal matrix with its non-negative entities $σ1,σ2,…,σKL$ in decreasing order. To exclude irrelevant results, the SVD can be truncated according to a gap scheme (Zheng et al., 2018). Consequently, the original NEP is transformed into a normal linear one of looking for the eigenpairs $(λj,ηj)$ of the matrix

Common eigensolvers such as the lapack routines can be used to compute the eigenpairs of $H3$⁠. After obtaining them, the eigenvectors $ψj$ of the original problem can be recovered by

where S is a matrix in the form of $[S0, S1, ⋯, SK−1]$⁠, and

In numerical analysis, the contour integrals contained in Eqs. (7) and (12) can be computed approximately by means of numerical quadrature like the trapezoidal rule.

It is known that the BEM suffers from the spurious eigenfrequency trouble in solving exterior acoustic problems (Burton and Miller, 1971; Zheng et al., 2015). Similarly, spurious eigenfrequencies also turn up in the BEM solution for interior acoustic problems (Zheng et al., 2019). However, such spurious eigenfrequencies are often complex since they are consistent with the eigenfrequencies of the corresponding exterior problems. Owing to the damping, the true eigenfrequencies of cavities with absorbing boundary conditions are also complex numbers for which imaginary parts have the same sign as the spurious ones, which leads to the difficulty in identifying the spurious eigenfrequencies directly. Therefore, a scheme based on the combined boundary integral formulation (Burton and Miller, 1971) is then used to identify the spurious eigenfrequencies. To obtain the combined formulation, we first take the normal derivative of Eq. (3) to acquire the normal derivative BIE,

The combined formulation for the problem with IBC can be written as

where α is a complex coupling parameter and its imaginary part should be nonzero. The real effect of the combined formulation has been found in the previous research (Zheng et al., 2015) to remove the spurious eigenfrequencies in the complex domain but not really free of them. Unlike the situation in the acoustic radiation analysis in which a fixed coupling parameter is usually utilized, it will be indicated in Sec. 4 that a variant-parameter scheme can be applied to identify and filter out the spurious eigenfrequencies in the BEM eigenanalysis.

It should be further mentioned that strongly- and hyper-singular boundary integrals exist in Eq. (14). In this letter, discontinuous quadratic BEs are used to achieve good accuracy of numerical results. To evaluate the singular integrals, two identities given by Liu and Rudolphi (1999) are employed to transform the strongly- and hyper-singular integrals into weakly- and strongly-singular ones, respectively. Cauchy's principal value of the strongly-singular integral is then calculated according to the method given by Paget (1981).

A circular model with IBC is used to verify the accuracy and effectiveness of the developed scheme. The model having radius $a=1 m$ is divided into 360 discontinuous quadratic BEs. The acoustic medium is air with mass density $ϱm=1.2 kg/m3$ and sound speed $C=340 m/s$⁠. From the general solution to the problem, the analytical eigenfrequencies are found to be the roots of

where J_n is the nth order Bessel function of the first kind, $Z0=Z/(ϱmC)$⁠, and Z is the acoustic impedance determined by the empirical Delany-Bazley model (Delany and Bazley, 1970), which represents a porous material by the following complex propagation constants:

where $X=ϱmf/Rf$⁠, and R_f denotes the flow resistivity of the porous material. The impedance from the porous layer backed by a sound-hard wall is calculated by $Z=iZ̃cot(k̃h)$⁠, where h is the thickness of the layer.

In the following calculation, $Rf=104 N s/m4, h=0.1 m$⁠. The contour path is specified as a circle with radius $ρ=2.2$ and center $γ=(2.9,−1)$⁠. The integrals contained in Eqs. (7) and (12) are computed by the 256-point trapezoidal rule. Numerical eigenfrequencies are given in wave number in Table 1, and some mode shapes are depicted in Fig. 1. In Table 1, the numbers underlined are spurious eigenfrequencies which are known to be consistent with the eigenfrequencies of the corresponding exterior problem and happen to be the roots of $Hn(1)(ka)=0$ in this example. Thus, the relative errors ϵ_real and ϵ_imag defined by

are given in Table 1 for both true and spurious eigenfrequencies. In Eqs. (18) and (19), k_ana and k_num denote the analytical and numerical frequencies, $ℜ(·)$ and $ℑ(·)$ indicate the real part and imaginary part of a complex number, respectively.

The FEM solver in the software COMSOL Multiphysics is further used to solve the problem. The acoustic domain is divided into 33 746 quadratic triangular elements (67 853 DOFs). The mean size of the FEM mesh is close to that of the BEM mesh and the number of elements is also 360 on the boundary. As IBC is a nonlinear function of frequency, the FEM discretization leads to a quadratic eigenproblem. COMSOL reformulates it as a linear one and solves it iteratively. It is found in Table 1 that the true eigenfrequencies calculated by the present method are much more accurate than the ones calculated by COMSOL, especially their imaginary parts. It is known that the imaginary parts of the eigenfrequencies represent damping caused by absorbing materials and are very important in engineering applications like the choice of porous materials in the cavity noise control of tires (Zhang et al., 2019). It should be mentioned that we change the signs in the imaginary parts of the FEM results since the time dependence $eiωt$ is used in COMSOL.

Moreover, it is found in Table 1 that both true and spurious eigenfrequencies have negative imaginary parts, which brings trouble on identifying the spurious ones. The combined formulation is then used in the BEM, and the numerical results are depicted in Fig. 2, where the coupling parameter α in the combined formulation is set to $α=β·i/k$⁠, with $β=−0.4,−0.3,−0.2,−0.1,0,0.1,0.2,0.3,0.4$⁠, respectively. It can be found in Fig. 2 that the true eigenfrequencies stay stationary, but the spurious ones move obviously when the coupling parameter changes. Therefore, a variant-parameter scheme can be suggested to filter out the spurious eigenfrequencies. In this scheme, the combined formulation with different coupling parameters is first used in the BEM eigenanalysis. For instance, when two coupling parameters $α1=β1·i/k$ and $α2=β2·i/k$ are used, two set of results can be obtained. These two sets of results are then compared and the frequencies which satisfy $|Δk|>κ|β2−β1|$ are dropped. Here, $|Δk|$ is the absolute difference of two relevant frequencies, and κ is a constant which is suggested to be 0.5 in this letter.

A numerical framework based on the BEM is presented in this letter for the analysis of acoustic resonances in cavities with IBCs. The NEP arising from the BEM is transformed into a standard linear one through a contour integral approach. A variant-parameter scheme based on the combined boundary integral formulation is given to filter out complex spurious eigenfrequencies which are similar to true eigenfrequencies. The high accuracy and effectiveness of the method is verified through an instructive example. Further studies are recommended to examine the engineering applications of the developed method.

This work is supported by the National Natural Science Foundation of China (Grant Nos. 11872168 and 11674082).