劉 強(qiáng)，梁夢輝，鄭昌軍，畢傳興

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

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

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

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

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

1 聲學(xué)非線性特征值問題

當(dāng)考慮非黏性理想聲介質(zhì)時，空腔Ω內(nèi)的時諧聲場的頻域控制微分方程為

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

(1)

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

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

(2)

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

(K+ikC-k2M)p=F

(3)

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

(4a)

(4b)

(4c)

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

T(λ)ψ=0

(5)

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

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

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

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

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

為實(shí)現(xiàn)非線性特征值問題的轉(zhuǎn)化求解，我們首先在復(fù)平面內(nèi)選取一個能夠包圍待求特征值區(qū)間的封閉曲線C，然后形成如下一類矩陣

(6)

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

(7)

(8)

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

(H2-λjH1)vj=0

(9)

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

ψj=Svj

(10)

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

(11)

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

2.2 單特征值的靈敏度分析

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

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

(12)

式中，()′=?()/??為對設(shè)計(jì)變量?的導(dǎo)數(shù)。此處的設(shè)計(jì)變量?可以為形狀參數(shù)、材料屬性參數(shù)、拓?fù)鋮?shù)等，本文側(cè)重于形狀參數(shù)的靈敏度分析。

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

?jv′j=0

(13)

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

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

(14)

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

2.3 重特征值的靈敏度分析

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

H2X=H1XΛ

(15)

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

(16)

(17)

WΓ=ΓΛ′

(18)

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

(19)

(20)

由于隨機(jī)矩陣U和V與設(shè)計(jì)變量無關(guān)，因此矩陣M對設(shè)計(jì)變量的導(dǎo)數(shù)為

(21)

式中，Tn(z)=-T(z)-1T′(z)T(z)-1，T′(z)為有限元系數(shù)矩陣T(z)對設(shè)計(jì)變量的導(dǎo)數(shù)。當(dāng)考慮的設(shè)計(jì)變量與頻率無關(guān)時，T′(z)為

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

(22)

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

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

(23)

的圍道積分。當(dāng)關(guān)注的特征值區(qū)間為[λmin,λmax]時，若選取橢圓作為積分路徑，且橢圓的長軸和短軸分別與復(fù)平面的實(shí)軸和虛軸平行，并且采用NS點(diǎn)梯形公式進(jìn)行數(shù)值積分，式(23)可以轉(zhuǎn)換為

(24)

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

λj=γ+ρτj

(25)

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

λ′j=ρτ′j

(26)

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

2.4 算法的實(shí)現(xiàn)步驟

基于2.1節(jié)～2.3節(jié)的內(nèi)容，我們可以構(gòu)建復(fù)雜頻變阻抗邊界下的聲學(xué)特征頻率靈敏度的計(jì)算算法，具體如下：

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

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

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

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

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

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

(5) 統(tǒng)計(jì)單特征值的數(shù)目，記為p，將所有的重特征值及其特征向量置后；

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

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

(8) 構(gòu)造權(quán)重向量?k，組裝成系統(tǒng)方程式(14)；

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

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

(11) end for

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

3 數(shù)值算例

本章將使用兩個數(shù)值算例來驗(yàn)證本文所構(gòu)建的聲學(xué)特征頻率靈敏度分析方法的準(zhǔn)確性和適用性：第1個算例為具有解析解的二維圓環(huán)模型；第2個算例為三維消聲器模型。所有計(jì)算程序均采用FORTRAN 90編寫，有限元系數(shù)矩陣采用壓縮稀疏行(compressed sparse row,CSR)格式存儲，并通過Intel MKL提供的PARDISO求解器求解有限元稀疏線性方程組以形成矩陣M，M′和S。算例中的聲介質(zhì)均為空氣，其密度ρm=1.2 kg/m3，聲速c=340 m/s。算例采用在多孔吸聲材料建模中廣泛應(yīng)用的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所示(為顯示起見，圖中虛部值為原阻抗虛部值的相反數(shù))。從圖1中可以看出，阻抗Z隨頻率f的改變而變化，在0～500 Hz尤為明顯，因此在聲模態(tài)分析中考慮頻變阻抗是很有必要的。

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

3.1 二維圓環(huán)模型

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

圖2 二維圓環(huán)模型示意圖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為波數(shù)，需注意以下特征頻率皆以波數(shù)形式給出，為了以示區(qū)別我們統(tǒng)一稱之為特征值)。式(30)中：Jn和Yn分別為n階第一類和第二類貝塞爾函數(shù)；J′n和Y′n分別為Jn和Yn的一階導(dǎo)數(shù)。

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

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

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

圖3 圓環(huán)形空腔部分振型圖Fig.3 Eigenmodes of the circular model

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

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

3.2 消聲器模型

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

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

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

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

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

4 結(jié) 論

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