李穎潔，王 燾，校金友，文立華

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

隨著噪聲污染受到人們越來越多的關(guān)注，在航空航天、汽車制造等領(lǐng)域中，為了降低內(nèi)部環(huán)境的噪聲，必須考慮相應(yīng)的降噪設(shè)計。特別是在航天領(lǐng)域，航天器在發(fā)射起飛階段通常要經(jīng)受嚴重的氣動力。氣動噪聲等惡劣的環(huán)境會導致航天器內(nèi)電子設(shè)備、有效載荷以及整體結(jié)構(gòu)的失效和破壞[1]。

工程中通常使用吸聲材料或吸聲結(jié)構(gòu)來吸收入射聲波的能量，從而達到降低噪聲的目的[2]。使用多孔材料降低噪聲效果顯著，已經(jīng)受到越來越多的關(guān)注[3–5]。目前，工程上一般將多孔材料作為常值阻抗邊界進行處理，然而由于多孔材料的結(jié)構(gòu)復雜，使得計算聲場時需要考慮由多孔材料引入的復雜頻變阻抗邊界，特別是在低頻范圍，阻抗頻變特性非常顯著，此時將阻抗假設(shè)為常值的方法是不準確的[6]。試驗證明許多吸聲材料的阻抗都有著明顯的頻變特性[7]。

通過模態(tài)分析計算聲模態(tài)和頻率，是指導吸聲結(jié)構(gòu)設(shè)計、考核吸聲性能的重要途徑。對聲阻抗為常數(shù)的情況，聲模態(tài)分析是求解常規(guī)二次特征值問題，有限元方法多采用LANCZOS迭代或子空間迭代方法，計算效率高，且可用于大規(guī)模工程問題的計算分析。然而當考慮阻抗的頻變性能時，聲模態(tài)分析變?yōu)榍蠼庖话愕姆蔷€性特征值問題，現(xiàn)有主流商用有限元分析軟件（如ANSYS，NASTRAN）還不具備這項功能。

事實上，由于工程結(jié)構(gòu)減振降噪技術(shù)的發(fā)展需求，非線性特征值問題的高效數(shù)值解法已成為近年來的一個國際研究熱點[8–12]。經(jīng)典的多項式近似法[8]以及改進方法[9]只能不斷近似逼近原問題的系數(shù)矩陣；牛頓迭代法[10]以及相應(yīng)的改進方法[11]都需要反復求解線性方程組，對于大規(guī)模問題，計算效率低。圍道積分法（Contour Integral Method）是近年來發(fā)展起來的一種特征值問題數(shù)值解法。與前面所提的迭代方法不同，它可以一次性計算出復平面上一個閉曲線（圍道）內(nèi)的所有特征值，并且便于并行化[12–13]，因而受到廣泛關(guān)注。但是，該方法計算結(jié)果受圍道積分參數(shù)影響很大，容易出現(xiàn)虛假特征根。

本文作者最近在圍道積分法的基礎(chǔ)上提出了一種新型特征值問題的數(shù)值解法，記為RSRR（Resolvent Sampling based Rayleigh-Ritz projection method）。RSRR不僅可以用于標準特征值和廣義特征值問題，也可用于一般的非線性特征值問題。它采用采樣法構(gòu)造特征空間，解決了圍道積分法中矩積分法構(gòu)造特征空間時高階矩失效的問題[14]，并且構(gòu)造特征空間的范圍比矩積分法的范圍更大，使特征解更精確，穩(wěn)定性更好。同時，它還保留了圍道積分法適于并行計算的優(yōu)勢[15]。

本文將RSRR法應(yīng)用于含復雜阻抗邊界的聲場模態(tài)分析?；谟邢拊治鲕浖嗀NSYS，實現(xiàn)了RSRR，從而可將其用于一般復雜工程問題的建模、分析和后處理；通過典型算例驗證了方法的正確性，并通過航天整流罩模型的大規(guī)模計算分析，驗證了對實際工程問題的求解效率。本文方法和阻抗表達式無關(guān)，是一種通用的聲模態(tài)工程分析方法。

1 頻變阻抗非線性特征值問題

1.1 吸聲材料的頻變阻抗模型

工程上常用的吸聲降噪的材料有蜂窩板、泡沫材料、玻璃纖維、玻璃棉、纖維毯等。不同材料的吸聲特性不同。進行聲學特性分析時，通常是將吸聲材料等效為聲場的阻抗邊界條件。目前工程應(yīng)用中普遍將由吸聲材料引起的阻抗等效為常值進行計算，而實際上，吸聲材料的阻抗是隨頻率變化而改變的。在低頻范圍，受頻率影響很大。文獻[16]通過試驗得到玻璃纖維的阻抗-頻率曲線如圖1所示?？梢钥闯?，將阻抗等效為常值在低頻范圍會給分析帶來很大誤差。

圖1 玻璃纖維的聲阻抗

下面分別介紹兩種吸聲材料的頻變聲阻抗模型。在第四節(jié)的算例中將考慮這兩個模型。

Kelvin-Voigt模型：

Kelvin-Voigt模型采用一個簡單的一維模型模擬聲場表面多孔材料的聲阻抗。阻抗隨頻率變化的表達式為

其中參數(shù)kI和dI分別表示吸聲材料的彈性和粘性特性，ω為角頻率，j為虛數(shù)單位。這兩個參數(shù)通?？梢赃M行實驗或通過Biot-Allard理論得到。例如，對于圖1所示的玻璃纖維的試驗結(jié)果，可確定參數(shù)dI=50 Pa?s/m ，kI=5×106Pa?/m 。

Delany-Bazley模型：

Delany和Bazley對多孔纖維材料的阻抗特性進行大量的測量和研究，得到表面阻抗的表達式為

其中d是材料厚度，?是孔隙率，Zc為材料特征阻抗，k為復波數(shù)，δ=ρ0ω/2πσ，σ是多孔材料的電阻率，ρ0是密度，c是聲速。

Delany-Bazley模型在多孔吸聲材料建模中應(yīng)用很廣。文獻[7]中則通過實驗驗證了洋麻、大麻、椰肉、稻草等多種纖維吸聲材料吸聲性能的頻變特性，且采用Delany-Bazley模型擬合了各種纖維材料的頻變阻抗。

1.2 聲學有限元模態(tài)分析

考慮三維聲腔V，聲腔內(nèi)介質(zhì)為理想的非粘性介質(zhì)。頻域內(nèi)現(xiàn)行聲學Helmholtz方程為[17]。

其中p(r)為聲壓函數(shù)，q(r)為體積速度函數(shù)，ρ0為聲場密度，ω為角頻率。聲場邊界為Ω。考慮一般情況，邊界Ω由Dirichlet邊界Ωp、Neumann邊界ΩV和Robin邊界ΩZ三部分組成，即Ω=ΩP?ΩV?ΩZ，相應(yīng)的三種邊界條件分別為

采用加權(quán)余量法，獲得Helmholtz方程的等效積分弱形式，再將聲場V進行單元離散，對場變量進行插值近似，則可得到純聲學有限元方程。

其中Ka為聲剛度矩陣，Ma為聲質(zhì)量矩陣，Qi為聲學激勵向量，Vni為輸入速度向量，Pi為輸入聲壓向量。Ca為聲阻尼矩陣，與邊界聲阻抗有關(guān)，i行j列的矩陣元素計算式為

其中Ωz為阻抗邊界。Ni、Nj分別為有限元形函數(shù)矩陣的第i行元素和第j列元素。模態(tài)分析時，式(9)中的激勵為零，形成特征值問題

式中λ表示角頻率ω，v表示聲壓向量{pi}

由阻尼矩陣的計算式可以看出，當阻抗為常值時，阻尼陣Ca則為常矩陣。聲學特征值問題為常規(guī)二次特征值問題。而當聲阻抗為頻率的函數(shù)(ω)時，聲學特征值問題的求解則轉(zhuǎn)化為非線性特征值問題。雖然這類特征值問題在吸聲降噪設(shè)計分析以及其他工程領(lǐng)域應(yīng)用很廣，但目前仍沒有十分有效的大規(guī)模數(shù)值解法，主流商業(yè)有限元軟件也不具備直接求解這類非線性特征值問題的模塊。

2 RSRR算法

基于預(yù)解矩陣采樣的瑞利里茲法(記為RSRR)是我們團隊最近發(fā)展的非線性特征值問題的通用解法[15]。本節(jié)將該方法應(yīng)用于特征值問題的數(shù)值求解，基于ANSYS聲場分析模塊，二次開發(fā)了復雜聲場模態(tài)分析功能。下面首先介紹RSRR算法的基本原理，然后討論在ANSYS軟件中的數(shù)值實現(xiàn)方法。本文的實現(xiàn)思路在其他聲學有限元分析軟件中同樣適用。

考慮非線性特征值問題T(λ)v=0，現(xiàn)欲求復平面內(nèi)由簡單閉曲線C圍成的區(qū)域里面的所有特征值及其對應(yīng)的特征向量。RSRR算法計算過程分為兩步：一是采用預(yù)解矩陣采樣法構(gòu)造包含所有特征向量的特征空間S，取S的一組正交基為Q；二是基于Q采用經(jīng)典Raylei-Ritz投影法，將原特征值問題轉(zhuǎn)化為小型縮減特征值問題TQ(λ)g=0，并進行求解，其中，為Q的列數(shù)?？s減特征值問題的維數(shù)一般為幾十至幾百，屬于小規(guī)模問題，現(xiàn)有許多方法都可以用于求解，如NEWTON法、修改的圍道積分法[15]等。縮減特征值問題與原問題的特征值相同，對應(yīng)的特征向量存在轉(zhuǎn)化關(guān)系v=Qg。

2.1 特征空間的構(gòu)造

RSRR算法中，特征空間的構(gòu)造采用預(yù)解矩陣采樣法，步驟如算法1所示。

算法1：預(yù)解矩陣采樣算法

在采樣點位置，采樣點數(shù)目N，以及隨機向量數(shù)目L選取合適的情況下，構(gòu)造的空間S近似包含圍道C內(nèi)所有特征向量，原因如下。

根據(jù)Keldysh理論[15]，矩陣T-1(z)可由圍道C內(nèi)T(z)的特征值和特征向量表示

其中nc為互不相同的特征根數(shù)目，ηk表示非零空間的維度，μl,k表示λk的第l階局部重數(shù)。和分別是T(λk)的右、左正規(guī)化廣義特征向量，RC(z)表示圍道C之外的所有特征值的貢獻。由式可知，矩陣T(z)的特征值即為T-1(z)的極點。

構(gòu)造特征空間的目的是，從式(13)中提取出特征向量VC。為此，首先將式寫為如下矩陣形式

函數(shù)矩陣RC(z)代表圍道C之外所有特征值的貢獻，在圍道C內(nèi)部是解析函數(shù)。因此，可在一組給定的基函數(shù)gj(z)(j=1,…J)下展開成如下形式

圖2 采樣圍道示意圖

式中Rj是展開系數(shù)矩陣。常用的基函數(shù)有Lagrange多項式、Chebyshev多項式等；根據(jù)多項式逼近理論，對于解析函數(shù)，當插值點和基函數(shù)選取合適時，能獲得指數(shù)收斂速度。而對于本文的特征空間構(gòu)造方法，并不需要明確計算基函數(shù)，基函數(shù)的影響體現(xiàn)在采樣點（即插值點）的選取上。比如，當C為一個實區(qū)間時，采樣點可取為該區(qū)間上的Chebyshev點，相當于基函數(shù)gj(z)為Chebyshev多項式，式(16)的逼近精度很高。

隨機矩陣U∈?n×L的作用是對T-1(z)進行探測，獲取特征向量矩陣VC的信息。由T(z)矩陣Jordan標準型的矩陣結(jié)構(gòu)可知[14]，U的列數(shù)L至少要等于矩陣T(z)在圍道C內(nèi)最大的特征根重數(shù)，即

將式(15)代入式(14)之后，T-1(z)U可表示為

而矩陣Z(z)的表達式為

其中Dj(z)為L維對角陣，對角線元素的值都是gj(z)。

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

為有效提取出?空間，本文選取N和采樣點，形成采樣解空間S

在瑞利里茲投影中需要S的正交基Q，可以通過S的截斷奇異值分解（SVD）計算得到，即S≈QΣVH。由于S的列數(shù)N?L通常很小，并且獨立于問題的尺寸n。所以，SVD分解的計算量隨n線性變化。令ks代表S經(jīng)過SVD分解后的數(shù)值秩，要得到圍道C內(nèi)的所有特征值，則要滿足ks≥rank(VC)。最后需要指出，當采樣點距離特征值很近時，T(zi)-1U的數(shù)值會很大，為保證特征空間構(gòu)造過程的穩(wěn)定性，在SVD分解之前應(yīng)先對矩陣S的列向量進行歸一化處理。

2.2 算法的實現(xiàn)

采用RSRR算法求解非線性特征值問題的步驟如下。

算法 2：RSRR算法

本文基于ANSYS的參數(shù)設(shè)計語言（APDL）實現(xiàn)了以上RSRR算法，最終以二次開發(fā)的形式在ANSYS軟件中實現(xiàn)了考慮頻變阻抗邊界的聲場模態(tài)分析功能。下面詳述實現(xiàn)中的關(guān)鍵技術(shù)問題和處理方法。

物理建模在ANSYS前處理環(huán)境下進行。具體包括建立聲場模型，劃分網(wǎng)格和施加聲場邊界條件。由于ANSYS中對聲場無法施加頻變阻抗邊界條件，故可將頻變阻抗用任一常值阻抗代替，而在后面的二次開發(fā)程序中將其替換為頻變表達式。ANSYS求解模塊設(shè)置求解類型為模態(tài)分析，在求解有限元方程組之前中止，將質(zhì)量矩陣Ma、剛度矩陣Ka和阻尼矩陣Ca導出到數(shù)據(jù)文件中。

RSRR算法參數(shù)的確定。實際應(yīng)用中通常是欲求某個區(qū)間或復平面上某個區(qū)域內(nèi)的特征值，據(jù)此可確定特征值搜索范圍，邊界圍道C取橢圓或矩形較為方便。確定采樣點數(shù)目N，需要對待求特征值的數(shù)目進行預(yù)估，由于此處對預(yù)估精度要求不高，可以通過常值阻抗時的特征值數(shù)目來估計。采樣點可取為圍道C上積分點或圍道內(nèi)Chebyshev點。L不小于最大特征根重數(shù)，這可以根據(jù)模型的對稱性進行估計。計算經(jīng)驗表明，為保證特征根的計算精度，參數(shù)N?L的選取應(yīng)大于特征值數(shù)目的2～3倍。生成隨機矩陣U∈?n×L之后，對其進行歸一化有利于方法的穩(wěn)定性。

構(gòu)造特征空間Q時，首先將頻變式(12)阻抗模型代入式(10)，從文件讀入Ma、Ka和Ca并按式通過APDL Math重新計算T矩陣。調(diào)用ANSYS直接求解器，在每個采樣點上計算T(zi)-1U。將求解結(jié)果組裝成矩陣S。ANSYS中多波前法求解器的計算效率高、穩(wěn)定性好，在ANSYS中進行這部分計算保證了本文特征值解法較高的總體計算效率。對S的截斷奇異值分解可將MATLAB中的奇異值分解函數(shù)編譯成.exe文件，在ANSYS中調(diào)用。分解之前對其列向量歸一化，對保證方法的穩(wěn)定性非常重要。奇異值截斷閾值一般取為δ=10-14。

聚縮后小型特征值問題系數(shù)矩陣TQ的計算，主要是矩陣乘法運算，在ANSYS中進行。本文采用改進的Block-SS算法[15]求解小規(guī)模非線性特征值問題，其中的圍道仍取為RSRR中的圍道C。得到縮減特征值問題的特征對(λj,gj)。這個求解過程可通過MATLAB語言編程實現(xiàn)，并編譯成.exe文件在ANSYS中調(diào)用。

3 算例

提供三個算例驗證本文方法和ANSYS實現(xiàn)的正確性和計算效率。前兩個算例中聲場邊界阻抗模型分別為第2.1節(jié)中的Kelvin-Voigt模型和Delany-Bazley模型，目的是和現(xiàn)有文獻結(jié)果對比，驗證本文方法的正確性。第三個算例對鋪設(shè)吸聲材料的火箭整流罩模型進行聲模態(tài)分析，考察本文方法對于大規(guī)模問題的適用性和計算效率。所有算例均在一臺32核（Xeon E5-2630 v3，2.40 GHz）64 GB內(nèi)存的工作站上進行。線性方程組的求解均采用ANSYS軟件中frontal直接解法，采用16核并行運算。

3.1 算例1

考慮文獻[19]中尺寸為1 m×0.1 m×1 m的立方體聲場。在z=1 m處根據(jù)式(1)設(shè)置阻抗邊界條件，其中dI=50 Pa?s/m ，kI=5×106Pa?/m 。其余各面設(shè)置為剛性邊界。聲場介質(zhì)為空氣。本文空氣密度均取為ρ0=1.2 kg/m3，聲速均取為c=340 m/s。在有限元模型中，將聲場劃分為12 500個8節(jié)點6面體聲場單元。

為了方便與文獻[19]的解析解對比，本文計算50 Hz到350 Hz范圍內(nèi)的5個特征解及誤差。RSRR方法采樣點取在橢圓圍道上，橢圓中心位于（200，0），橢圓長軸a=150，橢圓短軸b=0.1a。根據(jù)3.2節(jié)關(guān)于參數(shù)N和L取法的討論，采樣點N取10，L取2。

表1 常值阻抗與頻變阻抗計算結(jié)果對比

RSRR方法求得的特征值見表1中第2列。此表第3列為文獻[19]中提供的解析解，與本文結(jié)果基本相同，證明了本文數(shù)值計算結(jié)果的正確性。為考察特征值和特征向量的綜合誤差，還計算了5個相應(yīng)特征對的相對誤差||T(λ)v||2/||v||2，均在10-10以下。

進一步研究了考慮阻抗的頻變特性對模態(tài)分析結(jié)果的影響。對文獻[19]中的模型，在z=1 m的表面處施加Zs=104Pa s/m的常值阻抗進行計算，可得相應(yīng)的特征頻率如表1中第4列所示。與考慮頻變阻抗的結(jié)果對比可知，將阻抗用常值代替會帶來模態(tài)頻率上的較大差異，這是因為500 Hz以下該吸聲材料Kelvin-Voigt模型的阻抗變化很大，如圖1所示。因此，對于阻抗頻變較為嚴重且對特征值精度要求較高的情況，必須考慮邊界阻抗的頻變特性。

3.2 算例2

考慮文獻[6]中的0.5 m×0.4 m×0.3 m立方體聲場模型，在x=0.5 m處施加阻抗邊界條件，阻抗模型為式(2)的Delany-Bazley模型，其中有關(guān)參數(shù)取值為d=20 cm 、σ=2×104Ns/m4、?=1 和ρ0=1.2 kg/m3。其余邊界為剛性邊界。聲場介質(zhì)為空氣，劃分7 500個8節(jié)點6面體聲場單元。

通過RSRR方法計算300 Hz～750 Hz頻率范圍內(nèi)的八階特征值以及相應(yīng)的特征對誤差。在矩形區(qū)域[200+0j；500+200j]邊界取點。取N=30 ，L=1。計算所得結(jié)果與文獻[6]進行對比。

由表2結(jié)果可見，對于較為復雜的阻抗模型，文獻[6]中采用特殊積分法對阻抗模型進行2次近似時會帶來一定的誤差。而采用RSRR法求解復雜阻抗邊界條件的聲學特征值時，不需要對頻變阻抗模型進行近似，所以更接近真實特征值，特征對的求解誤差也均在10-11以下。因此RSRR方法求解此類復雜頻變阻抗有限元聲學特征值問題具有優(yōu)越性。

表2 Delany-Bazley阻抗模型的聲學特征值求解

3.3 算例3

火箭整流罩結(jié)構(gòu)的作用是保護航天器的有效載荷、降低噪聲和振動的影響，為此通常在前錐和圓筒段鋪設(shè)吸聲材料以達到減振降噪的效果。本算例考慮圖3所示的內(nèi)部有衛(wèi)星的整流罩內(nèi)聲場。整流罩軸向最大尺寸為2.288 m，徑向最大尺寸為1.325 m。內(nèi)部衛(wèi)星軸向最大尺寸為0.95 m，徑向最大尺寸為0.75 m。聲場介質(zhì)為空氣，在整流罩內(nèi)表面深色區(qū)域鋪設(shè)多孔纖維吸聲材料，阻抗模型采用Delany-Bazley模型，有關(guān)參數(shù)取值與算例2相同。其余邊界為剛性邊界。由于模型結(jié)構(gòu)較為復雜，劃分單元時混合使用4面體和6面體單元，有限元模型規(guī)模為1 129 819自由度。

圖3 考慮衛(wèi)星的整流罩內(nèi)部結(jié)構(gòu)示意圖

采用RSRR方法計算復平面矩形區(qū)域[200+0j，500+100j]內(nèi)的所有特征值和特征向量。采樣點取為該區(qū)域內(nèi)的Chebyshev點，設(shè)N=102、L=2；采樣點分布如圖4所示。

圖4 采用RSRR計算有限元離散的鋪設(shè)多孔材料的整流罩內(nèi)聲場特征頻率（o表示采樣點，+代表求得的特征值）

RSRR方法共求得50個特征對，整個求解過程耗時約9.4 h，占用內(nèi)存26.7 GB，特征值分布如圖4所示，特征對的相對誤差如圖5所示。

圖5 整流罩特征根的求解誤差

圖6給出了計算所得前4階聲壓模態(tài)及相應(yīng)的特征頻率。

該算例表明，采用RSRR方法求解含復雜阻抗邊界的大規(guī)模聲模態(tài)分析問題，仍具有較高的計算效率和精度。

圖6 整流罩前4階模態(tài)圖

4 結(jié)語

本文將新型非線性特征值問題數(shù)值解法RSRR應(yīng)用于有限元聲模態(tài)分析中，解決了現(xiàn)有商用有限元軟件不能進行復雜頻變阻抗邊界的聲模態(tài)分析問題。

基于RSRR算法，在ANSYS中采用樣本法構(gòu)造特征向量搜索空間，采用此特征空間對原大規(guī)模非線性特征值問題進行瑞利里茲投影，使原問題轉(zhuǎn)化為小型非線性特征值問題，采用改進的Block-SS算法進行求解。采用采樣法構(gòu)造特征空間，既簡化了計算過程，又使得結(jié)果更精確。本文的方法可以穩(wěn)定高效地求解復雜阻抗邊界的聲學特征值，并且可用于大規(guī)模工程問題的計算。最后通過算例驗證了此方法的準確性和通用性。

[1]陳釗.火箭整流罩內(nèi)聲場分析及降噪技術(shù)研究[D].哈爾濱：哈爾濱工業(yè)大學，2015.

[2]陳克安.有源噪聲控制[M].北京：國防工業(yè)出版社，2014.

[3]裴春明，周兵，李登科，等.多孔材料和微穿孔板復合吸聲結(jié)構(gòu)研究[J].噪聲與振動控制，2015，35(5)：35-38.

[4]王營，趙武，黃丹.多孔材料聲學模型及其應(yīng)用[J].材料導報，2015，29(5)：145-149.

[5]屈伸，陳浩然.敷設(shè)多孔吸聲材料聲腔特征值分析的徑向積分邊界元法[J].計算力學學報，2015(1)：123-128.

[6]LEBLANC A,LAVIE A.Numericalanalysisof eigenproblem for cavities by a particular integral method with a low frequency approximation of surface admittance[J].Journal of the Acoustical Society of America,2012,131(5):3876-3882.

[7]BERARDIU,IANNACEG.Predictingthesound absorption of natural materials:Best-fit inverse laws for the acoustic impedance and the propagation constant[J].AppliedAcoustics,2017,115:131-138.

[8]LARBI W,DEü J F,OHAYON R.A new finite element formulation for internal acoustic problems with dissipative walls[J].International Journal for Numerical Methods in Engineering,2010,68(3):381-399.

[9]EFFENBERGER C.Robust successive computation of eigenpairs for nonlinear eigenvalue problems[J].SIAM J MatrixAnalAppl,2013,34(3):1231-1256.

[10]NAKA Y,OBERAI A A,SHINNCUNNINGHAM B G.Acoustic eigenvalues of rectangular rooms with arbitrary wall impedances using the interval Newton/generalized bisection method[J].Journal of the Acoustical Society ofAmerica,2005,118(6):3662-3671.

[11]VAN BEEUMEN R,MEERBERGEN K,MICHIELS W(2015)Compact rational Krylov methods for nonlinear eigenvalue problems[J].SIAM J Matrix Anal Appl,36(2):820-838.

[12]LEBLANC A,LAVIE A.Solving acoustic nonlinear eigenvalue problems with a contour integral method[J].Engineering Analysis with Boundary Elements,2013,37(1):162-166.

[13]XIAO J,ZHOU H,ZHANG C,et al.Solving large-scale finite element nonlinear eigenvalue problems by resolvent sampling based Rayleigh- Ritz method[J].Computational Mechanics,2016:1-18.

[14]YOKOTA S,SAKURAI T.A projection method for nonlinear eigenvalue problems using contour integrals[J].Jsiam Letters,2013,5:41-44.

[15]XIAO J,MENG S,ZHANG C,et al.Resolvent sampling based Rayleigh-Ritz method for large-scale nonlinear eigenvalue problems[J].Computer Methods in Applied Mechanics&Engineering,2016,310:33-57.

[16]LADEVEZE P,ZIENKIEWICZ O.New advances in computational structural mechanics[M].1992.

[17]李增剛，詹福良.Virtual.Lab Acoustics聲學仿真計算高級應(yīng)用實例[M].北京：國防工業(yè)出版社，2010.

[18]XIAO J,ZHANG C,HUANG T M,et al.Solving largescale nonlinear eigenvalue problems by rational interpolation approach and resolvent sampling based Rayleigh-Ritz method[J].International Journal for Numerical Methods in Engineering,2016.

[19]BERMúDEZA,RODRíGUEZR.Modellingand numericalsolution ofelastoacoustic vibrationswith interface damping[J].InternationalJournalfor Numerical Methods in Engineering,1999,46(10):1763-1779.