光滑有限元的聲學(xué)研究：時域和頻域分析

何智成，李光耀，成艾國，鐘志華，周 澤

(1.湖南大學(xué) 汽車車身先進(jìn)設(shè)計制造國家重點(diǎn)實(shí)驗(yàn)室，長沙 410082;2.新加坡國立大學(xué) 機(jī)械工程系工程科學(xué)高級計算中心，新加坡 117576)

光滑有限元的聲學(xué)研究：時域和頻域分析

何智成1，2，李光耀1，成艾國1，鐘志華1，周 澤1

(1.湖南大學(xué) 汽車車身先進(jìn)設(shè)計制造國家重點(diǎn)實(shí)驗(yàn)室，長沙 410082;2.新加坡國立大學(xué) 機(jī)械工程系工程科學(xué)高級計算中心，新加坡 117576)

在使用有限元進(jìn)行聲場的數(shù)值模擬中，存在著兩個主要誤差，一個是數(shù)值方法中常規(guī)的插值誤差，另外一個是計算聲學(xué)中所特有的耗散誤差(dispersion error)，后者則是影響聲學(xué)模擬仿真置信度的最重要因素。產(chǎn)生耗散誤差的本質(zhì)原因是由于有限元的數(shù)值模型剛度“偏硬”造成的。為了控制耗散誤差，最重要的是使數(shù)值模型更好地反映真實(shí)模型。采用一種基于邊光滑的有限元方法(ES－FEM)來對聲場的時域和頻域進(jìn)行數(shù)值模擬研究。該方法只采用對復(fù)雜問題域適應(yīng)性很強(qiáng)的三角形網(wǎng)格，通過引進(jìn)基于邊的廣義梯度光滑技術(shù)，使得有限元系統(tǒng)得到適當(dāng)?shù)摹败浕薄r域和頻域的算例表明了在使用同樣網(wǎng)格的情況下，該方法在聲學(xué)模擬中的精度比有限元模型高。

邊光滑有限元;聲學(xué);時域;頻域

隨著人們對汽車乘員艙的聲學(xué)品質(zhì)要求越來越高，汽車的NVH性能研究已經(jīng)成了當(dāng)今的研究熱點(diǎn)。一般來說，汽車的振動噪聲分為低頻噪聲(1～200 Hz，主要為結(jié)構(gòu)噪聲)，中頻噪聲(200～500 Hz)以及高頻噪聲(＞200 Hz，主要為風(fēng)噪)［1］。有限元仍然是現(xiàn)在分析結(jié)構(gòu)噪聲或者聲場的主要方法，雖然在現(xiàn)有的CAE商業(yè)分析軟件中，有限元能夠分析各種各樣的聲學(xué)問題，但是由于有限元系統(tǒng)的剛度偏硬［2］，從而使得聲波在有限元數(shù)值模型中傳播的速度大于實(shí)際的聲速，從而導(dǎo)致了聲波的傳播模擬時存在很大的誤差。由于數(shù)值模型剛度的差別，在頻域方面則主要體現(xiàn)在有限元模擬產(chǎn)生一個除了常規(guī)插值誤差以外的耗散誤差(dispersion error)，這個誤差也是中高頻聲場計算的主要誤差［3］。因此有限元在聲波傳播問題的模擬技術(shù)以及中高頻聲學(xué)問題的仿真技術(shù)都還不太成熟。

為了提高聲學(xué)仿真模擬的精度，控制有限元的耗散誤差，最根本的方法是適當(dāng)軟化有限元系統(tǒng)的整體剛度，從而使得有限元聲學(xué)模型的剛度接近真實(shí)模型剛度。為了軟化有限元系統(tǒng)的剛度，一方面可以通過增加有限元網(wǎng)格的剖分密度，另一方面也可以通過提高有限元插值的階數(shù)［4－5］，使得有限元系統(tǒng)變軟，從而在計算聲場中得到令人滿意的結(jié)果。有些學(xué)者通過其他新的方法來控制計算聲場的耗散誤差，如無網(wǎng)格方法［6］，間斷有限元法［7］。但是這幾種方法實(shí)現(xiàn)比較復(fù)雜，而且這幾種方法都是通過增加計算時間來換取的比較好的精度。

為了在聲學(xué)數(shù)值模擬中得到較好的精度，本文采用了一種新型的光滑有限元技術(shù)。它采用光滑的迦遼金弱形式［8－9］來離散問題域，通過引入的梯度光滑技術(shù)，使得有限元系統(tǒng)得到適當(dāng)軟化，從而較好的模擬聲場。這種方法尤其能在三角形和四面體網(wǎng)格的計算中得到很好的精度，并且由于這兩種網(wǎng)格對任意復(fù)雜的模型具有很好的適應(yīng)性，所以這種方法能夠很方便的計算任何復(fù)雜模型。另外在效率上，這種方法在力學(xué)模型中的計算效率比有限元和無網(wǎng)格都要高［10］。因此本文采用這種光滑的有限元來技術(shù)對聲學(xué)的時域和頻域進(jìn)行研究，數(shù)值結(jié)果表明了這種新方法能夠在聲學(xué)模擬中得到很好的精度。

1 聲學(xué)計算的數(shù)學(xué)模型

1.1 聲學(xué)數(shù)值模擬的強(qiáng)形式

假定聲波傳播的介質(zhì)是理想流體，無粘滯性，在任意形狀的封閉聲場域Ω，聲場狀況可以用關(guān)于 p的波動方程描述：

其中：Δ為Laplace算子，c為聲波在介質(zhì)中傳播的速度，t為時間。在實(shí)際中，所計算的問題都是有邊界的，其邊界條件包含Dirichlet邊界條件、Neumann、Robin邊界，分別用ΓD、ΓN和ΓA表示，其邊界條件的控制方程可以表示為：

其中：ρ為聲音傳播介質(zhì)的密度，v為聲學(xué)流體粒子的運(yùn)動速度。

其中：Z為邊界的阻抗，一般來說Z=1/An，An為邊界的導(dǎo)納。

在線性聲學(xué)中，聲壓梯度與速度滿足如下的歐拉方程：

1.2 聲學(xué)計算的弱形式

聲場的弱形式可以通過光滑的迦遼金弱形式(Smoothed Galerkin weak form)對控制方程(1)進(jìn)行離散得到，光滑的迦遼金弱形式離散化的詳細(xì)過程可見文獻(xiàn)［9］。聲場的光滑迦遼金弱形式可以寫為：

引入基于邊的梯度光滑技術(shù)對協(xié)調(diào)梯度場進(jìn)行修正，并在光滑域內(nèi)進(jìn)行高斯積分，形成并組裝剛度矩陣，因此光滑有限元中的剛度矩陣不同于有限元的剛度矩陣，它是對于光滑域的積分，而不是對于單元的積分。

2 光滑有限元模型的建立

對于任意二維問題，首先用3節(jié)點(diǎn)的三角形網(wǎng)格(記為T3)對問題域進(jìn)行離散，得到了域內(nèi)的x個節(jié)點(diǎn)和N條邊，如圖1所示。則基于邊k的光滑域Ωk(k=1，2，…，N)的構(gòu)造過程為：如果邊在問題域的內(nèi)部，則圍繞邊k，用虛線分別聯(lián)接其邊的兩個端點(diǎn)和相鄰三角形網(wǎng)格的中心點(diǎn)，形成光滑有限元的基本單元即光滑域;如果三角形的邊在問題域的邊界上，則以全局邊界封閉光滑域。圖中n表示光滑域Ωk的邊界外法線向量。該光滑域作為ES－FEM形成和組裝剛度矩陣的基本單位，類似于有限元中的單元。ES－FEM的剛度矩陣是對這些光滑域Ωk進(jìn)行數(shù)值積分來形成和組裝剛度矩陣。通過類似的方法，光滑域的構(gòu)造可以拓展到三維的四面體網(wǎng)格中，在三維光滑域構(gòu)造中，光滑域模型的構(gòu)造基于四面體的每個面：通過連接相鄰四面體的中心及四面體的面的三個頂點(diǎn)，從而形成基于面光滑的光滑有限元的基本單元，具體細(xì)節(jié)可以參考文獻(xiàn)［9］。

圖1 基于三角形網(wǎng)格的二維積分光滑域Fig.1 The smoothing domain based on the background triangular mesh

在聲學(xué)中，采用廣義梯度光滑技術(shù)對聲場基本變量的梯度即速度(式(5))實(shí)施光滑操作，則光滑域Ωk的剛度矩陣可以寫為［9］：

由式(15)、式(16)可以看出，通過廣義梯度光滑技術(shù)，巧妙地將任意復(fù)雜的域內(nèi)積分轉(zhuǎn)換為簡單的邊界積分。對二維問題，則把面積積分轉(zhuǎn)換為線積分。因此，與傳統(tǒng)的有限元構(gòu)造的梯度場相比，廣義光滑梯度公式無需求解其形函數(shù)的導(dǎo)數(shù)，因而降低了對形函數(shù)連續(xù)可導(dǎo)的要求。

3 聲學(xué)模擬的時域和頻域模型

3.1 時域分析

在時域分析中，有許多的方法來模擬波的傳播問題，像Newmark方法，Crank－Nicholson方法。本文采用Newmark方法，即在t～t+Δt的時間區(qū)域內(nèi)，Newmark積分方法采用如下假設(shè)［11］：

其中：α和δ是按積分精度和穩(wěn)定性要求決定的參數(shù)，在這里α=0.25，δ=0.5，Newmark方法中的時間t+Δt的聲壓是通過滿足時間t+Δt的控制方程所得到的，所以有：

通過式(20)可以求解t+Δt時刻的聲壓，然后通過式(17)和(18)可以得到聲壓梯度(速度)。

3.2 頻域分析

在頻域中，聲壓及速度都可以表達(dá)成諧波的形式，因此p和v具有如下的表達(dá)式：

其中：Vn為邊界上的法向速度，對于頻率響應(yīng)分析，在給定問題域及邊界條件的情況下，都可以通過式(22)求出聲場某點(diǎn)的頻率響應(yīng)。

在剛性邊界條件下，自由模態(tài)的計算可以寫成如下的表達(dá)式：

其中：ωi為模態(tài)特征值，φi為模態(tài)特征向量。

4 數(shù)值算例

4.1 時域模型算例——聲波的傳播

波的傳播為例［13］，如圖(2)所示，在管的左端作用是Dirichlet邊界條件p=0，在管的右端作用v=H(t)階躍速度邊界條件，H(t)是Heaviside函數(shù)，即在小于或者等于零時，H(t)=0;在大于零時，H(t)=1。其他邊界都為剛性壁條件。為了便于計算，取c=1，介質(zhì)的密度為1.225 kg/m3。該模型被離散成的三角形網(wǎng)格的平均尺寸大小為0.05 m，具有207個節(jié)點(diǎn)和324個單元。該模型分別采用FEM和ES－FEM來計算聲波的傳播過程。取管子右端點(diǎn)A和中間點(diǎn)B(如圖2所示)的響應(yīng)來比較兩種算法的準(zhǔn)確度。

圖2 含Heaviside類型激勵的長方形管道Fig.2 Rectangle tube under Heaviside－type excitation

圖(3)、圖(4)分別為圖(2)中所示兩點(diǎn)的計算結(jié)果曲線。圖中都對有限元(FEM)，光滑有限元(ES－FEM)及理論解(Exact)的計算結(jié)果都進(jìn)行了標(biāo)明。通過比較發(fā)現(xiàn)：在節(jié)點(diǎn)和網(wǎng)格數(shù)相同的情況下，ES－FEM得到的聲壓結(jié)果比有限元的結(jié)果要精確得多，如圖3(a)，圖4(a)。在圖3(b)、圖4(b)中，隨著時間的變化，有限元計算結(jié)果的誤差不斷得到積累，越來越偏離理論解，而光滑有限元在整個時間歷程上一直和理論解都很吻合。為了進(jìn)一步說明方法的優(yōu)越性，本算例也采用有限元方法計算了兩種更密節(jié)點(diǎn)數(shù)的模型，計算模型和結(jié)果分別表示為FEM*(節(jié)點(diǎn)個數(shù)為729，單元個數(shù)為1 280)和FEM＊＊(節(jié)點(diǎn)個數(shù)為2 737，單元個數(shù)為5 120)。由圖中可以看出，隨著網(wǎng)格的加密，有限元的聲壓和速度計算結(jié)果越來越精確，從而越來越接近理論解;通過與基于粗網(wǎng)格(207個節(jié)點(diǎn)和324個單元)計算的光滑有限元結(jié)果比較，發(fā)現(xiàn)采用ES－FEM計算的結(jié)果好于FEM*的結(jié)果，甚至和FEM＊＊的計算結(jié)果差不多，誤差在整個時間歷程上一直都很接近理論解。由此可見ES－FEM在計算時域問題具有明顯的優(yōu)勢。這主要是由于基于邊光滑的有限元系統(tǒng)，能夠適當(dāng)?shù)剀浕^剛的有限元系統(tǒng)，從而更好地模擬真實(shí)的模型，提供比較精確的梯度解，因此與相同網(wǎng)格下的有限元計算結(jié)果相比較，光滑有限元得到的解比有限元的計算結(jié)果要精確得多。

圖3 A點(diǎn)響應(yīng)Fig.3 The response of point A

表1 計算成本和計算誤差的比較Tab.1 The comparison of the computational cost and accuracy

為了進(jìn)一步比較ES－FEM和FEM的計算成本，本文比較了在相同網(wǎng)格的前提下，光滑有限元和有限元的計算時間，另外同時也計算出了A，B點(diǎn)聲壓和速度在0～40 s整個時域內(nèi)的誤差，如表1所示。從表中可以看出，在相同的模型基礎(chǔ)上，ES－FEM的計算時間比FEM的計算時間稍長，但是ES－FEM得到的聲壓及速度結(jié)果比FEM的結(jié)果要精確得多。表1同時也列出了有限元采用兩種較密網(wǎng)格模型的計算時間及計算誤差。由表中可以看到，隨著有限元網(wǎng)格密度的加大，計算時間越長，并且計算精度越高。但是相對于光滑有限元而言，F(xiàn)EM得到相同精度所消耗的時間遠(yuǎn)遠(yuǎn)要大于ES－FEM消耗的時間。由此可見，光滑有限元的相對計算成本(同等時間獲得的精度)比有限元小。

4.2 頻域模型算例

一般來說，消聲器的內(nèi)部聲場比較復(fù)雜，平面波理論無法準(zhǔn)確地預(yù)測其消聲特性。為了計算消聲器的消聲性能，本文通過建立消聲器內(nèi)部聲場的二維數(shù)值模型，如圖5所示，并且對進(jìn)出口及壁面的邊界條件進(jìn)行合理的處理，來研究了消聲器的聲學(xué)模態(tài)和傳遞損失(Transmission Loss，TL)。介質(zhì)為空氣，密度為 ρ=1.225 kg/m2，聲速為c=340 m/s。該消聲器離散成294節(jié)點(diǎn)494個單元，分別采用有限元及光滑有限元對其進(jìn)行計算。由于此消聲器沒有理論解，理論證明［14－15］可以采用很密的網(wǎng)格，通過有限元求解得到問題的參考解。這里用來計算參考解的網(wǎng)格模型具有24 627個節(jié)點(diǎn)，48 273個四邊形單元。

圖4 B點(diǎn)響應(yīng)Fig.4 The response of point B

圖5 消聲器的二維圖Fig.5 The illustration of the two－dimensional muffler

假如消聲器的邊界及進(jìn)出口為剛性壁面，則消聲器的聲學(xué)模態(tài)可以通過有限元或者光滑的有限元求出。表2列出了有限元和光滑有限元求得的消聲器的前十階聲學(xué)模態(tài)，采用有限元計算的參考解也列在表2中。通過對消聲器模態(tài)的計算，發(fā)現(xiàn)有限元計算的模態(tài)特征值要比參考解要大，反映了有限元系統(tǒng)剛度偏硬的特性，而光滑有限元計算的模態(tài)特征值比有限元的結(jié)果小，反映了光滑有限元系統(tǒng)比有限元系統(tǒng)偏軟。數(shù)據(jù)結(jié)果也表明了有限元的誤差是光滑有限元誤差的3，4倍左右，說明了光滑有限元系統(tǒng)相對有限元系統(tǒng)來說得到了適當(dāng)?shù)能浕?，更適合求解聲學(xué)問題。

表2 消聲器的聲模態(tài)Tab.2 The eigenfrequencies of muffler

在計算消聲器傳遞損失時，設(shè)置的邊界條件要滿足標(biāo)準(zhǔn)的邊界條件定義，即在消聲器的左邊加一個單位的速度邊界條件，即vn=1 m/s。在消聲器的右邊加完全吸聲邊界條件，即An=1/(ρc)。對于進(jìn)出口面積相同的消聲器的傳遞損失可以由以下的方式計算［11］：

圖6中為同種網(wǎng)格模型情況下，有限元和光滑有限元計算的消聲器的傳遞損失曲線。為了便于比較，圖6也提供了有限元采用兩種不同網(wǎng)格模型的結(jié)算結(jié)果，分別表示 FEM*(節(jié)點(diǎn)個數(shù)為685，單元個數(shù)為1 222)和 FEM＊＊(節(jié)點(diǎn)個數(shù)為 1 519，單元個數(shù)為2 810)。另外在聲學(xué)數(shù)值計算中，采用非常精細(xì)的網(wǎng)格可以得到與真實(shí)值相差很小的參考解，這里采用有限元計算參考解的模型含有24 627個節(jié)點(diǎn)，48 273個四邊形單元。通過比較，可以看出，在低頻段，有限元，光滑的有限元模型和參考解差不多，都十分接近參考解，但是隨著頻率的升高，ES－FEM的結(jié)果比FEM的結(jié)果更接近參考解，甚至達(dá)到了和FEM＊＊同等的精度。通過該算例，說明在相同的離散模型中，光滑有限元不論在低頻或高頻，數(shù)值結(jié)果都比有限元好。

圖6 ES－FEM和FEM計算的消聲器的傳遞損失Fig.6 The transmission loss of the muffler using the present ES－FEM and FEM

5 結(jié)論

本文對基于邊光滑的有限元在聲學(xué)的時域和頻域進(jìn)行了研究，并且通過數(shù)值算例對該方法進(jìn)行了驗(yàn)證，得到了如下的結(jié)論：

(1)該方法采用三角形單元，對任何復(fù)雜的問題具有很好的適應(yīng)性，所以該方法適用于任何復(fù)雜的模型。

(2)該方法引入了梯度光滑的理論，能夠適當(dāng)軟化過剛的有限元系統(tǒng)，從而得到比較精確的梯度場，因此比有限元更適合計算聲波的時域問題。

(3)該方法計算的模態(tài)特征值要小于有限元得到的值，并且更接近參考解的值，因此該方法的數(shù)值模型系統(tǒng)剛度比有限元的軟，因而在頻域問題中，該方法比有限元具有更好的精度。

［1］龐 劍，諶 剛，何 華.汽車噪音與振動——理論與應(yīng)用［M］.北京：北京理工大學(xué)出版社，2006.

［2］Liu G R.Meshfree methods：Moving beyond the finite element method［M］.USA：CRC Press，2009.

［3］Deraemaeker A，Babuska I，Bouillard Ph.Dispersion and pollution of the FEM solution for the helmholtz equation in one，two and three dimension［J］.Int.J.Numer.Meth.Engrg，1999，46：471－499.

［4］Petersen S，Dreyer D，Estorff O V.Assessment of finite and spectralelementshape functions or efficient iterative simulations of interior acoustics［J］.Comput.Methods Appl.Mech.Engrg，2006，195：6463－6478.

［5］Biermann J，Estorff O V，Petersen S，et al.Higher order finite and infinite elements for the solution of Helmholtz problems［J］.Comput.Methods Appl.Mech.Engrg，2009，198：1171－1188.

［6］ Bouillard Ph，Suleau S.Element－free Garlekin solutions for Helmholtz problems：formulation and numerical assessment of the pollution effect［J］.Comput.Methods Appl.Mech.Engrg.，1998，162：317－335.

［7］Alvarez G B，Loula A F D，Carmo E G D D，et al.A discontinuous finite element formulation for Helmholtz eqation［J］.Comput.Methods Appl.Mech.Engrg，2006，195：4018－4035.

［8］Liu G R.A weakened weak(W2)form for a unified formulation of compatible and incompatible methods，part I：theory and part II：applications to solid mechanics problems［J］.Int.J.Numer.Meth.Engrg，2010，81：1093 －1126.

［9］ He Z C，Liu G R，Zhong Z H，et al.An edge－based smoothed finite element method(ES－FEM)for analyzing three－dimensional acoustic problems［J］.Comput Methods Appl Mech Eng，2009，199：20－33.

［10］ Chen L，Nguyen H X，Nguyen T T，et al.Assessment of smoothed point interpolation methods for elastic mechanics［J］.Commun.Numer.Meth.Engng，2009，DOI：10.1002/cnm，1251.

［11］ Zienkiewicz O C，Taylor R L.The finite element method［M］.Butterworth－Heinemann：Oxford，2000.

［12］ Munjal M L.Acoustics of ducts and mufflers with application to exhaust and ventilation system design［M］.New York：Wiley，1987.

［13］ Wu T W.Boundary element in acoustics：fundamentals and computer codes［M］.Southampton：WIT Press，2000.

［14］ Ihlenburg F.Finite element analysis of acoustic scattering［M］.Springer Press，1998.

［15］ Ihlenburg F，Babuska I.Finite element solution of the Helmholtz equation with high wave－number part I：the hversion of the FEM［J］.Comput.Math.Appl，1995，30(9)：9－37.

Acoustic analysis using edge-based smoothed finite element method：time and frequency domain analysis

HE Zhi－cheng1，2，LI Guang－yao1，CHENG Ai－guo1，ZHONG Zhi－h(huán)ua1，ZHOU Ze1
(1.State Key Lab.of Advanced Technology for Vehicle Body Design＆ Manufacture，Hunan University，Changsha 410082，China;2.Centre for Advanced Computations in Engineering Science(ACES)，Department of Mechanical Engineering，National University of Singapore，117576 Singapore)

The standard finite element method(FEM)encounters two kinds of errors in the computational acoustic problems：the conventional interpolation error and the dispersion error，the later of which plays an important role in the computational acoustics.The dispersion error is rooted at the"overly－stiff"property of the FEM.In order to control the dispersion error，the most effective measure is to make the numerical model reflect the exact model.An edge－based smoothed finite element method(ES－FEM)was adopted for acoustic problems both in frequency domain and time domain.By using only triangular mesh which is very adaptive for any complicated domain and introducing gradient smoothing operation performed over each edge－based smoothed domain，the present method can reduce the stiffness and provide proper“softness”to the model.The results demonstrate that the ES－FEM can provide more accurate solution both in time and frequency domain compared with the linear FEM using the same meshes.

edge－based smoothed element method(ES－FEM);acoustics;time domain;frequency domain

O42

A

教育部博士學(xué)術(shù)新人獎項(xiàng)目資助;教育部長江學(xué)者與創(chuàng)新團(tuán)隊(duì)發(fā)展計劃(531105050037)

2011－04－26 修改稿收到日期：2011－08－23

何智成 男，博士生，1983年生

成艾國 男，教授，1972年生