復雜阻抗吸聲結構的大規模有限元聲模態分析方法

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

（西北工業大學 航天學院，西安 710072）

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

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

通過模態分析計算聲模態和頻率，是指導吸聲結構設計、考核吸聲性能的重要途徑。對聲阻抗為常數的情況，聲模態分析是求解常規二次特征值問題，有限元方法多采用LANCZOS迭代或子空間迭代方法，計算效率高，且可用于大規模工程問題的計算分析。然而當考慮阻抗的頻變性能時，聲模態分析變為求解一般的非線性特征值問題，現有主流商用有限元分析軟件（如ANSYS，NASTRAN）還不具備這項功能。

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

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

本文將RSRR法應用于含復雜阻抗邊界的聲場模態分析。基于有限元分析軟件ANSYS，實現了RSRR，從而可將其用于一般復雜工程問題的建模、分析和后處理；通過典型算例驗證了方法的正確性，并通過航天整流罩模型的大規模計算分析，驗證了對實際工程問題的求解效率。本文方法和阻抗表達式無關，是一種通用的聲模態工程分析方法。

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

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

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

圖1 玻璃纖維的聲阻抗

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

Kelvin-Voigt模型：

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

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

Delany-Bazley模型：

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

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

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

1.2 聲學有限元模態分析

考慮三維聲腔V，聲腔內介質為理想的非粘性介質。頻域內現行聲學Helmholtz方程為[17]。

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

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

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

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

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

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

2 RSRR算法

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

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

2.1 特征空間的構造

RSRR算法中，特征空間的構造采用預解矩陣采樣法，步驟如算法1所示。

算法1：預解矩陣采樣算法

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

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

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

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

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

圖2 采樣圍道示意圖

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

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

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

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

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

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

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

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

2.2 算法的實現

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

算法 2：RSRR算法

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

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

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

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

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

3 算例

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

3.1 算例1

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

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

表1 常值阻抗與頻變阻抗計算結果對比

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

進一步研究了考慮阻抗的頻變特性對模態分析結果的影響。對文獻[19]中的模型，在z=1 m的表面處施加Zs=104Pa s/m的常值阻抗進行計算，可得相應的特征頻率如表1中第4列所示。與考慮頻變阻抗的結果對比可知，將阻抗用常值代替會帶來模態頻率上的較大差異，這是因為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模型，其中有關參數取值為d=20 cm 、σ=2×104Ns/m4、?=1 和ρ0=1.2 kg/m3。其余邊界為剛性邊界。聲場介質為空氣，劃分7 500個8節點6面體聲場單元。

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

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

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

3.3 算例3

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

圖3 考慮衛星的整流罩內部結構示意圖

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

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

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

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

圖6給出了計算所得前4階聲壓模態及相應的特征頻率。

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

圖6 整流罩前4階模態圖

4 結語

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

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

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

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

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

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

[5]屈伸，陳浩然.敷設多孔吸聲材料聲腔特征值分析的徑向積分邊界元法[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聲學仿真計算高級應用實例[M].北京：國防工業出版社，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.