AN EFFICIENT NUMERICAL METHOD FOR LARGE-SCALE MODAL ANALYSIS USING BOUNDARY ELEMENT METHOD
-
摘要: 随着大规模快速边界元计算技术的发展，在复杂结构的动态设计、振动与噪声分析中愈来愈多地采用边界元法，因此求解大规模边界元特征值问题、进行复杂结构和声场模态分析，成为工程应用中一个十分重要，但却极具挑战性的课题，目前国际上还没有十分有效的数值方法.本文针对边界元法中典型的非线性特征值问题，提出了一种通用、高效的数值解法，称为基于预解矩阵采样的Rayleigh-Ritz投影法，记为RSRR.首先，通过求解一系列频域边界元问题来构造特征向量搜索空间，进而可以采用Rayleigh-Ritz投影，将原问题转化为一个可以采用现有方法求解的小规模缩减特征值问题；其次，为了降低Rayleigh-Ritz投影过程的计算量，基于解析函数的Cauchy积分公式，构造了边界元系数矩阵的插值近似方法，以及缩减特征值问题系数矩阵的快速计算方法，给出了插值项数的估计策略；最后，将RSRR与声学快速边界元法结合，应用于大规模吸声结构的复模态分析.数值算例表明，RSRR方法能够可靠地求出给定频段内的全部特征值和特征向量，具有计算效率高、精度高、通用等优点.Abstract: Thanks to the great advances in fast boundary element method (BEM) achieved in the last two decades, the BEM has been increasingly used in the dynamic design of engineering structures, the analysis of noise and vibration. Consequently, solving large-scale eigenvalue problems and performing modal analysis for complicated structures and acoustic fields using the BEM becomes an very important but challenging task; so far there are no effective numerical methods for this purpose. This paper aims to extend the application of the newly-developed resolvent sampling based Rayleigh-Ritz projection method (RSRR) to the solution of the general nonlinear eigenvalue problems (NEP) in BEM. First, in order to generate reliable eigenvector search spaces, a series of BEM linear systerms in frequency domain are solved. Then the original NEP can be transformed to a reduced NEP based the classical Rayleigh-Ritz procedure, and the reduced NEP could be solved by those exiting NEP solvers easily. Second, to reduce the prohibitive computational burden involved in solving the projected NEP by the Rayleigh-Ritz procedure, a BEM matrix interpolation technique and a fast computation method for reduced NEP systerm matrix are proposed based on the discretized Cauchy integral formula of analytic functions. Then a simple rule for estimating the number of terms in the interpolation is proposed as well. Finally, the RSRR method is used to solve large-scale practical acoustic modal analysis problems using fast BEM with complicated sound absorbing boundary conditions. Numerical results indicate that the method can robustly dig out all the interested eigenvalues and the corresponding eigenvectors with good accuracy and high computational efficiency.
-
引言
边界元法[1](boundary element method，BEM)是一种重要的工程数值方法，在计算力学、声学等领域应用很广[2-3].尤其是对于中高频振动和波的传播问题[4]，边界元法不仅能有效降低求解规模，还能避免数值色散误差，具有很大的发展潜力.近年来，边界元快速算法逐步发展成熟，显著提高了边界元大规模分析的计算效率，将边界元法向实际工程应用推进了一大步[5].
边界元法用于结构减振降噪分析，常常需要进行模态分析，以获得系统的固有频率和模态.边界元模态分析的难点在于求解如下形式的非线性特征值问题
$$ {\pmb T}(\lambda )v = {\bf 0} $$ (1) 式中，$\lambda $和$v$分别为特征值和与之对应的特征向量，系统矩阵${\pmb T} (z) \in \mathbb{C}^{n\times n}$是声波数$z = \omega/ c$的函数，$\omega $是声波的圆频率，$c$是声速，$n$表示系统总自由度数.特别地，当${\pmb T} (z)$为$z$的线性或二次函数时，上述问题即为标准特征值问题或二次特征值问题；后者是最简单的一类非线性特征值问题，普遍存在于有限元结构分析领域，其数值解法已经较为成熟[6].
在边界元法中，矩阵${\pmb T} (z)$关于$z$的表达式十分复杂，甚至无法得到.对这类非线性特征值问题，数学上还没有一个通用、有效、适合大规模求解的数值解法[7-9].加之边界元矩阵${\pmb T} (z)$还是一个非对称的满阵，计算量和存储量都很大，过去采用的行列式搜索法[10]、域内离散法[11]、域积分转换法[12-14]等都无法用于大规模求解.
围道积分法(contour integral method，CIM)是近年来才提出的一种求解非线性特征值问题的数值解法[15-17]，最初被用于求解广义特征值问题[18].它不需要迭代，能一次计算出复平面给定区域内的全部特征值和相应特征向量；同时还适合于并行求解，计算效率很高，具有广阔的应用前景.目前该方法已被用于简单边界元特征值问题的求解[19-20].
但是，CIM可能导致漏根或是虚假特征根，求解精度难以保证.这是因为它采用矩阵${\pmb T}(z)$的高阶矩来构造特征向量空间[19-21].矩的阶数越高，这些矩阵越趋于线性相关，从而导致特征空间病态或构造失败.另外，按照CIM的思想，频率采样点只能落在包含寻根区域的某个围道上，而实际计算结果表明，采用域内更为靠近特征值的采样点，求解精度会更高[22].
为解决上述问题，最近我们发展了一种稳定、高效的非线性特征值问题数值新解法，称为基于预解矩阵采样的Rayleigh-Ritz投影法，简记为RSRR (resolvent sampling based rayleigh-ritz projection method)[22].它不仅可以构造出可靠的特征空间、自由布置频率采样点以提高求解精度，而且具备CIM能一次求出指定范围内所有特征值、适合并行计算的优点.目前，RSRR已被成功用于有限元和一类简单的边界元大规模非线性特征值问题的求解，在精度、计算效率和稳定性等方面均显示出较大优势[22-23].
本文把RSRR应用于工程边界元减振降噪分析中，解决Rayleigh-Ritz投影法中边界元矩阵${\pmb T}(z)$的快速计算问题.该问题是制约大规模边界元特征值问题求解的瓶颈[24]，因为求解Rayleigh-Ritz投影所得的缩减特征值问题，往往需要成百上千步的迭代或矩阵求逆，而每一步都需要计算一个稠密矩阵${\pmb T}(z)$.对于自由度几十万、上百万的大规模问题，其计算量是普通计算设备无法承受的.
本文提出一种基于Cauchy积分公式的边界元矩阵插值方法，可以在复平面内任意区域上快速构造Rayleigh-Ritz缩减矩阵的高精度插值表达式，从而避免了在求解缩减特征值问题中计算原矩阵${\pmb T}(z)$的难题, 大大降低了RSRR求解边界元特征值问题的计算量；还给出了一种简便的插值项数确定方法，可根据精度要求确定一个近似最少的插值项数，有效控制整个求解过程的计算量.与文献[22]中的一维Chebyshev插值法相比，本文的Cauchy插值法可在复平面内任意二维区域上构造插值近似，可用于解决高阻尼、高吸声减振降噪设计中常见的复特征值问题[25-28].
1. 边界元特征值问题
本文考虑吸声结构分析中的边界元特征值问题，但所提之方法具有通用性，也适用于弹性动力学、电磁学等问题.
声学问题的边界积分方程为
$$ c(x)p(x) + \int _{\partial D} {\dfrac{\partial F(x, y)}{\partial n_y }p(y)\text{d} y} =\\ \qquad \int _{\partial D} {F(x, y)\dfrac{\partial p(y)}{\partial n_y }\text{d} y} $$ (2) 式中，$D$为声场区域，$\partial D$为$D$的边界，$p$为声压，$F(x, y) = {\exp ({\rm i}z\vert x -y\vert)} / {\vert x -y\vert }$是三维Helmholtz方程的基本解，$c(x)$为自由项，当$\partial D$在$x$附近光滑时，$c(x) = 1 / 2$，$n_y $表示$y$点处的外法向.
采用文献[29]的Nyström边界元法离散式(2)，选择核无关快速定向压缩算法[30]对其进行加速.将声场边界$\partial D$离散为二次曲面三角形单元，用6节点二次协调元进行几何近似，用6节点二次非协调元进行函数近似.采用非协调元可以避免边界元角点问题[4]，方便几何建模，有益于方法的通用性.同时，我们的计算经验表明，相同自由度数目情况下，采用非协调元可以获得与协调元相同的精度.经Nyström方法离散边界积分方程(2)，最终可得如下标准形式的边界元离散方程
$$ {\pmb H}(z){\pmb p}(z) ={\pmb G}(z){\pmb q}(z) $$ (3) 式中，${\pmb H}$和${\pmb G}$是边界元离散矩阵，${\pmb p}$和${\pmb q}$分别是包含节点声压及其法向导数的列向量.矩阵${\pmb H}$和${\pmb G}$中的近场修正部分仍需要进行常规边界元法中的奇异和近奇异积分计算.本文奇异积分采用文献[31]中的数值计算方法，而近奇异积分采用逐步细分的方法来计算.
考虑Robin边界条件
$$ \alpha (x)p(x) + \beta (x)\dfrac{\partial p(x)}{\partial n_x } = \gamma (x) $$ (4) 当$\alpha \ne 0, \beta = 0$时，式(4) 退化成Dirichlet边界条件，当$\alpha = 0$, $\beta \ne 0$时，式(4) 退化成Neumann边界条件.
下面考虑吸声结构分析中常见的一类Robin边界条件，即阻抗边界条件
$$ \dfrac{\partial p(x)}{\partial n_x } = - {\rm i}\dfrac{\rho _0 \omega }{Z(x, z)}p(x) $$ (5) 式中，$Z$为边界声阻抗，与节点位置$x$和频率$z$有关；$\rho _0 $为声场介质密度，${\rm i} = \sqrt { -1} $.式(5) 在每个边界元节点$x_j $上都成立，将这些式子统一写成矩阵形式即为
$$ {\pmb q}(z) = {\pmb B}(z) {\pmb p}(z) $$ (6) 式中，${\pmb B} (z)$是维度为$n\times n$的对角矩阵，其对角线元素表示为$ -{\rm i}\rho _0 \omega / Z(x_j, z)$.将式(6) 代入式(3)，可得
$$ [{\pmb H}(z)-{\pmb G}(z) {\pmb B} (z)] {\pmb p}(z) = {\bf 0} $$ (7) 上式即为吸声结构模态分析中所要求解的特征方程，系统矩阵${\pmb T}(z) ={\pmb H}(z) -{\pmb G}(z) {\pmb B}(z)$. Dirichlet和Neumann边界条件所对应的特征方程可类似得到.
边界元法中，矩阵${\pmb T}(z)$的元素是由基本解或其法向导数经单元积分后得到的.由于基本解是频率$z$的复杂函数，加之单元积分通常只能靠数值方法计算，而无法获得解析表达式，因此${\pmb T}(z)$关于$z$的解析表达式通常无法得到，而只能对于给定的$z$通过数值方法来计算.这不仅增大了边界元特征值问题求解的难度，而且还导致其计算量很大，以至于目前仍没有有效的大规模数值解法，本文正是要解决这些问题.
2. 非线性特征值问题的RSRR解法
RSRR算法是基于Rayleigh-Ritz投影的特征值解法.给定一个包含待求解特征向量的搜索空间${\pmb S}$，设${\pmb S}$的一组正交基为${\pmb Q}$，则原非线性特征值问题(1) 可以通过Rayleigh-Ritz投影转化为如下的缩减形式
$$ {\pmb T}_Q (z) {\pmb g} = {\bf 0} $$ (8) 式中，${\pmb T}_Q (z) ={\pmb Q}^{\rm H}{\pmb T}(z) {\pmb Q} \in \mathbb{C}^{k_S \times k_S }$，$k_S $为${\pmb Q}$的列向量数目, 上标H表示共轭转置.式(8) 与式(1) 具有相同的特征值，设$(\lambda, {\pmb g})$为${\pmb T}_Q (z)$的任一特征对，则对应原问题的特征对为$(\lambda, {\pmb Q \pmb g})$.下面首先介绍搜索空间${\pmb S}$的构造方法，然后给出RSRR算法的基本步骤.
2.1 特征空间的构造
本文的预解矩阵采样法(resolvent sampling method，RES)是构造特征空间${\pmb S}$的一种稳定、高效的数值方法.设要计算复平面上由简单闭合曲线$C$围成的区域内的所有特征值和相应特征向量，用RES构造特征空间可分为以下两步：
(1) 在围道$C$(见图 1)内选择$N$个采样点$z_k $，$k = 1, 2, \cdots, N$，同时生成$L$个线性无关的$n$维随机向量，组成矩阵${\pmb U} \in \mathbb{\pmb C}^{n\times L}$；
![]()
图 1 围道$C$示意图.图中闭合轮廓线代表围道$C$，"+"代表$C$内的待求特征值，"$\circ$"代表$C$外的特征值，"$\bullet$"代表采样点Figure 1. A diagram showing the contour $C$ (closed contour). eigenvalues interested (+), eigenvalues not interested ($\circ$) and sampling points in the contour ($\bullet$)(2) 对$k = 1, 2, \cdots, N$，计算${\pmb T}(z_k)^{ -1}{\pmb U}$，并将其组装成矩阵${\pmb S} \in \mathbb{C}^{n\times NL}$
$$ {\pmb S} = \left[{\pmb T}(z_1 )^{-1}{\pmb U} \, , \ \ {\pmb T}(z_2 )^{-1}{\pmb U}\, , \ \ \cdots, \ \ {\pmb T}(z_N )^{-1}{\pmb U} \right] $$ (9) 矩阵${\pmb S}$的列空间$span(S)$即为RES所形成的特征向量搜索空间.
下面首先证明：当采样点布置、参数$L$和$N$选取都合理的情况下，空间$span({\pmb S})$将近似包含围道$C$内的所有特征向量，在此基础上讨论采样点和随机矩阵${\pmb U}$对RES方法效率的影响情况.
边界元矩阵${\pmb T}(z)$虽然形式复杂，但由于基本解是频率$z$的全纯函数(holomorphic function)，因此${\pmb T} (z)$也是$z$的全纯矩阵值函数(holomorphic matrix-value function).根据解析矩阵函数的Keldysh定理[16]，边界元矩阵${\pmb T} (z)$的逆矩阵${\pmb T}^{-1}(z)$可由$C$内的特征值和特征向量表示为如下形式
$$ {\pmb T}(z)^{-1}=\sum\limits^{n_C}_{k=1}\sum\limits^{\eta_k}_{l=1}\sum\limits^{\mu_{l, k}}_{j=1} (z-\lambda_k)^{-j}\cdot \\ \qquad \sum\limits^{\mu_{l, k}-j}_{m=0} {\pmb v}^{l, k}_m ({\pmb w}^{l, k}_{\mu_{l, k-j-m}})^{\rm H}+{\pmb R}_C(z) $$ (10) 式中，$n_C $为不计重根的互不相同特征值的数目，$\eta _k $为$T(\lambda _k)$的几何重数，$\mu _{l, k} $表示$\lambda _k $的局部重数，${\pmb v}_m^{l, k} $和${\pmb w}_{\mu _{l, k -j -m} }^{l, k} $分别是${\pmb T} (\lambda _k)$的右、左正规化广义特征向量，${\pmb R}_C (z)$表征$C$之外的所有特征值的贡献，详细理论可参考文献[16, 32].
将式(10) 的右端项写成矩阵形式
$$ {\pmb T}(z)^{ - 1} ={\pmb V}_C {\pmb\varPhi} (z){\pmb W}_C^{\rm H} + {\pmb R}_C (z) $$ (11) 其中，${\pmb\varPhi} (z)$是由${\pmb T}(z)$特征值的Jordon块组成的矩阵，${\pmb V}_C $和${\pmb W}_C $分别是矩阵${\pmb T}(z)$和${\pmb T}(z)^{\rm H}$在特征根$\lambda $处的广义特征向量标准系统(canonical system of generalized eigenvectors).为了从${\pmb T}(z)^{ -1}$中提取出特征向量空间$span({\pmb V}_C)$，必须尽可能地消除${\pmb R}_C (z)$的影响.由于${\pmb R}_C (z)$在特征值附近解析，因此可以用一组合适的基函数$\{g_j (z)\}$来逼近它
$$ {\pmb R}_C (z) \approx \sum\limits_{j = 1}^J {\pmb R}_j g_j (z) $$ (12) 式中，${\pmb R}_j $为展开系数矩阵，$J$为展开项数.当$\{g_j (z)\}$选择合适时，式(12) 的逼近过程可以很快收敛.例如，对于区间上的解析函数，用Chebyshev多项式逼近可以获得指数收敛速度.
现对${\pmb T} (z)^{-1}$右乘以包含$L$个线性无关向量的矩阵${\pmb U}$，并将式(12) 代入式(11)，整理可得
$$ {\pmb T} (z)^{ - 1}{\pmb U} = {\pmb V}_C \bar {\pmb \varPhi }(z) + \sum\limits_{j = 1}^J \bar{\pmb R}_j g_j (z) = \hat{\pmb S}{\pmb Z}(z) $$ (13) 式中，$\bar {\pmb\varPhi }(z) ={\pmb \varPhi} (z){\pmb W}_C^{\rm H} {\pmb U} \in \mathbb{C}^{\bar {n}_C \times L}$，$\bar {n}_C $表示待求特征值的数目，$\bar {\pmb R}_j ={\pmb R}_j {\pmb U} \in \mathbb{C}^{n\times L}$，$\hat {\pmb S}$和${\pmb Z}(z)$的表达式如下
$$ \hat {\pmb S} = \left[ {{\pmb V}_C } \ \ {\bar{\pmb R}_1 } \ \ {\bar {\pmb R}_2 } \ \ \cdots \ \ {\bar {\pmb R}_J } \right] \in \mathbb{C}^{n\times \left( {\bar {n}_C + J \cdot L} \right)} $$ (14) $$ {\pmb Z}(z) = \left[ {\bar {\pmb \varPhi }(z)} \ \ {{\pmb D}_1 (z)} \ \ {{\pmb D}_2 (z)} \ \cdots \ {{\pmb D}_J (z)} \right]^{\rm T} \in \mathbb{C}^{\left( {\bar {n}_C + J \cdot L} \right)\times L} $$ (15) 式中${\pmb D}_j (z) = {\rm diag} (g_j (z), g_j (z), \cdots, g_j (z)) \in \mathbb{C}^{L\times L}$.
由式(14) 可知，${\pmb V}_C $中的特征向量包含在$\hat{\pmb S}$的列空间中，即
$$ span({\pmb V}_C ) \subset span(\hat{\pmb S}) $$ 因此，$span(\hat {\pmb S})$是有效的特征空间，$\hat {\pmb S}$和式(9) 定义的${\pmb S}$有如下关系
$$ {\pmb S} \approx \hat {\pmb S} \cdot \left[{{\pmb Z}(z_1 ), \; {\pmb Z}(z_2 ), \; \cdots, \; {\pmb Z}(z_N )} \right] = \hat{\pmb S}\hat{\pmb Z} $$ (16) 式中，$\hat{\pmb Z} = \left[{{\pmb Z} (z_1), \; \; {\pmb Z} (z_2), \; \; \cdots, \; \; {\pmb Z} (z_N)} \right] \in \mathbb{C}^{\left({\bar {n}_C + J \cdot L} \right)\times N \cdot L}$.由式(16) 可知，当矩阵$\hat{\pmb Z}$的秩等于$\bar {n}_C + J \cdot L$时有
$$ span(\hat{\pmb S}) = span({\pmb S}) $$ (17) 此时，$N$和$L$应满足$\bar {n}_C + J \cdot L \leqslant N \cdot L$, 即
$$ (N - J)L \geqslant \bar {n}_C $$ (18) 以上推导表明，当采样点的布置、参数$L$、$N$和$J$取值合理的情况下，式(9) 定义的采样矩阵${\pmb S}$所张成的子空间会包含所有待求的特征向量.同时也表明，这几个参数对RES方法的效率有重要影响.但在实际应用中，这些参数的理想值事先很难确定.下面通过分析，给出这些参数的选取原则.
由上面推导过程可见，采样点$z_k $的作用是形成对未知矩阵${\pmb R}_C (z)$的插值，从而隐式地确定式(12) 中的基函数$\{g_j (z)\}$.据此，对于实特征值问题，可将采样点取为寻根区间上的Chebyshev点，这对应于Chebyshev插值，精度较高；而对复特征值问题，围道$C$通常为包含寻根区域的矩形，此时采样点可取为该矩形区域上的二维Chebyshev点.
采样点数目$N$和随机向量数目$L$共同决定着RES方法的精度和计算效率.就保证特征空间精度而言，大量计算经验表明，为满足式(18) 一般需要$N \cdot L > 2\bar {n}_C \sim 3\bar {n}_C$.同时，$U$中的$L$个随机向量的作用是探测多重特征值所对应的特征向量，因此$L$的取值应该不小于$T$在$C$内各特征值的最大局部重数[16]，即
$$ L \geqslant \mathop {\max }\limits_{k = 1, \cdots, n_C } \left( {\sum\limits_{l = 1}^{\eta _k } {\mu _{l, k} } } \right) $$ (19) 文献[21]研究表明，相比于$L$，增加$N$对提高精度的效果更显著.但是另一方面，由式(9) 可知，计算矩阵$S$需要求解$N$个线性方程组，每个包含$L$个右端向量，因此增大$N$和$L$都会大幅增加计算量.综上，一般在满足式(18) 的条件下，$N$和$L$应尽量小；而取$N = (2\bar {n}_C \sim 3\bar {n}_C) / L$.特征值数目$\bar {n}_C $并不要求准确地知道，因此可利用粗网格下的小规模问题进行预估，也可通过理论分析或参考相似结构的已知结果而获得.
2.2 RSRR算法
RSRR算法求解大规模边界元非线性特征值问题的步骤如下：
(1) 确定围道$C$、采样点数$N$、采样点$z_k, \; (k = 1, \; 2, \; \cdots \; , \; N)$、重数$L$、随机矩阵${\pmb U} \in \mathbb{C}^{n\times L}$、展开项数$d$；
(2) 计算${\pmb T}(z_k)^{ -1}{\pmb U}$, $k = 1, \; 2, \; \cdots, \; N$；
(3) 组装${\pmb S}$矩阵，并对其进行截断奇异值分解${\pmb S} \approx {\pmb Q \pmb\varSigma \pmb W}^{\rm H}$，得到一组正交基${\pmb Q}$；
(4) 计算缩减矩阵${\pmb T}_Q (z) = {\pmb Q}^{\rm H}{\pmb T} (z) {\pmb Q}$的插值近似，计算方法将在第3节详述；
(5) 求解缩减后的非线性特征值问题${\pmb T}_Q (z) {\pmb g} ={\bf 0}$，设所得的特征对为$(\lambda _j, {\pmb g}_j)$，$j = 1, 2, \cdots, \bar {n}_C $，则原问题的特征对为$(\lambda _j, {\pmb v}_j)$, $j = 1, 2, \cdots, \bar {n}_C $，${\pmb v}_j ={\pmb Q}{\pmb g}_j $.
在第三步的Rayleigh-Ritz投影中，计算正交基${\pmb Q}$的截断奇异值分解的奇异值舍去阈值一般设为$\delta = 10^{ -14}$，而${\pmb Q}$的列数$k_s $即为保留的奇异值数目.为了保证求出积分路径内的所有特征根，$k_s $应满足
$$ k_s \geqslant rank({\pmb V}_C ) \geqslant \bar {n}_C $$ (20) 奇异值分解的计算复杂度为$O(n)$，且一般$k_s $仅为几十到几百，因此奇异值分解的计算量在整个边界元特征值分析中可忽略不计.
实际应用中，会出现某些采样点$z_k $靠近特征值的情况，此时${\pmb T}(z_k)^{ -1}{\pmb U}$的数值远大于${\pmb S}$中其他向量，从而导致特征空间的总体精度下降，因此在计算正交基${\pmb Q}$之前需要对${\pmb S}$各列进行正规化.
3. 基于Cauchy积分的矩阵插值方法
经过Rayleigh-Ritz投影后的非线性特征值问题规模很小，见式(8)，现有的许多方法都可以对其进行求解[15-16, 22].本文采用文献[22]改进的Block Sakurai-Sugiura (Block SS)算法，该算法能够自动确定待求特征值的数目，并且计算出给定区域内的全部特征值及对应特征向量.但是，无论是Block SS算法还是其他方法，求解式(8) 都需要计算多个给定频率$z_i $处的系数矩阵${\pmb T}_Q (z_i) ={\pmb Q}^{\rm H}{\pmb T}(z_i){\pmb Q}$.由于边界元矩阵${\pmb T}(z_i)$是大规模稠密矩阵，这部分计算量将十分巨大，是制约边界元大规模非线性特征值问题求解的一个瓶颈.
3.1 矩阵近似
为有效降低多次计算${\pmb T}_Q (z_i)$的计算量，文献[22]通过Chebyshev插值方法获得了${\pmb T}_Q (z)$的Chebyshev插值近似格式，在求解${\pmb T}_Q (z) {\pmb g} = 0$时只需要调用小规模Chebyshev插值式计算${\pmb T}_Q (z)$，从而避免了反复计算边界元矩阵${\pmb T} (z)$.但是，由Chebyshev插值理论可知，此方法只有在特征值分布在直线附近时(如实特征值问题)才能保证足够的精度.而一般阻尼吸声结构的特征值分布在二维复平面内，该方法不再适用.
为克服Chebyshev插值方法的局限性，本文基于Cauchy积分法构造${\pmb T}_Q (z)$在复平面内任意单连通区域上的插值近似格式.由于边界元矩阵${\pmb T} (z)$为全纯函数，因此可由Cauchy积分公式表示为
$$ {\pmb T} (z) = \dfrac{1}{2{\rm{\pi }} {\rm i}}\oint_C \dfrac{{\pmb T} (\zeta )}{\zeta - z}\text{d}\zeta $$ (21) 式中，$C$为复平面内任意积分路径，本文中它与第2节的RSRR算法中的围道取为一致.采用数值积分计算上式的围道积分，可得
$$ {\pmb T}(z) = \dfrac{1}{2{\rm{\pi }} {\rm i}}\sum\limits_{i = 1}^d \dfrac{{\pmb T} (z_i )w_i }{z_i - z} $$ (22) 这里沿用文献[30]中对标量函数的叫法，将式(22) 称为矩阵函数${\pmb T} (z)$的第一类Cauchy插值公式.其中，总项数$d$决定了数值积分的精度，$w_i $为积分点$z_i $对应的权系数，它们由所采用的数值积分公式确定.对于光滑的闭合曲线，通常采用梯形公式可获得指数阶的计算精度.
由于是在复平面内进行矩阵插值，Cauchy插值公式比Chebyshev插值公式适用范围更广.但是，当$z$逼近$z_i $时，$1/{(z_i -z)}$趋于无穷大，式(22) 的精度会迅速下降，导致插值失败.针对这一现象，用${({\pmb T} (\zeta) -{\pmb T} (z))}/{(\zeta -z)}$代替式(21) 中的${{\pmb T} (\zeta)} / {\left({\zeta -z} \right)}$，前者在$z$逼近$\zeta $时分子分母同时趋于0，不存在奇异性，且其Cauchy积分值为0，即
$$ \dfrac{1}{2{\rm{\pi }} {\rm i}}\oint_C {\dfrac{{\pmb T} (\zeta ) - {\pmb T} (z)}{\zeta - z}\text{d} \zeta } = {\bf 0} $$ (23) 上式用数值积分公式离散后可得
$$ \dfrac{1}{2{\rm{\pi }} {\rm i}}\sum\limits_{i = 1}^d {\dfrac{\left[{{\pmb T} (z_i )-{\pmb T} (z)} \right]w_i }{z_i - z}} \approx {\bf 0} $$ (24) 整理上式，得到${\pmb T} (z)$的表达式如下
$$ {\pmb T}(z) \approx \sum\limits_{i = 1}^d \dfrac{{\pmb T} (z_i )w_i }{z_i - z} \Big/ \sum\limits_{i = 1}^d \dfrac{w_i }{z_i - z} =\\ \qquad \sum\limits_{i = 1}^d {\pmb T} (z_i )\varphi _i (z) $$ (25) 称上式为第二类Cauchy插值公式，式中
$$ \varphi _i \left( z \right) \approx {\dfrac{w_i }{z_i - z}} \Big / {\sum\limits_{i = 1}^d {\dfrac{w_i }{z_i - z}} } $$ (26) 将式(25) 代入${\pmb T}_Q (z_i) = {\pmb Q}^{\rm H}{\pmb T} (z_i) {\pmb Q}$即可得到缩减特征值问题系数矩阵的近似式
$$ {\pmb T}_Q (z) \approx \sum\limits_{i = 1}^d {\underbrace {{\pmb Q}^{\rm H}{\pmb T} (z_i ) {\pmb Q}}_{{\boldsymbol{C}}_i }\varphi _i (z)} = \sum\limits_{i = 1}^d {{\pmb C}_i \varphi _i (z)} $$ (27) 其中，系数矩阵${\pmb C}_i = {\pmb Q}^{\rm H}{\pmb T}(z_i) {\pmb Q}$为$k_S $维方阵，规模很小，可以直接计算并存储起来. ${\pmb C}_i $的主要计算量是在$d$个插值点$z_i $上计算边界元矩阵${\pmb T} (z_i)$，并通过矩阵向量乘法计算${\pmb T} (z_i) {\pmb Q}$.为有效控制计算量，需使$d$尽量小，第3.2节将提出$d$的估计方法.同时，由于$d$个插值点是预先给定的，因此$d$个系数矩阵${\pmb C}_i $很适合在不同处理器上并行计算.
相较第一类Cauchy插值公式(22)，第二类Cauchy插值公式(25) 不仅消除了当$z$逼近$z_i $时因奇异性导致的插值不稳定现象，还能显著提高插值精度.文献[33]称这一现象为"过收敛"现象，即式(25) 不仅在积分区域$C$内满足精度要求，在边界$C$附近也依然具有较高的精度.
下面通过一个边界元矩阵插值算例来展示"过收敛"现象，证明第二类Cauchy插值公式可以为本文的RSRR算法提供更充足的矩阵近似精度.选择算例4.2的边界元矩阵和频率区域，展开项数$d$取20，矩阵插值的相对误差按下式计算
$$ \varsigma = \log \dfrac{\left| {{\pmb T}_Q^\ast (z){\pmb y} - {\pmb T}_Q (z){\pmb y}} \right|}{\left| {{\pmb T}_Q^\ast (z){\pmb y}} \right|} $$ (28) 其中${\pmb T}_Q^\ast (z)$表示直接计算出的缩减边界元矩阵，${\pmb y}$为随机向量.第一类Cauchy插值和第二类Cauchy插值的相对误差云图见图 2，图中"$\circ$"代表积分路径上的20个插值点.可见，第一类Cauchy插值公式(22) 仅在区域中心能够保持一定精度，靠近积分边界时，其插值精度迅速恶化.而第二类Cauchy公式(25) 不仅在积分区域内能够保持高精度，在积分区域边界及其外部一定范围内也仍然呈现出很高的精度，明显优于第一类Cauchy公式，能够满足RSRR方法对矩阵近似精度的要求.
3.2 确定展开项数
在保证矩阵近似精度的前提下，展开项数$d$越小，近似矩阵的计算量越小.这里提出一种估计最小展开项数$d$的方法.根据3.1节的Cauchy矩阵插值方法，矩阵${\pmb T}(z)$的插值公式(25) 同样适用于对矩阵每个元素的插值近似，因此可以用矩阵元素插值所需的最小项数快速地估计$d$.
一般地，边界元法系数矩阵${\pmb T}(z)$的元素$t(z)$具有如下通式
$$ t(z) = f(z){\rm e}^{{\rm i}zr}\, , \ \ 0 \leqslant r \leqslant R $$ (29) 其中，$r$表示源点和场点的距离，$R$代表模型最大尺寸，在所关注频段内$f(z)$为$z$的解析函数.例如，采用配点法离散边界积分方程(2)，则${\pmb H}$矩阵的元素可写成
$$ t(z) = \int_\Delta {\dfrac{\partial }{\partial n_y }\dfrac{{\rm e}^{{\rm i}z\left| {x - y} \right|}}{\left| {x - y} \right|}w(y)} \text{d} y =\\ \qquad \int_\Delta {{\rm e}^{{\rm i}z\left| {x - y} \right|}\dfrac{{\rm i}z\left| {x - y} \right| - 1}{\left| {x - y} \right|^2}\dfrac{\partial \left| {x - y} \right|}{\partial n_y }w(y)} \text{d} y $$ (30) 式中，$\Delta $代表边界单元，$w(y)$代表边界元形函数，函数${\rm e}^{{\rm i}z\left| {x -y} \right|}$可展开成
$$ {\rm e}^{{\rm i}z\left| {x - y} \right|} \approx \sum\limits_k {{\rm e}^{{\rm i}zr_k }\phi _k (\vert x - y\vert )} $$ (31) 其中，$r_k=|x-y_k|$, $\phi_k$是展开基函数.故
$$ t(z) \approx \sum\limits_k {{\rm e}^{{\rm i}zr_k }\! \underbrace {\int_\Delta {\dfrac{{\rm i}z\left| {x - y} \right| - 1} {\left| {x - y} \right|^2}\dfrac{\partial \left| {x - y} \right|} {\partial n_y }\phi _k (\vert x - y\vert )w(y)} \text{d} y}_{: = f_k (z)}} =\\ \qquad \sum\limits_k {f_k (z) {\rm e}^{{\rm i}z\left| {x - y} \right|}} $$ (32) 当$r = R$时，函数$t(z)$振荡最为剧烈，对展开项数要求最为苛刻.因此，本文用函数$t(z) = f(z) {\rm e}^{{\rm i}zR}$的Cauchy插值来估计对${\pmb T}(z)$进行Cauchy插值所需的插值项数.
按式(25) 计算函数$t(z)$在围道$C$内的Cauchy插值近似$t(z) \approx \tilde {t}(z, d)$，定义近似函数$\tilde {t}(z, d)$插值的$L^2$误差为
$$ \varepsilon (d) = \sqrt {\dfrac{\int_C {\left[{t(z)-\tilde {t}(z, d)} \right]^2} \text{d} z}{\int_C {t(z)^2} \text{d} z}} $$ (33) 假设矩阵${\pmb T}(z)$的插值精度要求为$\varepsilon _0 $，则可以将满足$\varepsilon (d_0) \leqslant \varepsilon _0 $的最小整数$d_0 $作为${\pmb T}(z)$所需的插值项数的估计，即$d = d_0 $.
实际应用中，也可以通过以下3个途径获得更为可靠的$d$的估计：
(1) 选取合适的函数$f(z)$. $f(z)$通常取单项式$z^m$, $m = 0, 1, \cdots$即可，其中指数$m$的选择视基本解和边界条件而定.对于Nyström边界元法中的${\pmb G}$矩阵，$m$通常取为0，因为${\pmb G}$矩阵的近场矩阵元素就是形如$\int_\Delta {{{\rm e}^{{\rm i}z\left| {x -y} \right|}} / {\left| {x -y} \right|}w(y)} \text{d} y$的单元积分的线性组合，具有式(29) 的形式，其中的$f_k (z)$与频率无关；同理，近似${\pmb H}$矩阵时$m$取为1.
(2) 选取多个函数$t(z) = f(z) {\rm e}^{{\rm i}zr}$分别进行Cauchy插值，$d$的估计值应该对每个函数都满足插值精度要求.一般情况下，可以选择多个不同的$0 \leqslant r_k \leqslant R$来得到一系列函数$t_k (z) = f(z) {\rm e}^{{\rm i}zr_k }$.
(3) 在估计值$d_0 $的基础上给定一个安全系数$C_d $，即令$d = \left\lfloor {C_d \cdot d_0 } \right\rfloor $.由于Cauchy插值通常是指数收敛的，因此系数$C_d $通常不必太大，如$1 \leqslant C_d \leqslant 1.2$.
4. 算例
本文所有算例均在一台32核(Xeon E5-2630 v3，2.40GHz)64GB内存的工作站上进行，所有程序都采用单线程编写.边界元法快速算法的精度$\varepsilon $设定为$1.0\times 10^{ -5}$，线性方程组采用广义极小残量法(GMRES)迭代求解，收敛残差为$1.0\times 10^{ -6}$，并采用对角块预条件技术.
4.1 解析解对比算例
该算例旨在验证本文方法的正确性，并与现有方法进行计算精度和效率上的对比.
考虑单位立方体内声腔模型，将其表面离散为1200个三角形单元，边界元自由度数为7200，选择一组相互平行的表面设定为Dirichlet边界，其余4个表面设定为Neumann边界.
该问题的特征频率有解析解.文献[20]采用围道积分法计算前12个特征频率，而本文采用RSRR方法.给定采样点数$N= 20$，采样点取为实区间(150, 450) Hz上的Chebyshev点；因为该频段内特征值的最大重数为4，故取$L = 4$.矩阵插值项数$d = 20$，插值方法采用第3节构造的Cauchy插值，围道积分路径选为椭圆，其长轴与采样点所在实区间重合，短、长轴比为0.05，插值点取为积分路径上的高斯积分点.根据${\left| { {\rm Re}(\lambda _j) -\lambda _j^ * } \right|} / {\lambda _j^ * }$计算数值解$\lambda _j $与解析解$\lambda _j^\ast $的误差，见图 3.该算例共计耗时23.3 min，占用内存227.1 MB.
由图 4可以看出，RSRR方法计算精度明显高于文献[18]的围道积分法.同时本算例运用RSRR算法只需生成$N = 20$次BEM系数矩阵，求解$N\times L = 80$个BEM方程组.而文献[20]要生成256次BEM系数矩阵，求解$256\times 15 = 3\, 840$个BEM方程组，求解边界元方程的数目是RSRR的近50倍，但是求解精度却比RSRR低了近4个量级！
![]()
图 4 基于Cauchy插值的BEM矩阵和核函数最大误差曲线Figure 4. Max relative error of BEM matrix and Kernel function based on the Cauchy interpolation4.2 复杂阻抗算例
考虑一长方体内声腔模型，其长、宽、高分别为0.5 m、0.4 m和0.3 m，将其表面离散为1 256个三角形单元，边界元自由度数为7 536.在其中一个0.4 m$\times$0.3 m的表面上敷设多孔吸声材料，等效为Delany-Bazley阻抗边界条件[25]，阻抗表达式如下，其余表面为Neumann边界，关注频率区域为(300+0i, 1000+200i)Hz，围道$C$取为该矩形.
$$ \left.\begin{array}{l} Z_c (\omega ) = \rho _0 c\left[{1 + 0.057\delta ^{-0.754}-j0.087\delta ^{-0.732}} \right] \\ k(\omega ) = \dfrac{\omega }{c}\left[{1 + 0.0978\delta ^{-0.700}-j0.189\delta ^{-0.595}} \right] \\ Z = - {\rm i}\dfrac{Z_c }{\phi }\cot kh \end{array}\!\!\right\} $$ (34) 式中，$\delta = {\rho _0 \omega }/{2{\rm{\pi }} \sigma }$，$\sigma $表示流阻，$h$和$\phi $分别表示多孔吸声材料层的厚度和孔隙率，$\rho _0 $为多孔材料密度.本文取$\sigma = 2\times 10^4$N$\cdot$s/m$^4$, $h = 20$ cm, $\phi = 1$和$\rho _0 = 1.2$ kg/ m$^3$.
首先，确定最优插值项数$d$.考虑阻抗边界，给定边界元核函数的表达式$t(z) = Z \cdot z \cdot {\rm e}^{{\rm i}zr}$，分别取$r = R, \; 0.8R, \; 0.6R, \; 0.4R$，插值点取为围道上的高斯积分点，插值项数$d$依次取$10, \; 16, \; 20, \; 26, \; 30$，根据式(33) 计算不同插值项数下$t(z)$和${\pmb T}(z)$的插值误差.由于${\pmb T}(z)$为非对称满阵，直接对其进行Cauchy插值计算量太大，难以实现，此处采用了一个近似等效的方法：随机取200个边界元节点，计算这些节点在${\pmb T}(z)$中对应的$200\times 200$子矩阵块，用此子矩阵块代替${\pmb T}(z)$进行插值，取其最大插值误差代替${\pmb T}(z)$的最大插值误差，多次重复这一过程，每次只保留最大误差，最后以获得的子矩阵块的最大插值误差近似表示${\pmb T}(z)$的最大插值误差.将上述结果绘制成误差曲线，见图 4.可以看出，当$r = R$时，$t(z)$的Cauchy插值误差描述了对${\pmb T}(z)$进行Cauchy插值的最大误差，即可以通过函数$t(z) = f(z){\rm e}^{{\rm i}zR}$的Cauchy插值来估计对${\pmb T}(z)$进行Cauchy插值所需的插值项数.
观察图 4可知，当$d = 20$时，矩阵插值精度达到BEM快速算法精度，即插值项数满足要求. RSRR算法的其余参数分别设置为$N = 50$，$L = 2$，采样点根据二维Chebyshev插值规则分布于围道内，计算该模型的固有频率和振型.将特征值计算结果与文献[25]相同算例结果进行对比，见表 1.可以看出，本文RSRR方法的计算结果和文献[23]中的结果是一致的，但是本文的精度较高.该算例共计耗时41.4 min，占用内存344.3 MB.
表 1 关注频段内的全部固有频率和误差Table 1. Natural frequencies and error in the spectrum in attention
为验证$d$的估计方法针对RSRR算法的合理性，保持其他参数不变，对比$d = 10, \; 16, \; 20, \; 26, \; 30$时各特征对的误差，见图 5.可以看出当$d= 20$时，特征对的精度正好达到最优，即本文提出的确定$d$取值的方法针对RSRR算法是合理的.
4.3 大规模工程应用算例
将本文提出的RSRR非线性特征值解法应用于实际工程中大规模吸声结构的复模态分析.考虑火箭整流罩内声腔模型，其半剖边界元模型见图 6.该模型轴向最大尺寸为2.29 m，径向最大尺寸为1.324 m，将其离散为64 162个6节点非协调三角形单元，共计384 972自由度，全部表面设定为阻抗边界，阻抗模型与算例4.2相同.
本算例欲求矩形区域(200+0i, 500+100i)Hz内的所有特征解，取$C$为该矩形.考虑阻抗边界，给定边界元核函数$t(z) = Z \cdot z \cdot {\rm e}^{{\rm i}zR}$，插值点取为围道上的高斯积分点，$d$依次取$10, \; 16, \; 20, \; 26, \; 30$，根据式(33) 计算$t(z)$和${\pmb T}(z)$在不同插值项数下的Cauchy插值误差，见图 7，${\pmb T}(z)$的近似方法同4.2节.可以看出，$t(z) = f(z){\rm e}^{{\rm i}zR}$的Cauchy插值误差能够反映出对${\pmb T}(z)$每个元素进行Cauchy插值的最大误差，且当$d$取26时，矩阵近似精度已满足要求.本算例取$d = 30$. RSRR算法的其余参数分别设置为$N = 50$，$L = 4$，采样点根据二维Chebyshev插值规则分布于关注频段内，计算该模型的固有频率和振型.
![]()
图 7 基于Cauchy插值的BEM矩阵和核函数最大误差曲线Figure 7. Max relative error of BEM matrix and Kernel function based on the Cauchy interpolation该算例耗时122.3 h，占用内存21.023 GB.共计算出51组特征对，其相对误差分布见图 8，前四阶固有频率以及对应的声压模态见图 9.该算例显示了RSRR算法对大规模问题仍然保持了较高的计算精度.
5. 结论
本文提出了一种解决大规模边界元特征值问题的数值方法RSRR，为工程中的大规模边界元模态分析提供了有效计算手段.
(1) 基于解析矩阵函数的Keldysh定理，提出了一种简单易行、可靠性高的特征向量搜索空间构造方法，该方法直接采用一系列边界元频域解的一组基空间构造特征空间，无需对边界元系数矩阵取高阶矩，解决了传统围道积分方法存在的特征空间易出现病态甚至构造失败的问题.
(2) 基于解析函数的Cauchy积分公式，构造了复平面内边界元系数矩阵的插值方法，以及缩减特征值问题系数矩阵的快速计算方法，还根据核函数与系数矩阵的关系给出了插值项数的估计策略.该方法在降低RSRR方法求解边界元特征值问题计算量的同时，也将其扩展到了工程中常见的边界元复特征值分析中.
(3) 将所提出的RSRR方法与声学快速边界元结合，在普通工作站上实现了复杂边界条件下自由度接近40万的边界元声模态分析，显示了本文方法的大规模计算分析能力.
-
![]()
图 1 围道$C$示意图.图中闭合轮廓线代表围道$C$，"+"代表$C$内的待求特征值，"$\circ$"代表$C$外的特征值，"$\bullet$"代表采样点
Figure 1. A diagram showing the contour $C$ (closed contour). eigenvalues interested (+), eigenvalues not interested ($\circ$) and sampling points in the contour ($\bullet$)
![]()
图 4 基于Cauchy插值的BEM矩阵和核函数最大误差曲线
Figure 4. Max relative error of BEM matrix and Kernel function based on the Cauchy interpolation
![]()
图 7 基于Cauchy插值的BEM矩阵和核函数最大误差曲线
Figure 7. Max relative error of BEM matrix and Kernel function based on the Cauchy interpolation
-
[1] 姚振汉, 王海涛.边界元法.北京:高等教育出版社, 2010 Yao Zhenghan, Wang Haitao. Boundary Element Method. Beijing:Higher Education Press, 2010 (in Chinese)
[2] 高效伟, 刘健, 彭海峰.集成单元边界元法及其在主动冷却热防护系统分析中的应用.力学学报, 2016, 48(4):994-1003 http://lxxb.cstam.org.cn/CN/abstract/abstract145937.shtml Gao Xiaowei, Liu Jian, Peng Haifeng. Integrated unit BEM and its application in analysis of actively cooling TPS. Chinese Journal of Theoretical and Applied Mechnics, 2016, 48(4):994-1003 (in Chinese) http://lxxb.cstam.org.cn/CN/abstract/abstract145937.shtml
[3] 申建伟, 孙中奎, 詹世革等.第9届全国动力学与控制青年学者学术研讨会报告综述.力学学报, 2015, 47(6):1079-1083 doi: 10.6052/0459-1879-15-366 Shen Jianwei, Sun Zhongkui, Zhan Shige, et al. Review of the ninth national symposium on dynamics and control for young scholars. Chinese Journal of Theoretical and Applied Mechanics, 2015, 47(6):1079-1083 (in Chinese) doi: 10.6052/0459-1879-15-366
[4] 荣俊杰, 校金友, 文立华.弹性动力学高阶核无关快速多极边界元法.力学学报, 2014, 46(5):776-785 http://lxxb.cstam.org.cn/CN/abstract/abstract144819.shtml Rong Junjie, Xiao Jinyou, Wen Lihua. A high order kernel independent fast multipole boundary element method for elastodynamics. Chinese Journal of Theoretical and Applied Mechanics, 2014, 46(5):776-785 (in Chinese) http://lxxb.cstam.org.cn/CN/abstract/abstract144819.shtml
[5] Liu YJ, Mukherjee S, Nishimura N, et al. Recent advances and emerging applications of the boundary element method. Applied Mechanics Reviews, 2011, 64(3):030802 http://appliedmechanicsreviews.asmedigitalcollection.asme.org/article.aspx?articleid=1678921
[6] Tisseur F, Meerbergen K. The quadratic eigenvalue problem. SIAM Review, 2001, 43(2):235-286 doi: 10.1137/S0036144500381988
[7] Mehrmann V, Schröder C. Nonlinear eigenvalue and frequency response problems in industrial practice. Journal of Mathematics in Industry, 2011, 1(1):7 doi: 10.1186/2190-5983-1-7
[8] Effenberger C. Robust solution methods for nonlinear eigenvalue problems.[PhD Thesis]. École polytechnique fédérale de Lausanne, 2013 doi: 10.5075/epfl-thesis-5920
[9] Van Beeumen, Rational Krylov methods for nonlinear eigenvalue problems[PhD Thesis], KU Leuven, 2015 http://lirias.kuleuven.be/handle/123456789/487915
[10] Ali A, Rajakumar C, Yunus SM. Advances in acoustic eigenvalue analysis using boundary element method. Computers & Structures, 1995, 56(5):837-847 http://www.sciencedirect.com/science/article/pii/0045794995000126
[11] Bezine G. A mixed boundary integral-finite element approach to plate vibration problems. Mechanics Research Communications, 1980, 7(3):141-150 doi: 10.1016/0093-6413(80)90003-8
[12] Nardini D, Brebbia CA. A new approach to free vibration analysis using boundary elements. Applied Mathematical Modelling, 1983, 7(3):157-162 doi: 10.1016/0307-904X(83)90003-3
[13] Ali A, Rajakumar C, Yunus SM. On the formulation of the acoustic boundary element eigenvalue problems. International Journal for Numerical Methods in Engineering, 1991, 31(7):1271-1282 doi: 10.1002/(ISSN)1097-0207
[14] Kamiya N, Andoh E. Standard eigenvalue analysis by boundaryelement method. International Journal for Numerical Methods in Biomedical Engineering, 1993, 9(6):489-495 doi: 10.1002/cnm.1640090606/pdf
[15] Asakura J, Sakurai T, Tadano H, et al. A numerical method for nonlinear eigenvalue problems using contour integrals. JSIAM Letters, 2009, 1:52-55 doi: 10.14495/jsiaml.1.52
[16] Beyn WJ. An integral method for solving nonlinear eigenvalue problems. Linear Algebra and Its Applications, 2012, 436(10):3839-3863 doi: 10.1016/j.laa.2011.03.030
[17] Sakurai T, Asakura J, Tadano H, et al. Error analysis for a matrix pencil of Hankel matrices with perturbed complex moments. JSIAM Letters, 2009, 1:76-79 doi: 10.14495/jsiaml.1.76
[18] Yokota S, Sakurai T. A projection method for nonlinear eigenvalue problems using contour integrals. JSIAM Letters, 2013, 5:41-44 doi: 10.14495/jsiaml.5.41
[19] Zheng CJ, Chen HB, Gao HF, et al. Is the Burton-Miller formulation really free of fictitious eigenfrequencies? Engineering Analysis with Boundary Elements, 2015, 59:43-51 doi: 10.1016/j.enganabound.2015.04.014
[20] Gao HF, Matsumoto T, Takahashi T, et al. Eigenvalue analysis for acoustic problem in 3D by boundary element method with the block Sakurai-Sugiura method. Engineering Analysis with Boundary Elements, 2013, 37(6):914-923 doi: 10.1016/j.enganabound.2013.03.015
[21] Leblanc A, Lavie A. Solving acoustic nonlinear eigenvalue problems with a contour integral method. Engineering Analysis with Boundary Elements, 2013, 37(1):162-166 doi: 10.1016/j.enganabound.2012.09.007
[22] Xiao JY, Meng SS, Zhang CZ, et al. Resolvent sampling based Rayleigh-Ritz method for large-scale nonlinear eigenvalue problems. Computer Methods in Applied Mechanics and Engineering, 2016, 310:33-57 doi: 10.1016/j.cma.2016.06.018
[23] Xiao JY, Zhou H, Zhang CZ, et al. Solving large-scale finite element nonlinear eigenvalue problems by resolvent sampling based Rayleigh-Ritz method. Computational Mechanics, 2016:1-18 http://dl.acm.org/citation.cfm?id=3054532
[24] Mehrmann V, Voss H. Nonlinear eigenvalue problems:A challenge for modern eigenvalue methods. GAMM-Mitteilungen, 2004, 27(2):121-152 doi: 10.1002/gamm.v27.2
[25] Leblanc A, Lavie A. Numerical analysis of eigenproblem for cavities by a particular integral method with a low frequency approximation of surface admittance. The Journal of the Acoustical Society of America, 2012, 131(5):3876-3882 doi: 10.1121/1.3699270
[26] Du JT, Li WL, Liu ZG, et al. Acoustic analysis of a rectangular cavity with general impedance boundary conditions. The Journal of the Acoustical Society of America, 2011, 130(2):807-817 doi: 10.1121/1.3605534
[27] 屈伸, 陈浩然.敷设多孔吸声材料声腔特征值分析的径向积分边界元法.计算力学学报, 2015, 32(1):123-128 doi: 10.7511/jslx201501021 Qu Shen, Chen Haoran. Eigenvalue analysis for acoustical cavity covered with porous materials by using the radial integration boundary element method. Chinese Journal of Computational Mechanics, 2015, 32(1):123-128 (in Chinese) doi: 10.7511/jslx201501021
[28] 陈文炯, 刘书田.周期性吸声多孔材料微结构优化设计.计算力学学报, 2013, 30(1):45-50 doi: 10.7511/jslx201301008 Chen Wenjiong, Liu Shutian. Optimizing design of micro-structural configurations of periodic porous soundabsorbing materials. Chinese Journal of Computational Mechanics, 2013, 30(1):45-50(in Chinese) doi: 10.7511/jslx201301008
[29] Cao YC, Wen LH, Xiao JY, et al. A fast directional BEM for largescale acoustic problems based on the Burton-Miller formulation. Engineering Analysis with Boundary Elements, 2015, 50:47-58 doi: 10.1016/j.enganabound.2014.07.006
[30] Engquist B, Ying Lexing. Fast directional multilevel algorithms for oscillatory kernels. SIAM Journal on Scientific Computing, 2007, 29(4):1710-1737 doi: 10.1137/07068583X
[31] Rong J, Wen L, Xiao J. Efficiency improvement of the polar coordinate transformation for evaluating BEM singular integrals on curved elements. Engineering Analysis with Boundary Elements, 2014, 38:83-93 doi: 10.1016/j.enganabound.2013.10.014
[32] Gohberg I, Rodman L. Analytic matrix functions with prescribed local data. Journal d'Analyse Mathématique, 1981, 40(1):90-128 doi: 10.1007/BF02790157
[33] Austin AP, Kravanja P, Trefethen LN. Numerical algorithms based on analytic function values at roots of unity. SIAM Journal on Numerical Analysis, 2014, 52(4):1795-1821 doi: 10.1137/130931035
下载:
计量
- 文章访问数: 1288
- HTML全文浏览量: 205
- PDF下载量: 1031
