頻變阻抗邊界聲空間特征頻率靈敏度分析的有限元法

劉 強，梁夢輝，鄭昌軍，畢傳興

(合肥工業大學 噪聲振動工程研究所，合肥 230009)

隨著我國工業化進程的加速，噪聲問題正日益受到重視。聲學優化設計為結構的低噪聲設計提供了一種有效思路，得到了越來越多的關注[1-3]。通常把聲學優化設計按照設計變量的類型分為3個層次，即尺寸優化、形狀優化和拓撲優化。靈敏度分析為基于梯度的優化方法提供了優化方向和量化依據，在聲學優化設計過程中扮演著非常重要的角色[4-5]。

空腔噪聲不僅會影響汽車等交通工具的乘坐舒適性，而且還可能會致使航天器內部的電子設備和結構發生失效和破壞。為了降低空腔噪聲，工程中常采用在聲腔邊界鋪設多孔吸聲材料或吸聲結構等方法。多孔吸聲材料內部具有大量連通的孔隙，聲波進入材料內部傳播時能量會不斷損耗，從而起到吸聲作用。由于多孔吸聲材料成本低且降噪效果明顯，因此得到了廣泛應用。為了便于處理，聲學仿真分析中一般將多孔吸聲材料等效為聲場的阻抗邊界條件[6-7]。通過聲模態分析可以得到聲腔特征頻率和振型等信息，可以為吸聲性能的考核、吸聲結構的設計提供參考[8]。若將多孔吸聲材料簡化為常值阻抗邊界條件，聲模態分析求解的是常規二次特征值問題，可以采用子空間迭代或LANCZOS迭代等方法，計算效率高且適合大規模工程問題的求解。然而，由于多孔吸聲材料結構復雜，吸聲性能與材料的流阻率、孔隙率、厚度以及背板條件等諸多因素有關，因此通常具有頻變特性，尤其是在低頻段[9-10]。對于頻變阻抗邊界條件，聲模態分析求解的是典型的非線性特征值問題[11]，求解過程復雜且現有主流有限元分析軟件(如ANSYS、COMSOL Multiphysics等)尚不具有該功能。這也使在此基礎上進行聲模態的靈敏度分析變得更加困難。

對于非線性特征值問題的求解，近些年來基于圍道積分法發展起來一類新的數值解法[12-13]。該類方法不需迭代，能夠通過求解若干線性方程組將非線性特征值問題轉化成規模很小的廣義特征值問題，從而一次性地計算出選定復區間內的全部特征值和特征向量。該類方法的優勢是計算精度高，且適合大規模并行計算。在前期工作中，本課題組已基于圍道積分法發展了聲模態的邊界元分析方法[14]和聲振耦合模態的有限元-邊界元耦合分析方法[15]。然而，這些研究所求解的非線性特征值問題并不是由材料的頻變特性所產生的，而只是由邊界元系數矩陣的頻率相關性所導致的。

本文針對由頻變阻抗邊界條件所導致的非線性特征值問題，基于聲學有限元法構建了一種聲學特征頻率的靈敏度分析方法。首先，采用圍道積分法將非線性特征值問題轉換為規模很小的廣義特征值問題；然后，通過對廣義特征值問題的直接求導，構建聲學特征頻率的靈敏度分析方法。數值算例驗證了該方法的準確性和適用性。

1 聲學非線性特征值問題

當考慮非黏性理想聲介質時，空腔Ω內的時諧聲場的頻域控制微分方程為

?2p(x)+k2p(x)=0,x∈Ω

(1)

式中：?2為拉普拉斯算子；p(x)為點x處的聲壓；k=w/c=2πf/c為波數；w為角頻率；c為聲速；f為頻率。

為了降低空腔噪聲，工程中常采用在聲場邊界鋪設多孔吸聲材料等方法。常用的多孔吸聲材料有纖維板、泡沫鋁、玻璃棉、聚氨酯海綿、泡沫塑料等，不同的多孔材料具有不同的吸聲特性。多孔材料的吸聲性能與材料的流阻率、孔隙率、厚度以及背板條件等因素有關，并且通常具有頻變特性，特別是在低頻段。當進行聲學特性分析時，通常將多孔吸聲材料等效為聲場的阻抗邊界條件。目前常用的等效模型有Delany-Bazley模型和Miki模型等。由這些等效模型得到的聲阻抗均為頻率的函數，因此相應的阻抗邊界條件可以寫成

(2)

采用加權余量法獲得控制微分方程式(1)的等效積分弱形式，代入邊界條件，再將聲腔Ω進行單元離散，并對場變量進行插值近似，經整體組裝后就可以得到聲學有限元系統方程

(K+ikC-k2M)p=F

(3)

式中：K，C，M分別為聲剛度矩陣、聲阻尼矩陣和聲質量矩陣；p為節點聲壓向量；F為聲激勵向量。聲阻尼矩陣C與邊界S上的聲阻抗有關，K，M和C的單位矩陣具體形式分別為

(4a)

(4b)

(4c)

若令T(k)=K+ikC-k2M，則在聲模態分析時，當齊次系統方程

T(λ)ψ=0

(5)

具有非平凡解，即當ψ≠0時，λ為特征值，相應的ψ稱為特征向量。需注意本文中的特征值λ對應波數k，其與角頻率w和頻率f的關系為k=w/c=2πf/c。

由聲阻尼矩陣C的單元矩陣計算式(4c)可知：當假設吸聲材料的等效聲阻抗為常值時，聲阻尼矩陣C為常矩陣，此時式(5)所描述的聲學特征值問題為一個常規的二次特征值問題；而當考慮更接近實際的頻變聲阻抗時，例如采用Delany-Bazley模型獲得等效聲阻抗值，聲阻尼矩陣C則為頻率的非線性矩陣函數，這就導致式(5)所描述的聲學特征值問題成為一個典型的非線性特征值問題。相對于二次特征值問題，非線性特征值問題的求解要困難很多，目前仍缺少準確高效的通用數值解法，這也使在此基礎上進行聲學特征值的靈敏度分析變得更加困難。

2 非線性特征值的靈敏度分析

為了求解式(5)所描述的聲學非線性特征值的靈敏度，本章首先基于近一二十年來發展起來的圍道積分法將非線性特征值問題轉化為易于求解的小規模廣義特征值問題；然后，在此基礎上構建一種聲學特征值的靈敏度分析方法。該方法不需迭代，規避了傳統迭代方法在求解非線性特征值問題時可能遇到的解的發散問題和計算成本過高等不足[16]。

2.1 非線性特征值問題的線性化

為實現非線性特征值問題的轉化求解，我們首先在復平面內選取一個能夠包圍待求特征值區間的封閉曲線C，然后形成如下一類矩陣

(6)

通過矩陣M可以構造出塊對稱的Hankel矩陣H1和H2，即

(7)

(8)

由Asakura等的研究可知，矩陣H2相對于H1的廣義特征值等價于原始非線性特征值問題式(5)在封閉曲線C內的特征值。若用λj表示封閉曲線C內的第j個特征值，vj為特征值λj對應的特征向量(若無指出，以下特征向量均指右特征向量)，那么H2相對于H1的廣義特征值問題可以表示為

(H2-λjH1)vj=0

(9)

式(9)描述的廣義特征值問題可以很容易地轉換為標準線性特征值問題進行求解，在求得特征值和對應的特征向量后，式(5)的原始特征向量ψj可由特征向量vj通過

ψj=Svj

(10)

計算得到。這里的矩陣S=[S0,S1,…,SK-1]，其中的子矩陣S為

(11)

對比矩陣H1，H2和T的階數，可以發現轉換后的廣義特征值問題的規模要遠小于原非線性特征值問題的規模，因此式(9)描述的廣義特征值問題能夠快速地求解。另外，在計算矩陣M和S時，參數K和L的選取以及特征值數目m的確定可以參考Zheng等的研究。

2.2 單特征值的靈敏度分析

對于單特征值(即代數重數為1的特征值)，為了求其對任意設計變量的靈敏度值，我們將轉換后的小規模廣義特征值問題，即式(9)，直接對設計變量求導得到

λ′jH1vj=(H′2-λjH′1)vj+(H2-λjH1)v′j

(12)

式中，()′=?()/??為對設計變量?的導數。此處的設計變量?可以為形狀參數、材料屬性參數、拓撲參數等，本文側重于形狀參數的靈敏度分析。

由于式(12)中的λ′j和v′j皆未知，因此無法直接進行求解。為此在求解之前，對其增加約束條件，若采用最大歸一化方法將特征向量vj正則化后再對設計變量求導[17]，可以得到

?jv′j=0

(13)

式中，?j={0,…,σ,…,0}為KL階權重向量，σ為第e個位置的非零任意常數，e為原特征向量vj絕對值最大的元素的位置。

聯立式(12)和式(13)可以得到

(14)

可以證明，式(14)的系數矩陣為滿秩矩陣。因此，單特征值λj的靈敏度值λ′j可以通過求解式(14)得到。

2.3 重特征值的靈敏度分析

對于重特征值λ，若假設其代數重數為r(r>1)，那么該重特征值及其對應的右特征向量應滿足

H2X=H1XΛ

(15)

式中：X為λ對應的所有右特征向量組成的KL×r階矩陣；Λ=λI為r×r階對角矩陣，I為單位矩陣。顯然，X的列向量張成的向量空間中的任意向量都是λ對應的特征向量，但并不是所有特征向量都滿足隨設計變量連續變化的要求，因此，若仍采用式(14)計算重特征值λ對應的靈敏度值，將可能會導致計算結果的不準確。

(16)

(17)

WΓ=ΓΛ′

(18)

注意到式(14)和式(18)中包含Hankel矩陣H1和H2對設計變量的導數，根據式(7)和式(8)可知

(19)

(20)

由于隨機矩陣U和V與設計變量無關，因此矩陣M對設計變量的導數為

(21)

式中，Tn(z)=-T(z)-1T′(z)T(z)-1，T′(z)為有限元系數矩陣T(z)對設計變量的導數。當考慮的設計變量與頻率無關時，T′(z)為

T′(z)=K′+izC′-z2M′

(22)

此處的K′=?K/??，C′=?C/??，M′=?M/??，為了簡便一般可以通過差分方法計算得到。然而，差分方法的計算精度受差分步長影響，且由于需要多次計算系數矩陣也影響了計算效率。為了更準確、高效地計算這些導數矩陣，本文采用直接微分方法，直接對聲剛度矩陣、聲阻尼矩陣和聲質量矩陣求設計變量的導數，從而得到導數矩陣的解析表達式，類似過程可以參考Wang等[18]的研究中關于結構剛度矩陣對設計變量導數的計算過程，這里不再贅述。

注意到式(6)、式(11)和式(21)中需要計算一系列形如

(23)

的圍道積分。當關注的特征值區間為[λmin,λmax]時，若選取橢圓作為積分路徑，且橢圓的長軸和短軸分別與復平面的實軸和虛軸平行，并且采用NS點梯形公式進行數值積分，式(23)可以轉換為

(24)

式中：ρ=(λmax-λmin)/2為橢圓長半徑；γ=(λmax+λmin)/2為橢圓圓心；ζ為橢圓的短軸和長軸的比值；zt=γ+ρ(cosθt+iζsinθt)；θt=2π(t-1/2)/NS。當ζ=1時，橢圓即成為圓形；當關注的特征值為實數或離實軸很近時，可以采用較小的ζ，例如ζ=0.1。需要注意的是：式(24)對積分路徑進行了平移和縮放，基于其得到的廣義特征值問題的特征值τj還需通過

λj=γ+ρτj

(25)

變換成原始特征值，同理需要通過

λ′j=ρτ′j

(26)

得到原始特征值的靈敏度值。

2.4 算法的實現步驟

基于2.1節～2.3節的內容，我們可以構建復雜頻變阻抗邊界下的聲學特征頻率靈敏度的計算算法，具體如下：

輸入：NS，K，L，γ，ρ，ζ

輸出：λj，λ′j，ψj，j=1,2,…,m

(1) 初始化隨機矩陣U和V；

(2) 計算矩陣M，M′和S(=0,1,2,…,2K-1)；

(3) 組裝Hankel矩陣H1，H2，H′1，H′2和轉換矩陣S；

(4) 求解H2相對于H1的廣義特征值問題，得到λj，uj和vj(j=1,2,…,m)；

(5) 統計單特征值的數目，記為p，將所有的重特征值及其特征向量置后；

(6) fork=1,2,…,p

(7) 按最大值歸一化法對vk進行正則化處理，并記錄最大值的位置e；

(8) 構造權重向量?k，組裝成系統方程式(14)；

(9) 求解式(14)得到λ′k；

(10) 通過式(10)計算出ψk；

(11) end for

(12) 對每個不同的重特征值，分別求解式(18)和式(10)，重特征值對應的靈敏度及其特征向量。

3 數值算例

本章將使用兩個數值算例來驗證本文所構建的聲學特征頻率靈敏度分析方法的準確性和適用性：第1個算例為具有解析解的二維圓環模型；第2個算例為三維消聲器模型。所有計算程序均采用FORTRAN 90編寫，有限元系數矩陣采用壓縮稀疏行(compressed sparse row,CSR)格式存儲，并通過Intel MKL提供的PARDISO求解器求解有限元稀疏線性方程組以形成矩陣M，M′和S。算例中的聲介質均為空氣，其密度ρm=1.2 kg/m3，聲速c=340 m/s。算例采用在多孔吸聲材料建模中廣泛應用的Delany-Bazley模型獲得等效的聲場阻抗邊界條件，如對于厚度為h且具有剛性背板的多孔吸聲材料，由Delany-Bazley模型得到的表面聲阻抗為

(27)

其中，

(28)

(29)

式中：δ=ρmf/Rf，Rf為多孔材料的流阻率。若h=0.1 m，Rf=10 000 Pa·s/m2，由式(27)得到阻抗Z隨頻率f的變化曲線如圖1所示(為顯示起見，圖中虛部值為原阻抗虛部值的相反數)。從圖1中可以看出，阻抗Z隨頻率f的改變而變化，在0～500 Hz尤為明顯，因此在聲模態分析中考慮頻變阻抗是很有必要的。

圖1 阻抗曲線圖Fig.1 The curve of acoustic impedance

3.1 二維圓環模型

二維圓環模型的示意圖如圖2所示，圖2中的r和R分別為圓環的內徑和外徑。該模型具有解析解，常用作汽車輪胎內聲腔的二維簡化模型。汽車輪胎內聲腔的共振會致使車內噪聲在200～250 Hz出現一個明顯峰值[19-20]，工程中常采用在輪胎壁面鋪設多孔吸聲材料的方法加以解決。圓環模型的內邊界代表輪輞，設置為剛性邊界；外邊界代表輪胎壁面，由于鋪設多孔吸聲材料的緣故，設置為阻抗邊界。該模型的解析特征值為滿足方程

圖2 二維圓環模型示意圖Fig.2 The 2D circular model

J′n(kr)[Yn(kR)-iZ0Y′n(kR)]-Y′n(kr)[Jn(kR)-iZ0J′n(kR)]=0

(30)

的k值(此處的k為波數，需注意以下特征頻率皆以波數形式給出，為了以示區別我們統一稱之為特征值)。式(30)中：Jn和Yn分別為n階第一類和第二類貝塞爾函數；J′n和Y′n分別為Jn和Yn的一階導數。

在數值仿真中，圓環的外徑R=1 m，內徑r=0.5 m；多孔吸聲材料的厚度h=0.1 m，流阻率Rf=10 000 Pa·s/m2。圓環域共離散成1 465個四邊形二次聲學單元。圍道積分路徑選取參數為γ=(3.9,0.5)，ρ=1.8和ζ=0.5的橢圓，梯形積分點數NS=150，參數K和L分別設為7和4。

特征值的計算結果及其實部和虛部的相對誤差，如表1所示。從表1中可以看出，除了第11個特征值為單特征值外，其余皆為重特征值。第1、第3、第5、第7、第11和第12個特征值對應的振型如圖3所示。從表1中還可以看出，本文所構建的聲學非線性特征值分析方法對單特征值和重特征值均可以給出非常準確的計算結果，且所有特征值都具有正虛部，該虛部反映了多孔吸聲材料對聲場的衰減作用，可以用作多孔吸聲材料流阻率等參數的選取判據。特征值靈敏度的計算結果及其實部和虛部的相對誤差，如表2和表3所示。其中：表2是以內徑r作為設計變量；表3是以外徑R作為設計變量。從表2和表3中可以看出，隨著特征值的增大，計算結果實部和虛部的相對誤差逐漸變大，但整體誤差水平都很低，這表明本文所構建的分析方法可以準確地計算出聲學特征頻率的靈敏度值。

表1 數值特征值及其相對誤差Tab.1 Numerical eigenvalues and their relative errors

圖3 圓環形空腔部分振型圖Fig.3 Eigenmodes of the circular model

表2 特征值靈敏度及其相對誤差(內徑r為設計變量)Tab.2 Eigenvalue sensitivities and their relative errors (the design parameter is the inner radius r)

表3 特征值靈敏度及其相對誤差(外徑R為設計變量)Tab.3 Eigenvalue sensitivities and their relative errors (the design parameter is the outer radius R)

3.2 消聲器模型

本算例通過如圖4所示的消聲器模型進一步驗證本文構建的聲學特征頻率靈敏度分析方法在工程問題中的適用性。消聲器軸向的最大尺寸為0.9 m，截面寬和高的最大尺寸分別為0.30 m和0.15 m，沿軸線方向的兩個排氣管的半徑為0.04 m，長度為0.15 m。為了提高降噪性能，在消聲器共振腔沿排氣管軸線方向的周向內壁鋪設有多孔吸聲材料。不同于算例1采用聲腔的形狀參數作為設計變量，本算例將采用多孔吸聲材料的厚度作為設計變量。

圖4 消聲器模型示意圖Fig.4 The muffler model

多孔吸聲材料的厚度h=0.015 m，流阻率Rf=1 424.2 Pa·s/m2，在數值仿真中，同樣采用Delany-Bazley模型將鋪設多孔吸聲材料的壁面等效為聲學阻抗邊界，而其余壁面設置為聲學剛性邊界。經有限元網格劃分后，聲場域共離散成4 575個四面體二次聲學單元。圍道積分路徑選取圓心γ=(9.1,0.2)、半徑ρ=5.0的圓形，梯形積分點數NS=150，參數K和L分別設為5和4。在靈敏度分析中，取多孔材料的厚度h為設計變量。由于該算例沒有解析解，因此采用差分法計算得到靈敏度的參考值，差分步長設置為h/105。特征值及其靈敏度的計算結果，以及靈敏度的差分參考值，如表4所示。部分振型圖如圖5所示。從表4中可以看出，本文構造的聲學特征頻率靈敏度分析方法的計算結果與差分方法得到的參考值非常接近。然而相比差分法，本文方法不存在步長敏感和反復計算等問題，具有更好的準確性和適用性。

表4 消聲器模型的特征值及其靈敏度Tab.4 Eigenvalues and their sensitivities of the muffler model

圖5 消聲器模型的部分振型圖Fig.5 Eigenmodes of the muffler model

4 結 論

針對鋪設多孔吸聲材料聲空間的優化設計，本文基于聲學有限元法，構建了一種聲模態特征頻率靈敏度分析的數值方法，為計算聲學非線性特征值靈敏度提供了一種新思路。該方法首先通過圍道積分法將復雜的非線性特征值靈敏度問題轉化為規模很小的廣義特征值靈敏度問題；然后一方面通過對右特征向量正則化處理實現單特征值靈敏度的計算，提高了方法的計算效率；另一方面采用重構右特征向量空間的方法實現重特征值靈敏度的計算，拓展了方法的應用范圍。數值算例驗證了該方法的求解精度和適用性，以及在工程問題中的應用潛力。另外，該方法還可以很容易地推廣到聲振耦合問題[21-22]的特征頻率靈敏度分析中去。