方春慧，劉鳳景，周 予，方 智

（1.煙臺汽車工程職業(yè)學(xué)院，山東 煙臺 265500；2.盛寶油壓機械（煙臺）有限公司，山東 煙臺 265500；3.華中科技大學(xué) 船舶與海洋工程學(xué)院，武漢 430074）

求解消聲器橫向模態(tài)常用的方法有解析方法和數(shù)值方法。解析方法[1–2]求解簡單、快速、精確，但是僅局限于有解析方程的橫截面，比如圓形截面，矩形截面以及橢圓形截面。對于任意形狀的橫截面橫向模態(tài)的求解，只能使用數(shù)值方法[3–4]，現(xiàn)在較為成熟的常用的數(shù)值方法為二維有限元方法（FEM）和邊界元方法(BEM)。該類方法均需要依賴于網(wǎng)格進(jìn)行求解，是基于網(wǎng)格的求解方法。無網(wǎng)格方法(MFM)[5]是一種相對較新的數(shù)值方法，求解時節(jié)點是相對于網(wǎng)格自由存在的，甚至不需要網(wǎng)格，本文所用到的徑向基函數(shù)點插值配點方法是一種真正的無網(wǎng)格方法，場節(jié)點自由任意分布在問題域和邊界上。無網(wǎng)格方法最先在力學(xué)領(lǐng)域應(yīng)用較多，近十幾年開始使用在聲學(xué)領(lǐng)域。1998年，Bouillard應(yīng)用基于移動最小二乘的無網(wǎng)格伽遼金法計算了二維內(nèi)部聲學(xué)問題，數(shù)值算例比較結(jié)果表明該方法比有限元法的計算精度高，計算速度快[6]。2002年，Chen等應(yīng)用邊界配點的無網(wǎng)格方法計算了二維聲腔的聲學(xué)特征值問題，并指出該方法較有限元法簡單，容易實施[7]。2010年，大連理工大學(xué)的張宏偉等應(yīng)用基于點插值的配點型無網(wǎng)格法求解了Helmholtz方程，通過數(shù)值算例比較證明了該方法具有較高的精度和良好的收斂性[8]。2011年，湖南大學(xué)的姚凌云等將分區(qū)光滑徑向點插值無網(wǎng)格法應(yīng)用于二維聲學(xué)分析中，對管道和二維轎車聲學(xué)問題的數(shù)值分析結(jié)果表明，該方法比有限元法具有更高的精度[9]。2012年，海軍工程大學(xué)的胡海等將邊界無網(wǎng)格法應(yīng)用到結(jié)構(gòu)聲輻射計算中，通過與解析結(jié)果對比表明該方法具有更高的插值和計算精度[10]。華中科技大學(xué)的黃其柏課題組將徑向配點無網(wǎng)格法應(yīng)用到氣動聲學(xué)和腔體結(jié)構(gòu)聲學(xué)的研究中，取得很好的成果[11]。Christina研究了使用徑向點插值方法（RPIM）求解二維Helmholtz方程時的色散效應(yīng)[12]。但是當(dāng)導(dǎo)數(shù)邊界條件存在時，使用RPIM求解會出現(xiàn)不穩(wěn)定。鑒于此，一種新型的處理方法：Hermite RPIM方法被提出處理導(dǎo)數(shù)邊界條件。

本文使用多項式Hermite徑向基函數(shù)點插值法求解問題域內(nèi)計算點的形函數(shù)，使用配點方法離散二維橫向模態(tài)方程，進(jìn)而求解消聲器橫向本征波數(shù)以及模態(tài)形狀，研究形狀參數(shù)對計算精度和速度的影響。

1 Hermite型RPIM求解形函數(shù)

對于膨脹腔消聲器，二維聲壓控制方程可以表示為

其中pxy為橫向聲壓，為二維笛卡爾坐標(biāo)系的拉普拉斯算子，kxy為橫向方向的本征波數(shù)。

對于剛性壁，橫向邊界條件可以表示為

使用Hermite徑向基函數(shù)點插值法構(gòu)造形函數(shù)。場節(jié)點隨機分布在問題域Ω和邊界Γ上。對于任意一個計算點，可以形成一個影響域。影響域可以是任意形狀，本文選取圓形。問題域和影響域的示意圖如圖1所示。假設(shè)在影響域內(nèi)有n個場節(jié)點（包括導(dǎo)數(shù)邊界上的節(jié)點），導(dǎo)數(shù)邊界上的節(jié)點數(shù)目為nDB。計算點的聲壓ph可以寫成如下形式

圖1 問題域和影響域

其中Ri(x)是徑向基函數(shù)，dc為位于計算點附近的節(jié)點間距。αc和q是無量綱的形狀參數(shù)。xi是影響域內(nèi)每個節(jié)點的坐標(biāo)，ai是相應(yīng)的系數(shù)。是位于導(dǎo)數(shù)邊界上的點的坐標(biāo)，bj是相應(yīng)的系數(shù)。qk是多項式，m為多項式的項數(shù)，ck是多項式的系數(shù)。n為導(dǎo)數(shù)邊界上第j個計算點的單位外法線向量，lxj和lyj為方向余弦。

將方程(3)寫成如下矩陣形式

方程式(3)中的系數(shù)可采用通過影響域中的所有n個節(jié)點的函數(shù)值的插值以及等于導(dǎo)數(shù)邊界節(jié)點的函數(shù)導(dǎo)數(shù)值而確定。

對局部支持域中的所有n個節(jié)點（包括導(dǎo)數(shù)邊界上的節(jié)點）可得到

式中l(wèi)=1,2…,n。

對所有導(dǎo)數(shù)邊界上的節(jié)點可得到

式中l(wèi)=1,2…,nDB。

為獲得唯一解，可以施加以下約束條件

聯(lián)立方程式(10)至式(12)可得到如下矩陣方程

其中，廣義力矩矩陣G為對稱矩陣，其中每個矩陣可以表示為

求解方程式(13)可以得到聲壓表達(dá)式為

其中，形函數(shù)向量φ可以寫成

進(jìn)而，可以得到聲壓近似函數(shù)及其導(dǎo)數(shù)表達(dá)式

2 配點法求解橫向模態(tài)

一般問題域內(nèi)場節(jié)點是任意分布的，為了方便起見，配點取與場節(jié)點相同的點。假設(shè)問題域內(nèi)部節(jié)點數(shù)為Nd，邊界節(jié)點數(shù)為Nb，總的節(jié)點數(shù)為N=Nb+Nd。若所研究的影響域內(nèi)包含導(dǎo)數(shù)邊界上的節(jié)點，對于影響域內(nèi)部節(jié)點xI，將式(20)代入到式(1)中可以得到以下矩陣方程式

KI和MI分別是影響域內(nèi)部節(jié)點xI的聲剛度和聲質(zhì)量矩陣。值得注意的是，在所有內(nèi)部節(jié)點都需建立方程式(25)。

導(dǎo)數(shù)邊界上的節(jié)點需要滿足兩個方程，一個為方程式(25)，另一個為如下形式

值得注意的是，在所有導(dǎo)數(shù)邊界節(jié)點上均需要建立方程式(28)。

若所研究的影響域內(nèi)不包含邊界節(jié)點，則方程式(25)和式(28)均可簡化。

組裝所有節(jié)點聲剛度和聲質(zhì)量矩陣，可得到如下系統(tǒng)本征方程。

與傳統(tǒng)有限元方法不同的是，在無網(wǎng)格配點方法中，矩陣組裝的具體方法是將其節(jié)點矩陣按行擺放在一起組成總體矩陣。

求解方程式(31)可以得到相應(yīng)的本征波數(shù)和本征向量矩陣。

3 結(jié)果與討論

由理論分析可知，本征波數(shù)計算誤差主要來自三部分：徑向基函數(shù)形參的選擇，導(dǎo)數(shù)邊界條件的處理以及配點布置。鑒于圓形截面的本征波數(shù)有解析解，本文采用Hermite型RPIM配點法求解圖2所示的圓形截面的本征波數(shù)，并進(jìn)行其誤差分析。圓形橫截面的半徑取r=0.024 3 m。本征波數(shù)數(shù)值結(jié)果與解析結(jié)果的相對誤差定義為

如圖2所示，圓形截面內(nèi)的節(jié)點分布是任意的，由MATLAB編寫簡單程序得到。

圖2 圓形橫截面示意圖

3.1 場節(jié)點數(shù)目的影響

場節(jié)點數(shù)目對前6階橫向本征波數(shù)的計算誤差的影響如圖3和圖4所示。

圖3 問題域內(nèi)部節(jié)點數(shù)目的影響

圖4 問題域邊界節(jié)點數(shù)目的影響

從圖中可以看出，本征波數(shù)的計算精度并沒有隨著節(jié)點數(shù)目的增加而增大。對于本文研究的圓形橫截面，當(dāng)內(nèi)部節(jié)點數(shù)目超過100，邊界節(jié)點數(shù)目超過60，本征波數(shù)的計算誤差反而會增大，這是因為隨著節(jié)點增多，廣義力矩陣G的條件數(shù)也會隨之增大，會造成計算結(jié)果的誤差增大。所以在計橫截面本征波數(shù)時，需要選擇合適的節(jié)點數(shù)目。對于圖2所示的圓形橫截面，內(nèi)部節(jié)點和邊界節(jié)點分別取100和40最為合適。

3.2 影響域尺寸的影響

圖5給出了影響域內(nèi)包含不同個數(shù)節(jié)點時的本征波數(shù)的計算誤差值。與問題域內(nèi)場節(jié)點數(shù)目的影響趨勢大致相同，本征波數(shù)的計算精度并不是隨著影響域增大而增大。

圖5 影響域尺寸的影響

由圖5可以看出，當(dāng)影響域內(nèi)節(jié)點數(shù)目超過35，計算誤差反而增大，這也是由于隨著影響域尺寸的增大，廣義力矩陣條件數(shù)增大的原因。所以在選擇計算點的影響域尺寸時，需要平衡計算精度和計算速度的關(guān)系。

3.3 徑向基函數(shù)形參的影響

圖6和圖7給出了徑向基函數(shù)形參對本征波數(shù)計算誤差的影響。形參的選擇對計算誤差的影響比較大，本文通過多次多參數(shù)嘗試，選擇了幾個常用的有相互關(guān)聯(lián)的形狀參數(shù)q和ac，并且研究了它們的影響。

由圖6可以看出，當(dāng)q=1時，相對誤差基本上為100%，數(shù)值計算結(jié)果不正確。當(dāng)q=±0.5時，復(fù)合二次徑向基函數(shù)為標(biāo)準(zhǔn)徑向基函數(shù)，但是數(shù)值計算結(jié)果并沒達(dá)到最大的精度，反而當(dāng)q=0.98和q=1.03時本征函數(shù)的計算誤差最小。

由圖7可以看出，計算誤差隨著ac的增大而減小，但是存在一個計算誤差最小的形參值ac=3.5?？梢缘贸鼋Y(jié)論，對于本文研究的二維聲學(xué)問題，最優(yōu)的徑向基函數(shù)形參可以取為q=0.8或者q=1.03和ac=0.35。但是需要注意的是，對于不同的問題，甚至是相同的問題不同的橫截面形狀和尺寸，最優(yōu)的形狀參數(shù)的取值會發(fā)生變化。

圖6 形參q的影響

圖7 形參ac的影響

3.4 計算速度的比較

為了驗證本文方法和編寫的計算程序的正確性和適用性，本文分別使用二維有限元方法和本文的無網(wǎng)格方法計算了如圖8所示的不規(guī)則結(jié)構(gòu)跑道圓橫截面的本征頻率和模態(tài)形狀。跑道圓的尺寸為長短軸分別為a=0.25和b=0.15 m。橫截面的本征波數(shù)和模態(tài)形狀分別如圖9和圖10所示。

圖8 跑道圓橫截面示意圖

由圖可以看出，有限元方法和無網(wǎng)格方法計算的本征頻率的相對誤差在5%以內(nèi)，再次驗證了本文方法的正確性。

圖9 跑道圓截面本征頻率

圖10 跑道圓截面模態(tài)形狀

為了比較本文方法的計算速度，表1給出了分別使用有限元方法和無網(wǎng)格方法計算圖2和圖8所示橫截面的時間比較。為了得到相同的精度，有限元方法和無網(wǎng)格方法計算本征模態(tài)時選取相同的節(jié)點數(shù)，本文中圓形截面節(jié)點數(shù)設(shè)置為140，跑道圓截面節(jié)點數(shù)為1 140。表1中所示的計算時間為求解本征波數(shù)的時間，不包括前處理時間，比如有限元方法中生成網(wǎng)格的時間，無網(wǎng)格方法中前處理時間即節(jié)點生成的時間可以忽略不計，因為節(jié)點可以通過簡單的MATLAB語句生成。由表1可以知道，無網(wǎng)格方法比有限元方法節(jié)省時間，尤其是計算較大復(fù)雜的橫截面時無網(wǎng)格方法的快速性更為突出。

但是，由前面的分析可知，影響無網(wǎng)格方法計算精度的因素很多，尤其是徑向基函數(shù)形狀參數(shù)的選擇，不同的計算結(jié)構(gòu)最優(yōu)的形參取值不同，域內(nèi)節(jié)點數(shù)目以及每個計算點影響域的尺寸的選取也將不同，這也是與有限元方法相比無網(wǎng)格方法的不足之處。

表1 本征波數(shù)求解速度比較/min

3.5 跑道圓形雙出口膨脹腔消聲器橫向模態(tài)計算

為了體現(xiàn)本文方法的有效性，使用該方法計算了如圖11所示的跑道圓形雙出口插管膨脹腔消聲器每個特征橫截面的本征模態(tài)。膨脹腔跑道圓橫截面的長短軸長度分別為a=0.218 m和b=0.109 m，進(jìn)出口管直徑為d1=d2=d3=0.048 6 m，兩出口偏移量均為δ1=δ2=0.054 5 m 。該消聲器共有SA、SB、SC、SD、SE、SF、SG7個特征橫截面，由于對稱性以及進(jìn)出口管尺寸相同，所以特征橫截面SA、SE、SG具有相同的本征模態(tài)。表2列出了本文無網(wǎng)格方法計算的本征頻率和模態(tài)形狀。

圖11 跑道圓形雙出口插管膨脹腔消聲器

表2 消聲器特征截面橫向模態(tài)

4 結(jié)語

本文提出Hermite型徑向基函數(shù)點插值配點方法求解消聲器的橫向模態(tài)。通過比較圓形和跑道圓形橫向本征波數(shù)，驗證了本文提出的方法的計算精度和速度。進(jìn)而分析了問題域內(nèi)場節(jié)點數(shù)目，影響域的尺寸以及徑向基函數(shù)的形狀參數(shù)對本征波數(shù)計算誤差的影響。對于每一個影響因素的選取，都存在一個最優(yōu)值。雖然相比有限元方法，本文提出的徑向基函數(shù)插值配點法計算精度和計算速度均具有優(yōu)勢，但是對于每一個影響因素最優(yōu)取值的選取，則需要重復(fù)性地嘗試，這也是該方法不足之處。

[1]SELAMET A,XU M B,LEE I-J.Analytical approach for sound attenuation in perforated dissipative silencers with inlet/outlet extension[J].J.Acoust.Soc.Am.,2005,117(4):2078-2089.

[2]KIRBY R,DENIA F D.Analytic mode matching for a circular dissipative silencer containing mean flow and a perforated pipe[J].J.Acoust.Soc.Am.,2007,122(6),3471-3482.

[3]ZHAO S.On the spurious solutions in the high-order finite element methods for eigenvalue problems[J].Computer Methods in Applied Mechanics and Engineering,2007,196(49-52):5031-5036.

[4]PETYT M,KOOPMAN G H,PINNINGTON R J.The acoustic modes of a rectangular cavity with a rigid incomplete partition[J].J.Sound.Vib,1997,53(1):71-82.

[5]LIU G R,GU Y T.An introduction to meshfree methods and their programming[M].New York:Springer,2005.

[6]BOUILLARD P,SULEAUB S.Element-free galerkin solutions for helmholtz problems: formulation and numerical assessment of the pollution effect[J].Computer Methods in Applied Mechanics and Engineering,1998,162,317-335，

[7]CHEN J T,CHANG M H,CHEN K H,et al.Boundary collocation method for acoustic eigen analysis of threedimensionalcavities using radialbasis function[J].Computational Mechanics,2002,29(4):392-408.

[8]李美香，張宏偉，李衛(wèi)國，基于點插值的配點型無網(wǎng)格方法解 Helmholtz問題[J].計算力學(xué)學(xué)報，2010，27（3）：533-536.

[9]姚凌云，于德介，臧獻(xiàn)國.聲學(xué)數(shù)值計算的分區(qū)光滑徑向點插值無網(wǎng)格法[J].振動與沖擊，2011，30(10)：188-192.

[10]胡海，郭文勇，馬龍.結(jié)構(gòu)聲輻射計算的邊界無網(wǎng)格法[J].噪聲與振動控制，2012，32(5)：37-41.

[11]LI K,HUANG Q B,WANG J L,et al.An improved localized radial basis function meshless method for computational aeroacoustics[J].Engineering Analysis with Boundary Elements,2011,35(1):47-55.

[12]CHRISTINA W,OTTO V E.Dispersion analysis of the meshfree radialpointinterpolation method forthe Helmholtz equation[J].Int J Numer Meth Eng,2009,77(12):1670-89.