廣義有限差分法在含阻抗邊界空腔聲學分析中的應用1)

陳增濤 王發杰,?,2) 王 超

?(青島大學動力集成及儲能系統工程技術中心,青島大學機電工程學院,山東青島266071)

?(青島大學多功能材料與結構力學研究院,山東青島266071)

引言

目前,有限元、有限差分和邊界元法作為空腔聲場模擬計算中的主要方法,在汽車發動機傳熱或車內噪聲控制、封閉管道噪聲聲場等方面均有廣泛應用.然而對于一些多維復雜幾何域問題,合理有效的模擬計算很難準確實現.不同于以上依賴網格的數值計算方法,無網格法通過節點信息建立插值基函數,插值基函數不依賴于節點之間的有序拓撲鏈接,因此相對于傳統方法不會受到網格劃分的約束.近年來無網格方法的研究受到了高度重視,成為目前科學與工程計算中的研究熱點之一[1-6].

根據插值基函數的類型不同,無網格法可分為邊界型[7-8]和區域型[9-10]兩大類.邊界型無網格法保留了降低所求問題維度的優點,并且避免了復雜的奇異和近奇異積分的計算.近年來,邊界型無網格方法引起了廣泛的關注,主導方法有邊界節點法[11-12]、基本解法[13-14]、邊界無單元法[15]、奇異邊界法[16-18]、正則化無網格方法[19]等.以上方法均在多種領域成功運用并取得了不同程度的成功,但是仍存在一些決定性的問題.例如,該類方法在實際工程中的應用范圍[20]很大程度上會受到插值基函數的影響.

相對于邊界型無網格法,區域型無網格方法擺脫了插值基函數的約束,所以對各種偏微分方程初邊值問題具有普遍適用性,應用領域更為廣泛.現如今,比較常用的區域型無網格法包括光滑粒子水動力學法[21]、擴散單元法[22]、無單元伽遼金法[23-24]和廣義有限差分法(Generalized Finite Difference Method,GFDM)[25]等.其中,廣義有限差分法是一種最近發展起來的新型無網格方法,已被成功應用于生物傳熱分析[26]、薄板彎曲問題[27]、分數階反常擴散[28]、熱彈性[29]和反問題[30]等領域.

廣義有限差分法不僅保留了有限差分法把微分方程轉變成差分方程的優點,而且引入了點簇的概念,在模擬復雜的幾何形狀問題時十分具有優勢.除此之外,廣義有限差分法只需要在物理域內隨機布置一組節點,通過相鄰節點之間的關聯得到中心節點的函數值,節省了有限元和邊界元法中費時費力的網格劃分及數值積分,大大節省了數據準備的時間,在人工越來越昂貴的情況下,這也是非常重要的優勢.

本文首次將廣義有限差分法應用于含阻抗邊界空腔聲學問題的分析和模擬,建立了空腔聲場問題的廣義有限差分法數值離散格式.數值試驗測試了經典的二維車腔噪聲模型和三維空腔聲學模型,并與有限元法的數值結果對比來驗證該方法的有效性、精確度和穩定性.本工作為含阻抗邊界空腔聲學問題的仿真提供了新的、簡單高效的計算工具,也拓展了廣義有限差分法的應用領域.

1 空腔聲場問題的基本理論

聲波在傳播過程中會不斷被拉伸和壓縮,動作非常快,質點間來不及進行熱交換.因此聲波的傳播過程可被看成是絕熱過程[31],其熱力學方程為

其中,p0是空氣靜態壓強,ρ0為微元體初始密度,γ 是空氣的比熱比,將式(1)對時間求導可得

推導可以得到在x方向上的一維聲學方程

結合x、y和z三個方向的一維聲學方程,可以得到三維系統的聲學運動方程

在上式中引入拉普拉斯算子?2

則式(4)可以寫成

對簡諧激勵穩態聲場,利用變換P(x,t)=p(x)eiωt,可將式(3)中的時域聲場方程變換為頻域聲場方程,即Helmholtz 波動方程

其中,K=ω/c稱為波數,ω 為圓頻率,c為聲速.上述方程描述了頻域內穩態聲場在空間的分布規律.對于二維和三維聲場分布,同理可得到與(7) 式相同形式的控制方程.聲壓與簡諧聲波的振動速度關系式為

在已知邊界條件下,可定量對其進行求解.聲場有如下幾種邊界條件:

(1) 狄利克雷(Dirichlet) 邊界條件,即聲壓邊界條件

(2)黎曼(Neumann)邊界條件,即速度邊界條件

(3)洛平(Robin)邊界條件,即聲阻抗邊界條件

在式(9)～式(11)中,pD是邊界處聲壓值,An是聲導納系數.

2 廣義有限差分法的基本理論

根據廣義有限差分法(GFDM)計算空腔聲場問題的基本理論,對求解區域布置節點進行離散,其方法是基于移動最小二乘法與泰勒展開.首先在整個計算區域內布置N=ni+nd+nn+nr(ni表示內部節點數目,nd,nn,nr分別表示Dirichlet,Neumann 和Robin 邊界所對應的節點數目) 個節點,再將每個節點上的偏微分項轉換成由子區域內各節點物理量與權重系數乘積的線性累加.以二維問題為例,對區域內的第i個節點而言,選擇m個支撐點(鄰近點),形成一個子區域,如圖1 所示.

圖1 廣義有限差分法配點布置示意圖Fig.1 Schematic diagram of GFDM

通過觀察上述各式可以看出中心節點處的微分項可以變換為由子區域中各點上的物理量與權重系數乘積的線性累加.由此可知,廣義有限差分法的核心思想是將第i點上的微分量用對應子區域內支撐點上的物理量線性權重累加表示,使所有內部點滿足控制方程式,所有邊界點滿足邊界條件,形成一個大型的稀疏線性方程組,對該方程組求解即可得到計算區域內全部點上的物理量.

(1)對于計算區域內部節點(xi,yi),i=1,...,ni,其上的物理量必定滿足聲學控制方程.由式(20)可知,只需提取相應的微分項,然后將線性組合的方程式代入控制方程,則

從而得到最終的計算結果.

從上述數值實施過程中可以看出,GFDM 與傳統的有限差分法十分類似,方程(32) 中的系數矩陣C為稀疏矩陣.另一方面,從形成最終線性方程組(32)的過程可以觀察到,系數矩陣C是一個N行N列的方陣,而每一行最多僅有m個非零元素,這里N為計算區域內全部節點的個數,m為局部支撐節點數目.顯然,m遠小于N,因此方程組(32)是一個稀疏線性方程組.由此可見,與其他具有稠密矩陣方程的數值方法相比,GFDM 能夠非常有效地分析高維和復雜幾何形狀的問題.

3 數值算例

為了驗證廣義有限差分法求解含阻抗邊界空腔聲學問題的有效性和可行性,現針對二維簡單聲腔結構和二維汽車駕駛艙以及三維聲腔模型進行研究,重點探討無網格法和有限元方法求解不同頻率下聲腔的聲場分布情況.需要指出的是,在下面所有算例中,腔體內的介質均為空氣,其密度ρ 為1.225 kg/m3,聲波在空氣中傳播的速度c為343m/s,為節省篇幅,這些參數在下述算例中將不再具體定義.除特殊說明外,數值結果均在i5-5200 CPU@2.20GHz,內存4GB的計算機上運算所得.文中有限元結果均來自于有限元仿真軟件COMSOL Multiphysics 5.4.

為了評估數值誤差,采用最大相對誤差(MRE)和均方根誤差(RMS E)

其中NT是測試點的數量,pexa和pnum分別表示第i個測試點的解析解和數值解.這里,除特殊說明外測試點均為區域內部節點.

3.1 混合邊界管道聲場分析

現求解如圖2 所示的二維矩形聲腔結構,其長度L=1.0 m,寬度H=0.2 m,在A 端口邊界給定Neumann 速度邊界,法向速度vn=0.01 m/s,管道B端口是剛性壁,其他邊界為聲學硬邊界條件.

圖2 管道聲學模型Fig.2 Pipe acoustic model

管道聲壓的精確解為

按照GFDM 求解此類問題的前處理過程,首先將該模型離散成17 996 個點,即總節點數N=17 996,計算時選用的內部支撐點m=12.圖3 給出了當頻率分別為400 Hz 和600 Hz 時,計算得到管道底部軸線上的聲壓分布.分析圖3 聲壓曲線圖可知,不同頻率下GFDM 所得數值解與解析解均非常吻合.

圖3 不同頻率下管道底部軸線聲壓分布Fig.3 Distributions of sound pressure at the lower boundary under different frequencies

為了檢驗GFDM 中總節點數目和內部支撐點數目對計算結果的影響,研究了頻率為400 Hz 時,不同內部支撐點數目和總節點數目對計算方法的影響.分析表1 計算所得的數據可知,總節點數N=17 996時,隨著內部支撐點數目m的增大,GFDM 的計算最大誤差和均方根誤差隨之變小,也就是此方法的計算誤差隨之收斂.觀察表2 數據可知,在內部支撐點數m=12 時,隨著模型的總節點數目的增加,計算結果的最大誤差和均方根誤差都會隨之減小.可見,模型離散的點數越多,GFDM 的計算結果越精確,表明算法具有較好的計算精度和收斂性.

表1 總節點數N=17 996 時,不同支撐點數下的最大相對誤差、均方根誤差和計算時間(CT)Table 1 When N=17 996,the MREs,RMS Es and calculation time(CT)under different number of supporting nodes

表2 局部支撐點數為m=12 時,不同總節點數下的最大相對誤差、均方根誤差和計算時間(CT)Table 2 When m=12,the MREs,RMS Es and calculation time(CT)under different number of total nodes

為從數值上分析最大計算頻率與節點間距(對應于離散點數)之間的關系,我們選取不同的頻率和節點間距?h=,將GFDM 計算所得的RMS E誤差曲面繪制于圖4(a)中.觀察圖中誤差分布情況,我們可以清楚地看到計算誤差隨著頻率和節點間距的同時增大而增大.為更好地分析最大計算頻率與節點間距之間的規律,圖4(b) 給出了對應的二維誤差分布情況,圖中紅色圓點表示誤差大小為RMS E=1.0×10?2的位置.通過曲線擬合的思想,得到了誤差為RMS E=1.0×10?2所在點構成的一條擬合曲線(圖中綠色曲線):f=35·?h(?0.65).至此,我們便得到了一個估計最大計算頻率的經驗公式

需要指出的是,根據經驗公式(38),在保證計算精度為RMS E1.0×10?2的情況下,我們可以輕松地估計出能夠計算的最大頻率.這里僅從數值的角度獲得了一個計算頻率與節點間距關系的經驗公式,理論上獲得準確的估計公式目前仍面臨巨大挑戰,還需進一步研究和探索.

圖4 不同頻率和節點間距下GFDM 所得計算誤差RMS E 的分布Fig.4 Distribution of RMS Es under different frequencies and node intervals

為了考察GFDM 對高頻聲學問題的計算精度和計算穩定性,我們選取計算頻率為20 000 Hz,分析其計算結果.通過經驗公式(38),我們能夠估計出計算頻率為20 000 Hz 時取得良好精度(RMS E<1.0×10?2)所需的節點間距為?h=0.000 055 m,此時,該模型大概需要66 000 000 個離散節點.受工作條件和計算機內存的限制,想要模擬如此大規模的聲學問題,仍然面臨巨大挑戰,即便使用另外一臺較大內存的計算機(i7-6700 CPU@3.40GHz,內存16 GB)進行模擬,4 999 996 個節點(節點間距?h=0.000 2 m)就已經達到了計算機的內存極限.在該離散節點數目下的計算結果如圖5 所示,可以看出GFDM 所得結果與精確解(exact)十分吻合.計算過程花費了3990.333 s的計算時間,計算誤差為RMS E=0.087 362.另外,從圖5 中可以看出,高頻問題對應的聲壓波動非常劇烈,這類問題的精確模擬對大多數數值算法來說都是比較棘手的,上述數值結果表明GFDM 能夠有效模擬高頻問題,對該類問題具有較高的計算精度和計算穩定性.

圖5 頻率為20 000 Hz 時管道底部軸線聲壓分布Fig.5 Distributions of sound pressure at the lower boundary under f=20 000 Hz

3.2 二維車腔聲場分析

在整車的概念設計階段,車內聲腔可近似地被看作為棱柱體,因此可將其簡化為二維模型進行噪聲分析.圖6 所示為某車簡化之后的車內聲腔模型,其水平方向的長度L=3.41 m,豎直方向的高度H=1.21 m,相關條件已在圖中標示[35].腔體內部介質為空氣,由于發動機振動激勵引起的噪聲占整車噪聲的絕大部分,因而在汽車的前圍板上施加一Neumann 邊界條件,其振動速度為0.01 m/s.此外,我們也在乘員艙的頂部cd 處施加一Robin 邊界條件以模擬內飾材料的阻尼,這里聲阻抗值為413 Pa·s/m3.在算例分析中,選取底部軸線ab 上的聲壓數據進行研究.圖7 給出了模型的GFDM 節點分布示意圖.

圖6 具有阻抗邊界條件的二維車腔模型Fig.6 A two-dimensional car cavity with impedance boundary

圖7 車腔的GFDM 節點分布Fig.7 Distribution of nodes on the car cavity

車內噪聲按其頻率高低一般可分為低頻噪聲(1 ～200 Hz)、中頻噪聲(200 ～500 Hz) 和高頻噪聲(>500 Hz).在實際駕乘過程中,接觸最多的是中低頻噪聲.因此,我們首先分析了低頻(f=200 Hz)時的聲壓分布,為驗證算法的有效性和準確度,我們把細密網格下有限元(FEM)的數值解作為參考解.采用有限元法進行聲學數值計算時,為了保證其結果的準確性,需滿足聲學“經驗法則”,即最小波長內至少有6 個網格單元.在接下來的計算分析中,首先將汽車腔體模型在有限元軟件中劃分為43 480 個三角形單元,網格尺寸為0.01 m,即每波長內包含171個網格單元.使用GFDM 計算時,將算例模型離散成11 558 個離散節點,即N=11 558,節點間距為0.02 m.圖8 給出了頻率較低f=200 Hz 時,內部支撐點數m=12,總節點數N=11 558 時聲腔模型內聲壓分布云圖.從圖8 聲壓分布云圖可以看出,即使在較低頻率下,FEM 和GFDM 計算的結果整體上吻合良好,表明了GFDM 使用低于有限元求解經驗法則的節點數就可以獲得比較理想的求解精度,也體現了該方法在降低計算內存消耗和提高數值計算效率方面的潛能.

圖8 有限元和廣義有限差分法計算所得聲壓分布云圖Fig.8 Contours of sound pressure obtained by the FEM and GFDM

為驗證GFDM 在實際問題中的計算收斂性和穩定性,在GFDM 計算所用點數遠小于FEM 離散節點的情況下,將三組不同的總節點數N計算所得數值結果與FEM 參考解對比分析,圖9 和圖10 分別給出了400 Hz 和600 Hz 頻率下底部邊界上的GFDM 和FEM 的對比情況.圖9 展示了不同頻率下的聲壓曲線,圖10 是不同頻率下的聲壓級曲線.

從圖9 和圖10 可以看出,GFDM 在不同節點數下計算所得聲壓和聲壓級均與有限元結果很大程度上吻合,即使采用較少的節點數目(N=11 558),GFDM 的計算結果也與有限元結果非常接近,表明算法是準確有效的.另外,從圖9 和圖10 中可以看出,GFDM 的計算精度隨著計算節點數目的增加而提高,說明具有較好的收斂性.綜上,GFDM 在求解含阻抗邊界車腔聲學問題時具有較好的計算穩定性和收斂性.相比于FEM,GFDM 無需網格劃分等繁瑣的前處理過程,僅需離散節點,所以可作為聲腔分析的一種有效工具.

圖9 頻率為400 Hz 和600 Hz 時有限元法與廣義有限差分法所得聲壓曲線的對比Fig.9 Comparisons of sound pressure obtained by the FEM and GFDM,under f=400 Hz and 600 Hz

圖10 頻率為400 Hz 和600 Hz 時有限元法與廣義有限差分法所得聲壓級曲線的對比Fig.10 Comparisons of sound pressure level obtained by the FEM and GFDM,under f=400 Hz and 600 Hz

3.3 三維空腔聲場分析

為了檢驗GFDM 計算三維復雜聲場問題的能力,下面對三維空腔模型的聲學特性仿真計算[36],該聲腔的介質為空氣,模型的A 端面處帶有初始速度條件,B 端面處是吸聲材料,數值分析中所施加的Neumann 及Robin 邊界條件也同上一算例一致.在有限元分析軟件中,將聲腔模型較細化劃分為51 469 個單元,單元平均尺寸為0.011 m,將頻率為400 Hz 下(根據聲學“經驗法則”,每波長內分布78 個網格單元)的有限元結果作為參考值.三維空腔結構的尺寸和布點方式如圖11 所示.

圖11 三維空腔及其GFDM 節點分布Fig.11 Three-dimensional cavity and the distribution of nodes

使用FEM 和GFDM 分別計算模型邊界點上的聲壓,圖12 和圖13 分別給出了兩種方法計算400 Hz時的聲壓和聲壓級分布云圖.從圖中發現,GFDM 計算的數值結果和FEM 參考值十分吻合,可見GFDM在三維聲腔模型中也有較好的計算精度.

圖12 400 Hz 下有限元和廣義有限差分法所得聲壓分布云圖Fig.12 Contours of sound pressure obtained by the FEM and GFDM under f=400 Hz

圖13 400 Hz 下有限元和廣義有限差分法所得聲壓級分布云圖Fig.13 Contours of sound pressure level obtained by the FEM and GFDM under f=400 Hz

選取圖11 中的F(0.05,0.05,0.12)點作為監測點來模擬人的耳朵所在位置,分析該處聲壓級隨頻率增加的變化情況.圖14 給出了GFDM 和FEM 所得監測點F 處聲壓級頻率響應曲線的對比.盡管在中頻(700 Hz)附近和較高頻率(>1000 Hz)附近出現了微小的偏差,但GFDM 的計算結果和FEM 參考解仍然比較吻合,充分表明了算法的有效性和準確性.需要指出的是,增加GFDM 的離散節點數目,可以使得兩種算法的數值結果更加接近.

圖14 監測點F 處的頻率響應曲線Fig.14 Frequency response curves monitoring point F

3.4 三維S 型圓柱管道聲場分析

圓柱管道是生活中比較常見的空腔模型,為了進一步驗證本文方法的計算準確性,下面對S 型管道分別采用FEM 和GFDM 計算并對其聲壓分析.三維管道模型直徑為D=0.1 m,長度L=1.0 m,a端面處施加Neumann 速度邊界vn=0.01 m/s,b端面處有吸聲材料,阻抗值Z=413 Pa·s/m3.其余邊界條件為剛性壁,=0.模型尺寸和布點格式如圖15 所示.

圖15 S 型聲腔及其節點分布Fig.15 S-shaped acoustic cavity and nodal distribution

為了保證GFDM 對高頻問題的有效性,我們將模型離散成 126 905 個節點,此時的節點間距0.003 5 m.同時,將具有729 224 個網格單元(包含86 624 個邊界單元和642 600 個域單元)和網格平均大小為0.003 m 的FEM 結果作為參考解.圖16 和圖17 分別給出了低頻400 Hz 和高頻1600 Hz 頻率下管道表面聲壓分布云圖以及絕對偏差云圖.分析聲壓云圖中可以發現,GFDM 計算的結果和FEM 參考解幾乎完全吻合,并且從絕對偏差云圖可以看出,400 Hz 下的偏差均小于0.06,1600 Hz 下的偏差均小于0.45,這充分說明了方法的可行性和有效性.因此,GFDM 在模擬含阻抗邊界空腔問題時的計算精度和數值穩定性是比較理想的.

圖16 400 Hz 下有限元和廣義有限差分法所得聲壓和絕對偏差的分布云圖Fig.16 Contours of sound pressure and absolute deviation obtained by the FEM and GFDM under f=400 Hz

圖17 1600 Hz 下有限元和廣義有限差分法所得聲壓和絕對偏差分布云圖Fig.17 Contours of sound pressure and absolute deviation obtained by the FEM and GFDM under f=1600 Hz

4 結論

廣義有限差分法(GFDM)相比于傳統網格類數值方法,避免了網格劃分的前處理過程,并且繼承了數值算法中局部化的優點,在模擬復雜高維幾何域問題時有較好的準確性,并且保證計算精度的同時也提高了計算效率.本文探討了GFDM 在含阻抗邊界空腔聲學問題中的應用.以具有精確解的二維管道邊值問題為例,討論了GFDM 的穩定性和計算高效性,并獲得了最大計算頻率與節點間距關系的經驗公式.通過無精確解的復雜聲腔模型聲學分析,檢驗了GFDM 在二維和三維不同求解域內模擬計算的可行性、計算準確性和計算穩定性.

綜上所述,廣義有限差分法在空腔聲學問題中的應用已得到初步驗證,數值結果表明該方法是一種高效、精確、穩定、收斂的空腔聲學分析方法.然而,應用廣義有限差分法在更加實際的復雜聲場中還需要進一步研究,進行聲固耦合問題的求解和時域聲學仿真將是未來工作的重點.