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

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

(1.湖南大學 汽車車身先進設計制造國家重點實驗室，長沙 410082;2.新加坡國立大學 機械工程系工程科學高級計算中心，新加坡 117576)

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

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

(1.湖南大學 汽車車身先進設計制造國家重點實驗室，長沙 410082;2.新加坡國立大學 機械工程系工程科學高級計算中心，新加坡 117576)

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

邊光滑有限元;聲學;時域;頻域

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

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

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

1 聲學計算的數學模型

1.1 聲學數值模擬的強形式

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

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

其中：ρ為聲音傳播介質的密度，v為聲學流體粒子的運動速度。

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

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

1.2 聲學計算的弱形式

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

引入基于邊的梯度光滑技術對協調梯度場進行修正，并在光滑域內進行高斯積分，形成并組裝剛度矩陣，因此光滑有限元中的剛度矩陣不同于有限元的剛度矩陣，它是對于光滑域的積分，而不是對于單元的積分。

2 光滑有限元模型的建立

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

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

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

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

3 聲學模擬的時域和頻域模型

3.1 時域分析

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

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

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

3.2 頻域分析

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

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

在剛性邊界條件下，自由模態的計算可以寫成如下的表達式：

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

4 數值算例

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

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

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

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

圖3 A點響應Fig.3 The response of point A

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

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

4.2 頻域模型算例

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

圖4 B點響應Fig.4 The response of point B

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

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

表2 消聲器的聲模態Tab.2 The eigenfrequencies of muffler

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

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

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

5 結論

本文對基于邊光滑的有限元在聲學的時域和頻域進行了研究，并且通過數值算例對該方法進行了驗證，得到了如下的結論：

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

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

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

［1］龐 劍，諶 剛，何 華.汽車噪音與振動——理論與應用［M］.北京：北京理工大學出版社，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－hua1，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

教育部博士學術新人獎項目資助;教育部長江學者與創新團隊發展計劃(531105050037)

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

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

成艾國 男，教授，1972年生