王俊鵬 校金友 文立華

(西北工業(yè)大學航天學院，西安710072)

大規(guī)模邊界元模態(tài)分析的高效數值方法1)

王俊鵬 校金友2)文立華

(西北工業(yè)大學航天學院，西安710072)

隨著大規(guī)?？焖龠吔缭嬎慵夹g的發(fā)展，在復雜結構的動態(tài)設計、振動與噪聲分析中愈來愈多地采用邊界元法，因此求解大規(guī)模邊界元特征值問題、進行復雜結構和聲場模態(tài)分析，成為工程應用中一個十分重要，但卻極具挑戰(zhàn)性的課題，目前國際上還沒有十分有效的數值方法.本文針對邊界元法中典型的非線性特征值問題，提出了一種通用、高效的數值解法，稱為基于預解矩陣采樣的Rayleigh-Ritz投影法，記為RSRR.首先，通過求解一系列頻域邊界元問題來構造特征向量搜索空間，進而可以采用Rayleigh-Ritz投影，將原問題轉化為一個可以采用現(xiàn)有方法求解的小規(guī)?？s減特征值問題;其次，為了降低Rayleigh-Ritz投影過程的計算量，基于解析函數的Cauchy積分公式，構造了邊界元系數矩陣的插值近似方法，以及縮減特征值問題系數矩陣的快速計算方法，給出了插值項數的估計策略;最后，將RSRR與聲學快速邊界元法結合，應用于大規(guī)模吸聲結構的復模態(tài)分析.數值算例表明，RSRR方法能夠可靠地求出給定頻段內的全部特征值和特征向量，具有計算效率高、精度高、通用等優(yōu)點.

邊界元法，非線性特征值問題，瑞利--里茨投影，柯西積分，模態(tài)分析

引言

邊界元法[1](boundary element method，BEM)是一種重要的工程數值方法，在計算力學、聲學等領域應用很廣[23].尤其是對于中高頻振動和波的傳播問題[4]，邊界元法不僅能有效降低求解規(guī)模，還能避免數值色散誤差，具有很大的發(fā)展?jié)摿?近年來，邊界元快速算法逐步發(fā)展成熟，顯著提高了邊界元大規(guī)模分析的計算效率，將邊界元法向實際工程應用推進了一大步[5].

邊界元法用于結構減振降噪分析，常常需要進行模態(tài)分析，以獲得系統(tǒng)的固有頻率和模態(tài).邊界元模態(tài)分析的難點在于求解如下形式的非線性特征值問題

式中，λ和v分別為特征值和與之對應的特征向量，系統(tǒng)矩陣T(z)∈Cn×n是聲波數z=ω/c的函數，ω是聲波的圓頻率，c是聲速，n表示系統(tǒng)總自由度數.特別地，當T(z)為z的線性或二次函數時，上述問題即為標準特征值問題或二次特征值問題；后者是最簡單的一類非線性特征值問題，普遍存在于有限元結構分析領域，其數值解法已經較為成熟[6].

在邊界元法中，矩陣T(z)關于z的表達式十分復雜，甚至無法得到.對這類非線性特征值問題，數學上還沒有一個通用、有效、適合大規(guī)模求解的數值解法[79].加之邊界元矩陣T(z)還是一個非對稱的滿陣，計算量和存儲量都很大，過去采用的行列式搜索法[10]、域內離散法[11]、域積分轉換法[1214]等都無法用于大規(guī)模求解.

圍道積分法(contour integral method，CIM)是近年來才提出的一種求解非線性特征值問題的數值解法[1517]，最初被用于求解廣義特征值問題[18].它不需要迭代，能一次計算出復平面給定區(qū)域內的全部特征值和相應特征向量；同時還適合于并行求解，計算效率很高，具有廣闊的應用前景.目前該方法已被用于簡單邊界元特征值問題的求解[1920].

但是，CIM可能導致漏根或是虛假特征根，求解精度難以保證.這是因為它采用矩陣T(z)的高階矩來構造特征向量空間[1921].矩的階數越高，這些矩陣越趨于線性相關，從而導致特征空間病態(tài)或構造失敗.另外，按照CIM的思想，頻率采樣點只能落在包含尋根區(qū)域的某個圍道上，而實際計算結果表明，采用域內更為靠近特征值的采樣點，求解精度會更高[22].

為解決上述問題，最近我們發(fā)展了一種穩(wěn)定、高效的非線性特征值問題數值新解法，稱為基于預解矩陣采樣的 Rayleigh-Ritz投影法，簡記為RSRR(resolvent sampling based rayleigh-ritz projection method)[22].它不僅可以構造出可靠的特征空間、自由布置頻率采樣點以提高求解精度，而且具備CIM能一次求出指定范圍內所有特征值、適合并行計算的優(yōu)點.目前，RSRR已被成功用于有限元和一類簡單的邊界元大規(guī)模非線性特征值問題的求解，在精度、計算效率和穩(wěn)定性等方面均顯示出較大優(yōu)勢[22-23].

本文把RSRR應用于工程邊界元減振降噪分析中，解決Rayleigh-Ritz投影法中邊界元矩陣T(z)的快速計算問題.該問題是制約大規(guī)模邊界元特征值問題求解的瓶頸[24]，因為求解Rayleigh-Ritz投影所得的縮減特征值問題，往往需要成百上千步的迭代或矩陣求逆，而每一步都需要計算一個稠密矩陣T(z).對于自由度幾十萬、上百萬的大規(guī)模問題，其計算量是普通計算設備無法承受的.

本文提出一種基于Cauchy積分公式的邊界元矩陣插值方法，可以在復平面內任意區(qū)域上快速構造Rayleigh-Ritz縮減矩陣的高精度插值表達式，從而避免了在求解縮減特征值問題中計算原矩陣T(z)的難題,大大降低了RSRR求解邊界元特征值問題的計算量；還給出了一種簡便的插值項數確定方法，可根據精度要求確定一個近似最少的插值項數，有效控制整個求解過程的計算量.與文獻[22]中的一維Chebyshev插值法相比，本文的Cauchy插值法可在復平面內任意二維區(qū)域上構造插值近似，可用于解決高阻尼、高吸聲減振降噪設計中常見的復特征值問題[25-28].

1 邊界元特征值問題

本文考慮吸聲結構分析中的邊界元特征值問題，但所提之方法具有通用性，也適用于彈性動力學、電磁學等問題.

聲學問題的邊界積分方程為

式中，D為聲場區(qū)域，?D為D的邊界，p為聲壓，F(xiàn)(x,y)=exp(iz|x?y|)/|x?y|是三維 Helmholtz方程的基本解，c(x)為自由項，當?D在x附近光滑時，c(x)=1/2，ny表示y點處的外法向.

采用文獻[29]的Nystr¨om邊界元法離散式(2)，選擇核無關快速定向壓縮算法[30]對其進行加速.將聲場邊界?D離散為二次曲面三角形單元，用6節(jié)點二次協(xié)調元進行幾何近似，用6節(jié)點二次非協(xié)調元進行函數近似.采用非協(xié)調元可以避免邊界元角點問題[4]，方便幾何建模，有益于方法的通用性.同時，我們的計算經驗表明，相同自由度數目情況下，采用非協(xié)調元可以獲得與協(xié)調元相同的精度.經Nystr¨om方法離散邊界積分方程(2)，最終可得如下標準形式的邊界元離散方程

式中，H和G是邊界元離散矩陣，p和q分別是包含節(jié)點聲壓及其法向導數的列向量.矩陣H和G中的近場修正部分仍需要進行常規(guī)邊界元法中的奇異和近奇異積分計算.本文奇異積分采用文獻[31]中的數值計算方法，而近奇異積分采用逐步細分的方法來計算.

考慮Robin邊界條件

當α≠0,β=0時，式(4)退化成Dirichlet邊界條件，當α=0,β≠0時，式(4)退化成Neumann邊界條件.

下面考慮吸聲結構分析中常見的一類Robin邊界條件，即阻抗邊界條件

式中，Z為邊界聲阻抗，與節(jié)點位置 x和頻率z有關；ρ0為聲場介質密度，式(5)在每個邊界元節(jié)點xj上都成立，將這些式子統(tǒng)一寫成矩陣形式即為

式中，B(z)是維度為n×n的對角矩陣，其對角線元素表示為 ?iρ0ω/Z(xj,z).將式 (6)代入式 (3)，可得

上式即為吸聲結構模態(tài)分析中所要求解的特征方程，系統(tǒng)矩陣T(z)=H(z)?G(z)B(z).Dirichlet和Neumann邊界條件所對應的特征方程可類似得到.

邊界元法中，矩陣T(z)的元素是由基本解或其法向導數經單元積分后得到的.由于基本解是頻率z的復雜函數，加之單元積分通常只能靠數值方法計算，而無法獲得解析表達式，因此T(z)關于z的解析表達式通常無法得到，而只能對于給定的z通過數值方法來計算.這不僅增大了邊界元特征值問題求解的難度，而且還導致其計算量很大，以至于目前仍沒有有效的大規(guī)模數值解法，本文正是要解決這些問題.

2 非線性特征值問題的RSRR解法

RSRR算法是基于Rayleigh-Ritz投影的特征值解法.給定一個包含待求解特征向量的搜索空間S，設S的一組正交基為Q，則原非線性特征值問題(1)可以通過Rayleigh-Ritz投影轉化為如下的縮減形式

式中，TQ(z)=QHT(z)Q ∈CkS×kS，kS為Q的列向量數目,上標H表示共軛轉置.式(8)與式(1)具有相同的特征值，設(λ,g)為TQ(z)的任一特征對，則對應原問題的特征對為(λ,Qg).下面首先介紹搜索空間S的構造方法，然后給出RSRR算法的基本步驟.

2.1 特征空間的構造

本文的預解矩陣采樣法 (resolventsampling method，RES)是構造特征空間S的一種穩(wěn)定、高效的數值方法.設要計算復平面上由簡單閉合曲線C圍成的區(qū)域內的所有特征值和相應特征向量，用RES構造特征空間可分為以下兩步：

(1)在圍道C(見圖 1)內選擇 N個采樣點 zk，k=1,2,···,N，同時生成L個線性無關的n維隨機向量，組成矩陣U∈Cn×L；

(2)對 k=1,2,···,N，計算 T(zk)?1U，并將其組裝成矩陣S∈Cn×NL

矩陣S的列空間span(S)即為RES所形成的特征向量搜索空間.

圖1 圍道C示意圖.圖中閉合輪廓線代表圍道C，“+”代表C內的待求特征值，“°”代表C外的特征值，“·”代表采樣點Fig.1 A diagram showing the contour C(closed contour).eigenvalues interested(+),eigenvalues not interested(°)and sampling points in the contour(·)

下面首先證明：當采樣點布置、參數L和N選取都合理的情況下，空間span(S)將近似包含圍道C內的所有特征向量，在此基礎上討論采樣點和隨機矩陣U對RES方法效率的影響情況.

邊界元矩陣T(z)雖然形式復雜，但由于基本解是頻率z的全純函數(holomorphicfunction)，因此T(z)也是z的全純矩陣值函數(holomorphic matrix-value function).根據解析矩陣函數的Keldysh定理[16]，邊界元矩陣T(z)的逆矩陣T?1(z)可由C內的特征值和特征向量表示為如下形式

式中，nC為不計重根的互不相同特征值的數目，ηk為T(λk)的幾何重數，μl,k表示λk的局部重數，和分別是T(λk)的右、左正規(guī)化廣義特征向量，RC(z)表征C之外的所有特征值的貢獻，詳細理論可參考文獻[16,32].

將式(10)的右端項寫成矩陣形式

其中，Φ(z)是由 T(z)特征值的 Jordon塊組成的矩陣，VC和WC分別是矩陣T(z)和T(z)H在特征根λ處的廣義特征向量標準系統(tǒng)(canonical system of generalized eigenvectors).為了從T(z)?1中提取出特征向量空間span(VC)，必須盡可能地消除RC(z)的影響.由于RC(z)在特征值附近解析，因此可以用一組合適的基函數{gj(z)}來逼近它

式中，Rj為展開系數矩陣，J為展開項數.當{gj(z)}選擇合適時，式(12)的逼近過程可以很快收斂.例如，對于區(qū)間上的解析函數，用Chebyshev多項式逼近可以獲得指數收斂速度.

現(xiàn)對T(z)?1右乘以包含L個線性無關向量的矩陣U，并將式(12)代入式(11)，整理可得

式中 Dj(z)=diag(gj(z),gj(z),···,gj(z))∈ CL×L.

由式(14)可知，VC中的特征向量包含在的列空間中，即

以上推導表明，當采樣點的布置、參數L、N和J取值合理的情況下，式(9)定義的采樣矩陣S所張成的子空間會包含所有待求的特征向量.同時也表明，這幾個參數對RES方法的效率有重要影響.但在實際應用中，這些參數的理想值事先很難確定.下面通過分析，給出這些參數的選取原則.

由上面推導過程可見，采樣點zk的作用是形成對未知矩陣 RC(z)的插值，從而隱式地確定式(12)中的基函數{gj(z)}.據此，對于實特征值問題，可將采樣點取為尋根區(qū)間上的Chebyshev點，這對應于Chebyshev插值，精度較高；而對復特征值問題，圍道C通常為包含尋根區(qū)域的矩形，此時采樣點可取為該矩形區(qū)域上的二維Chebyshev點.

采樣點數目N和隨機向量數目L共同決定著RES方法的精度和計算效率.就保證特征空間精度而言，大量計算經驗表明，為滿足式(18)一般需要N·L＞2ˉnC～3ˉnC.同時，U中的L個隨機向量的作用是探測多重特征值所對應的特征向量，因此L的取值應該不小于T在C內各特征值的最大局部重數[16]，即

文獻[21]研究表明，相比于L，增加N對提高精度的效果更顯著.但是另一方面，由式(9)可知，計算矩陣S需要求解N個線性方程組，每個包含L個右端向量，因此增大N和L都會大幅增加計算量.綜上，一般在滿足式(18)的條件下，N和L應盡量?。欢=(2ˉnC～3ˉnC)/L.特征值數目ˉnC并不要求準確地知道，因此可利用粗網格下的小規(guī)模問題進行預估，也可通過理論分析或參考相似結構的已知結果而獲得.

2.2 RSRR算法

RSRR算法求解大規(guī)模邊界元非線性特征值問題的步驟如下：

(1)確定圍道C、采樣點數 N、采樣點zk,(k=1,2,···,N)、重數 L、隨機矩陣 U ∈ Cn×L、展開項數d；

(2)計算 T(zk)?1U,k=1,2,···,N；

(3)組裝S矩陣，并對其進行截斷奇異值分解S≈QΣWH，得到一組正交基Q；

(4)計算縮減矩陣 TQ(z)=QHT(z)Q的插值近似，計算方法將在第3節(jié)詳述；

(5)求解縮減后的非線性特征值問題TQ(z)g=0，設所得的特征對為(λj,gj)，j=1,2,···,ˉnC，則原問題的特征對為(λj,vj),j=1,2,···,ˉnC，vj=Qgj.

在第三步的Rayleigh-Ritz投影中，計算正交基Q的截斷奇異值分解的奇異值舍去閾值一般設為δ=10?14，而Q的列數ks即為保留的奇異值數目.為了保證求出積分路徑內的所有特征根，ks應滿足

奇異值分解的計算復雜度為O(n)，且一般ks僅為幾十到幾百，因此奇異值分解的計算量在整個邊界元特征值分析中可忽略不計.

實際應用中，會出現(xiàn)某些采樣點zk靠近特征值的情況，此時T(zk)?1U的數值遠大于S中其他向量，從而導致特征空間的總體精度下降，因此在計算正交基Q之前需要對S各列進行正規(guī)化.

3 基于Cauchy積分的矩陣插值方法

經過Rayleigh-Ritz投影后的非線性特征值問題規(guī)模很小，見式(8)，現(xiàn)有的許多方法都可以對其進行求解[1516,22].本文采用文獻 [22]改進的 Block Sakurai-Sugiura(Block SS)算法，該算法能夠自動確定待求特征值的數目，并且計算出給定區(qū)域內的全部特征值及對應特征向量.但是，無論是Block SS算法還是其他方法，求解式(8)都需要計算多個給定頻率zi處的系數矩陣TQ(zi)=QHT(zi)Q.由于邊界元矩陣T(zi)是大規(guī)模稠密矩陣，這部分計算量將十分巨大，是制約邊界元大規(guī)模非線性特征值問題求解的一個瓶頸.

3.1 矩陣近似

為有效降低多次計算TQ(zi)的計算量，文獻[22]通過Chebyshev插值方法獲得了TQ(z)的Chebyshev插值近似格式，在求解TQ(z)g=0時只需要調用小規(guī)模Chebyshev插值式計算TQ(z)，從而避免了反復計算邊界元矩陣T(z).但是，由Chebyshev插值理論可知，此方法只有在特征值分布在直線附近時(如實特征值問題)才能保證足夠的精度.而一般阻尼吸聲結構的特征值分布在二維復平面內，該方法不再適用.

為克服Chebyshev插值方法的局限性，本文基于Cauchy積分法構造TQ(z)在復平面內任意單連通區(qū)域上的插值近似格式.由于邊界元矩陣T(z)為全純函數，因此可由Cauchy積分公式表示為

式中，C為復平面內任意積分路徑，本文中它與第2節(jié)的RSRR算法中的圍道取為一致.采用數值積分計算上式的圍道積分，可得

這里沿用文獻[30]中對標量函數的叫法，將式(22)稱為矩陣函數T(z)的第一類Cauchy插值公式.其中，總項數d決定了數值積分的精度，wi為積分點zi對應的權系數，它們由所采用的數值積分公式確定.對于光滑的閉合曲線，通常采用梯形公式可獲得指數階的計算精度.

由于是在復平面內進行矩陣插值，Cauchy插值公式比Chebyshev插值公式適用范圍更廣.但是，當z逼近zi時，1/(zi?z)趨于無窮大，式(22)的精度會迅速下降，導致插值失敗.針對這一現(xiàn)象，用(T(ζ)? T(z))/(ζ? z)代替式 (21)中的 T(ζ)/(ζ? z)，前者在z逼近ζ時分子分母同時趨于0，不存在奇異性，且其Cauchy積分值為0，即

上式用數值積分公式離散后可得

整理上式，得到T(z)的表達式如下

稱上式為第二類Cauchy插值公式，式中

將式(25)代入TQ(zi)=QHT(zi)Q即可得到縮減特征值問題系數矩陣的近似式

其中，系數矩陣Ci=QHT(zi)Q為kS維方陣，規(guī)模很小，可以直接計算并存儲起來.Ci的主要計算量是在d個插值點zi上計算邊界元矩陣T(zi)，并通過矩陣向量乘法計算T(zi)Q.為有效控制計算量，需使d盡量小，第3.2節(jié)將提出d的估計方法.同時，由于d個插值點是預先給定的，因此d個系數矩陣Ci很適合在不同處理器上并行計算.

相較第一類 Cauchy插值公式 (22)，第二類Cauchy插值公式(25)不僅消除了當z逼近zi時因奇異性導致的插值不穩(wěn)定現(xiàn)象，還能顯著提高插值精度.文獻[33]稱這一現(xiàn)象為“過收斂”現(xiàn)象，即式(25)不僅在積分區(qū)域C內滿足精度要求，在邊界C附近也依然具有較高的精度.

下面通過一個邊界元矩陣插值算例來展示“過收斂”現(xiàn)象，證明第二類Cauchy插值公式可以為本文的RSRR算法提供更充足的矩陣近似精度.選擇算例4.2的邊界元矩陣和頻率區(qū)域，展開項數d取20，矩陣插值的相對誤差按下式計算

圖2 Cauchy插值相對誤差Fig.2 Error of the Cauchy interpolation

其中T?Q(z)表示直接計算出的縮減邊界元矩陣，y為隨機向量.第一類Cauchy插值和第二類Cauchy插值的相對誤差云圖見圖2，圖中“°”代表積分路徑上的20個插值點.可見，第一類Cauchy插值公式(22)僅在區(qū)域中心能夠保持一定精度，靠近積分邊界時，其插值精度迅速惡化.而第二類Cauchy公式(25)不僅在積分區(qū)域內能夠保持高精度，在積分區(qū)域邊界及其外部一定范圍內也仍然呈現(xiàn)出很高的精度，明顯優(yōu)于第一類Cauchy公式，能夠滿足RSRR方法對矩陣近似精度的要求.

圖2 Cauchy插值相對誤差(續(xù))Fig.2 Error of the Cauchy interpolation(continued)

3.2 確定展開項數

在保證矩陣近似精度的前提下，展開項數d越小，近似矩陣的計算量越小.這里提出一種估計最小展開項數d的方法.根據3.1節(jié)的Cauchy矩陣插值方法，矩陣T(z)的插值公式(25)同樣適用于對矩陣每個元素的插值近似，因此可以用矩陣元素插值所需的最小項數快速地估計d.

一般地，邊界元法系數矩陣T(z)的元素t(z)具有如下通式

其中，r表示源點和場點的距離，R代表模型最大尺寸，在所關注頻段內 f(z)為z的解析函數.例如，采用配點法離散邊界積分方程(2)，則H矩陣的元素可寫成

式中，?代表邊界單元，w(y)代表邊界元形函數，函數 eiz|x?y|可展開成

其中，rk=|x?yk|,φk是展開基函數.故

當r=R時，函數t(z)振蕩最為劇烈，對展開項數要求最為苛刻.因此，本文用函數t(z)=f(z)eizR的Cauchy插值來估計對T(z)進行Cauchy插值所需的插值項數.

按式(25)計算函數t(z)在圍道C內的Cauchy插值近似，定義近似函數?t(z,d)插值的L2誤差為

假設矩陣T(z)的插值精度要求為ε0，則可以將滿足ε(d0)≤ε0的最小整數d0作為T(z)所需的插值項數的估計，即d=d0.

實際應用中，也可以通過以下3個途徑獲得更為可靠的d的估計：

(2)選取多個函數t(z)=f(z)eizr分別進行Cauchy插值，d的估計值應該對每個函數都滿足插值精度要求.一般情況下，可以選擇多個不同的0≤rk≤R來得到一系列函數tk(z)=f(z)eizrk.

(3)在估計值 d0的基礎上給定一個安全系數Cd，即令由于 Cauchy 插值通常是指數收斂的，因此系數Cd通常不必太大，如1≤Cd≤1.2.

4 算例

本文所有算例均在一臺 32核 (Xeon E5-2630 v3，2.40GHz)64GB內存的工作站上進行，所有程序都采用單線程編寫.邊界元法快速算法的精度ε設定為1.0×10?5，線性方程組采用廣義極小殘量法(GMRES)迭代求解，收斂殘差為1.0×10?6，并采用對角塊預條件技術.

4.1 解析解對比算例

該算例旨在驗證本文方法的正確性，并與現(xiàn)有方法進行計算精度和效率上的對比.

考慮單位立方體內聲腔模型，將其表面離散為1200個三角形單元，邊界元自由度數為7200，選擇一組相互平行的表面設定為Dirichlet邊界，其余4個表面設定為Neumann邊界.

圖3 特征值誤差圖Fig.3 Relative error of eigenvalues

由圖4可以看出，RSRR方法計算精度明顯高于文獻[18]的圍道積分法.同時本算例運用RSRR算法只需生成N=20次BEM系數矩陣，求解N×L=80個BEM方程組.而文獻[20]要生成256次BEM系數矩陣，求解256×15=3840個BEM方程組，求解邊界元方程的數目是RSRR的近50倍，但是求解精度卻比RSRR低了近4個量級！

圖4 基于Cauchy插值的BEM矩陣和核函數最大誤差曲線Fig.4 Max relative error of BEM matrix and Kernel function based on the Cauchy interpolation

4.2 復雜阻抗算例

考慮一長方體內聲腔模型，其長、寬、高分別為 0.5m、0.4m和 0.3m，將其表面離散為1256個三角形單元，邊界元自由度數為 7536.在其中一個 0.4m×0.3m的表面上敷設多孔吸聲材料，等效為Delany-Bazley阻抗邊界條件[25]，阻抗表達式如下，其余表面為 Neumann邊界，關注頻率區(qū)域為(300+0i,1000+200i)Hz，圍道C取為該矩形.

式中，δ=ρ0ω/2πσ，σ表示流阻，h和φ分別表示多孔吸聲材料層的厚度和孔隙率，ρ0為多孔材料密度.本文取 σ =2×104N·s/m4,h=20cm,φ=1和ρ0=1.2kg/m3.

首先，確定最優(yōu)插值項數d.考慮阻抗邊界，給定邊界元核函數的表達式t(z)=Z·z·eizr，分別取r=R,0.8R,0.6R,0.4R，插值點取為圍道上的高斯積分點，插值項數d依次取10,16,20,26,30，根據式(33)計算不同插值項數下t(z)和T(z)的插值誤差.由于T(z)為非對稱滿陣，直接對其進行Cauchy插值計算量太大，難以實現(xiàn)，此處采用了一個近似等效的方法：隨機取200個邊界元節(jié)點，計算這些節(jié)點在T(z)中對應的200×200子矩陣塊，用此子矩陣塊代替T(z)進行插值，取其最大插值誤差代替T(z)的最大插值誤差，多次重復這一過程，每次只保留最大誤差，最后以獲得的子矩陣塊的最大插值誤差近似表示T(z)的最大插值誤差.將上述結果繪制成誤差曲線，見圖4.可以看出，當r=R時，t(z)的Cauchy插值誤差描述了對T(z)進行Cauchy插值的最大誤差，即可以通過函數t(z)=f(z)eizR的Cauchy插值來估計對T(z)進行Cauchy插值所需的插值項數.

觀察圖4可知，當d=20時，矩陣插值精度達到BEM快速算法精度，即插值項數滿足要求.RSRR算法的其余參數分別設置為N=50，L=2，采樣點根據二維Chebyshev插值規(guī)則分布于圍道內，計算該模型的固有頻率和振型.將特征值計算結果與文獻[25]相同算例結果進行對比，見表1.可以看出，本文RSRR方法的計算結果和文獻[23]中的結果是一致的，但是本文的精度較高.該算例共計耗時41.4min，占用內存344.3MB.

表1 關注頻段內的全部固有頻率和誤差Table 1 Natural frequencies and error in the spectrum in attention

圖5 d取不同值時特征對誤差圖Fig.5 Relative error of eigenpairs to di ff erent d

為驗證 d的估計方法針對 RSRR算法的合理性，保持其他參數不變，對比d=10,16,20,26,30時各特征對的誤差，見圖5.可以看出當d=20時，特征對的精度正好達到最優(yōu)，即本文提出的確定d取值的方法針對RSRR算法是合理的.

4.3 大規(guī)模工程應用算例

將本文提出的 RSRR非線性特征值解法應用于實際工程中大規(guī)模吸聲結構的復模態(tài)分析.考慮火箭整流罩內聲腔模型，其半剖邊界元模型見圖6.該模型軸向最大尺寸為2.29m，徑向最大尺寸為1.324m，將其離散為64162個6節(jié)點非協(xié)調三角形單元，共計384972自由度，全部表面設定為阻抗邊界，阻抗模型與算例4.2相同.

圖6 整流罩邊界元模型Fig.6 BEM model of the fairing

本算例欲求矩形區(qū)域(200+0i,500+100i)Hz內的所有特征解，取C為該矩形.考慮阻抗邊界，給定邊界元核函數t(z)=Z·z·eizR，插值點取為圍道上的高斯積分點，d依次取10,16,20,26,30，根據式(33)計算t(z)和T(z)在不同插值項數下的Cauchy插值誤差，見圖7，T(z)的近似方法同4.2節(jié).可以看出，t(z)=f(z)eizR的Cauchy插值誤差能夠反映出對T(z)每個元素進行Cauchy插值的最大誤差，且當d取26時，矩陣近似精度已滿足要求.本算例取d=30.RSRR算法的其余參數分別設置為N=50，L=4，采樣點根據二維Chebyshev插值規(guī)則分布于關注頻段內，計算該模型的固有頻率和振型.

圖7 基于Cauchy插值的BEM矩陣和核函數最大誤差曲線Fig.7 Max relative error of BEM matrix and Kernel function based on the Cauchy interpolation

該算例耗時122.3h，占用內存21.023GB.共計算出51組特征對，其相對誤差分布見圖8，前四階固有頻率以及對應的聲壓模態(tài)見圖9.該算例顯示了RSRR算法對大規(guī)模問題仍然保持了較高的計算精度.

圖8 特征對誤差圖Fig.8 Relative error of eigenpairs

圖9 前四階聲壓模態(tài)Fig.9 The fi rst four modal shapes of the sound pressure

5 結論

本文提出了一種解決大規(guī)模邊界元特征值問題的數值方法RSRR，為工程中的大規(guī)模邊界元模態(tài)分析提供了有效計算手段.

(1)基于解析矩陣函數的Keldysh定理，提出了一種簡單易行、可靠性高的特征向量搜索空間構造方法，該方法直接采用一系列邊界元頻域解的一組基空間構造特征空間，無需對邊界元系數矩陣取高階矩，解決了傳統(tǒng)圍道積分方法存在的特征空間易出現(xiàn)病態(tài)甚至構造失敗的問題.

(2)基于解析函數的Cauchy積分公式，構造了復平面內邊界元系數矩陣的插值方法，以及縮減特征值問題系數矩陣的快速計算方法，還根據核函數與系數矩陣的關系給出了插值項數的估計策略.該方法在降低RSRR方法求解邊界元特征值問題計算量的同時，也將其擴展到了工程中常見的邊界元復特征值分析中.

(3)將所提出的RSRR方法與聲學快速邊界元結合，在普通工作站上實現(xiàn)了復雜邊界條件下自由度接近40萬的邊界元聲模態(tài)分析，顯示了本文方法的大規(guī)模計算分析能力.

1姚振漢,王海濤.邊界元法.北京：高等教育出版社,2010(Yao Zhenghan，Wang Haitao.Boundary Element Method.Beijing:Higher Education Press,2010(in Chinese))

2高效偉,劉健,彭海峰.集成單元邊界元法及其在主動冷卻熱防護系統(tǒng)分析中的應用.力學學報,2016,48(4):994-1003(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))

3申建偉,孫中奎,詹世革等.第9屆全國動力學與控制青年學者學術研討會報告綜述.力學學報,2015,47(6):1079-1083(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))

4榮俊杰,校金友,文立華.彈性動力學高階核無關快速多極邊界元法.力學學報,2014，46(5):776-785(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))

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

6 Tisseur F,Meerbergen K.The quadratic eigenvalue problem.SIAM Review,2001,43(2):235-286

7 Mehrmann V,Schr¨oder C.Nonlinear eigenvalue and frequency response problems in industrial practice.Journal of Mathematics in Industry,2011,1(1):7

8 E ff enberger C.Robust solution methods for nonlinear eigenvalue problems.[PhD Thesis].′Ecole polytechnique f′ed′erale de Lausanne,2013

9 Van Beeumen,Rational Krylov methods for nonlinear eigenvalue problems[PhD Thesis],KU Leuven,2015

10 Ali A,Rajakumar C,Yunus SM.Advances in acoustic eigenvalue analysis using boundary element method.Computers&Structures,1995,56(5):837-847

11 Bezine G.A mixed boundary integral— fi nite element approach to plate vibration problems.Mechanics Research Communications,1980,7(3):141-150

12 Nardini D,Brebbia CA.A new approach to free vibration analysis using boundary elements.Applied Mathematical Modelling,1983,7(3):157-162

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

14 Kamiya N,Andoh E.Standard eigenvalue analysis by boundaryelement method.International Journal for Numerical Methods in Biomedical Engineering,1993,9(6):489-495

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

16 Beyn WJ.An integral method for solving nonlinear eigenvalue problems.Linear Algebra and Its Applications,2012,436(10):3839-3863

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

18 Yokota S,Sakurai T.A projection method for nonlinear eigenvalue problems using contour integrals.JSIAM Letters,2013,5:41-44

19 Zheng CJ,Chen HB,Gao HF,et al.Is the Burton–Miller formulation really free of fi ctitious eigenfrequencies?Engineering Analysis with Boundary Elements,2015,59:43-51

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

21 LeblancA,LavieA.Solvingacousticnonlineareigenvalueproblems with a contour integral method.Engineering Analysis with Boundary Elements,2013,37(1):162-166

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

23 Xiao JY,Zhou H,Zhang CZ,et al.Solving large-scale fi nite element nonlinear eigenvalue problems by resolvent sampling based Rayleigh-Ritz method.Computational Mechanics,2016:1-18

24 Mehrmann V,Voss H.Nonlinear eigenvalue problems:A challenge for modern eigenvalue methods.GAMM-Mitteilungen,2004,27(2):121-152

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

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

27屈伸,陳浩然.敷設多孔吸聲材料聲腔特征值分析的徑向積分邊界元法.計算力學學報,2015,32(1):123-128(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))

28陳文炯,劉書田.周期性吸聲多孔材料微結構優(yōu)化設計.計算力學學報,2013,30(1):45-50(Chen Wenjiong,Liu Shutian.Optimizing design of micro-structural con fi gurations of periodic porous soundabsorbing materials.Chinese Journal of Computational Mechanics,2013,30(1):45-50(in Chinese))

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

30 Engquist B,Ying Lexing.Fast directional multilevel algorithms for oscillatory kernels.SIAM Journal on Scienti fi c Computing,2007,29(4):1710-1737

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

32 Gohberg I,Rodman L.Analytic matrix functions with prescribed local data.Journal d’Analyse Math′ematique,1981,40(1):90-128

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

AN EFFICIENT NUMERICAL METHOD FOR LARGE-SCALE MODAL ANALYSIS USING BOUNDARY ELEMENT METHOD1)

Wang Junpeng Xiao Jinyou2)Wen Lihua
(College of Astronautics，Northwestern Polytechnical University，Xi’an 710072，China)

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 fi elds using the BEM becomes an very important but challenging task;so far there are no e ff ective 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.

BEM,nonlinear eigenvalue problems,Rayleigh-Ritz projection,Cauchy integral formula,modal analysis

O302,TB122

A

10.6052/0459-1879-17-040

2017–02–15收稿，2017–07–17 錄用,2017–07–20 網絡版發(fā)表.

1)國家自然科學基金(11102154，11472217)和中央高?；究蒲袠I(yè)務費專項資金資助項目.

2)校金友，教授，主要研究方向：計算力學，邊界元法，結構振動與噪聲.E-mail:xiaojy@nwpu.edu.cn

王俊鵬,校金友,文立華.大規(guī)模邊界元模態(tài)分析的高效數值方法.力學學報,2017,49(5):1070-1080

Wang Junpeng,Xiao Jinyou,Wen Lihua.An efficient numerical method for large-scale modal analysis using boundary element method.Chinese Journal of Theoretical and Applied Mechanics,2017,49(5):1070-1080