無網格FEM-BEM法計算艙室空間聲傳遞函數

劉延善， 曾向陽， 王海濤

(西北工業大學 航海學院，西安　710072)

無網格FEM-BEM法計算艙室空間聲傳遞函數

劉延善， 曾向陽， 王海濤

(西北工業大學 航海學院，西安710072)

摘要：對于艙室空間低頻聲場計算問題，通常采用有限元法(FEM)或者邊界元法(BEM)，但是當封閉空間內部存在其他物體，又需要考慮外部激勵對內場的影響時，直接采用單一的前述方法不能滿足建模計算要求，而且這些基于網格的方法在前處理、光順化處理等方面都存在不足。針對這些問題，將有限元和邊界元法結合，并將其無網格化，這對于封閉環境的聲場預測是一種新的嘗試。首先推導了無網格FEM-BEM計算模型，對微分方程的離散、形函數構造等細節進行了詳細闡述，然后以實際算例對提出的方法進行了驗證，結果表明，無網格FEM-BEM方法和SYSNOISE的計算結果保持一致；進一步與實測結果的對比表明，該方法在低頻率范圍內，各測點平均聲壓級相對誤差在5.26%以內。這說明，該方法不僅適用于復雜問題的計算，同時還具有良好的計算精度。

關鍵詞：有限元法；邊界元法；無網格法；聲傳遞函數

對于飛機艙室空間低頻聲場的數值計算，現有的研究多集中于有限元法(FEM)、邊界元(BEM)法等[1]。對于簡化的封閉空間，單一的有限元法或者邊界元法可以滿足計算要求[2-3]，但是當封閉空間內部存在座椅等物體時，一般只能采用有限元法進行分析[4]，而如果此時還需要考慮聲空間外部的擾動對內部的聲場影響時，則單一方法已無法滿足需求[5]，需要結合有限元法和邊界元法來進行計算。

無論是有限元、邊界元，還是二者結合的方法，有效計算頻率都受到網格密度的影響，對于復雜結構，還可能出現網格畸變問題。筆者所在課題組將力學領域的無網格法引入聲學[6]，建立的聲學無網格模型可以較好地克服這方面的問題。該方法去除了定義在求解域上的網格單元，前處理只要節點位置信息，可以方便地在求解域內布置節點，從而可以在局部求解問題上提升計算精度。不過，現有模型只考慮了聲源位于艙內的情況[7]，且未考慮艙內存在其他物體的情況。

實際艙室聲場研究中更關心外部聲源對于復雜艙室結構內部聲場的作用，因此，本文首先研究基于傳統聲學有限元法和聲學邊界元法相結合的準則來進行聲場數值計算，在此基礎上，將混合模型無網格化，提出一種無網格FEM-BEM模型。然后以一個矩形結構和一個圓柱結構為例，對其計算準確性進行初步檢驗。最后對一個小型艙段結構進行實驗測試，進一步驗證該模型的有效性。

1無網格FEM-BEM模型原理

在工程實際中，數值模型的計算本質上是偏微分方程的求解。首先針對描述問題的方程，通過加權余量法或者變分法將其轉化為等效積分形式。其次通過選擇合適的離散方式并施加相應的物理邊界條件，可以將數值積分問題轉化為線性方程組的求解。最后通過分配滿足計算要求的積分點并選擇合適的線性方程組求解方法，可得到各離散點處的聲學參量信息。

本文中的問題是聲激勵處于結構外部，需要求解艙室結構內部的聲場，這需要借鑒傳統的FEM-BEM結合方法。圖1為本文混合方法的原理示意圖。首先，無網格化BEM用于計算結構外部聲場，得到內部聲場的邊界信息，具體理論實現在后續第三節說明。其次，利用前一步所得邊界信息，借助無網格化的FEM法可獲得結構內部的任意位置的聲場信息，具體理論實現在后續第二節說明。基于這種結合思路，可以求解結構在外部激勵作用下的速度邊界，進而可以得到內部任意位置的聲學響應，由此可進一步求解其他的聲學參量。為說明無網格混合方法，首先闡述無網格法的基本原理，其次分別對無網格FEM和無網格BEM算法實現進行說明。

圖1　無網格FEM-BEM原理示意圖Fig.1 Schematic diagram of Meshless FEM-BEM

1.1無網格法

無網格法的基本思想是將待求解問題的區域離散為多個有限數目的節點，采用節點的權函數或者核函數來表征節點及其鄰域內的聲學參量，即利用權函數來近似地表示影響域內的聲場函數，進而形成與節點聲場相關的系統方程，之后再離散求解。

無網格法目前有十余種，區別體現在所用的微分方程的等效積分形式和離散近似方案。基于微分方程等效積分的加權余量法是求解偏微分方程的一種十分有效的方法，無網格法也可以用加權余量法或變分法來建立整體求解的系統方程[8]。加權權函數選取的不同可以給出不同的加權余量方法。常見的權函數有配點法、最小二乘法、Galerkin法等。后續計算主要采用Galerkin加權余量法來構造系統方程[9]。

形函數[10]直接決定了微分方程的離散方式。與一般的方法不同，在無網格法中，形函數也稱試探函數，它是定義在局部區域(支撐域)中的函數，只在其支撐域中有定義，而在其支撐域外為零，所以也稱緊支試探函數。在二維問題中，支撐域一般取為圓形域或矩形域。如圖2所示，在求解域Ω上，排布若干節點，各節點具有圓或矩形等幾何外形構成的權函數影響域ΩI，節點處的聲壓大小可通過權函數對求解域中任意點的聲參量作用產生不同程度的貢獻。

圖2　節點的圓形支撐域和矩形支撐域Fig.2 Nodal compacted supported domain

權函數是緊支函數，它只在對于節點xI生成的w(x-xI)≠0的局部區域有定義，而在其他區域為零。權函數有定義的局部區域一般稱為緊支撐域。權函數與空間距離有關系，設

w(x-xI)=w(d)

(1)

在問題域內，某位置的待求聲學參量p(x)可以由一組離散節點xI(I=1,2,…,n)與相對應的緊支函數NI(x)的線性組合近似表示為

(2)

式中，n是求解域中節點的總數，pI是函數p(x)在節點xI處的值，即p(xI)=pI。N和p分別為

N=[N1(x),N2(x),…,NN(x)]

(3)

p=[p1,p2,…,pN]T

(4)

對于任意位置處待求參量的求解，無網格法中需要對所有節點的貢獻量進行一次計算，而不像有限元法只是疊加該位置所包圍單元各節點的聲學參量的貢獻，從這一點看，無網格比FEM法計算量略多。但是在求解系統方程不同矩陣時，無網格不需要進行單元矩陣的組裝過程。

本文采用移動最小二乘法近似來構造形函數[11]。利用上述原理，可以分別將有限元算法和邊界元算法與無網格法結合，最后再實現二者的融合，下面分別敘述相應的理論實現。

1.2無網格有限元模型

如圖3所示，假設一封閉空間體積為V，內壁總面積為S，其中S1為剛性邊界面，S2為吸聲邊界面，S3為振動速度邊界面，n為外法向導數。

圖3　封閉空間示意圖Fig.3 Schematic of the enclosed space

邊界條件分別為：

S1邊界上：

(5)

S2邊界上：

(6)

S3邊界上：

(7)

封閉聲場內聲壓的計算屬于非齊次邊界條件的二階偏微分方程的求解，使用Galerkin法，經過推導[6]可得到系統的有限元方程為：

(K+jρ0ωC-k2M)p=F

(8)

式中,K為剛度矩陣，C為阻尼矩陣，M為質量矩陣，F為載荷矩陣。對模型使用無網格法離散化，根據上文的理論可實現形函數的計算，代入系統方程后，各個矩陣計算形式為：

(9)

(10)

(11)

(12)

式中，n為模型總的域內積分點數目，l為代表阻抗邊界的總積分點數，m為速度邊界面所包含的積分點數，ωi為體積分和面積分離散化時的各積分點的權系數，N為無網格法構造的形函數[11]，在給定幾何模型離散節點分布時，可直接求解系統方程中的不同矩陣。而前兩種邊界條件式(5)、式(6)可以通過阻尼矩陣C來添加，而速度邊界條件式(7)則通過載荷向量F來實現，但這里的速度需要借助無網格BEM來求解，最后求解式(8)就可獲得封閉空間內部任意點在給定邊界條件和速度邊界條件下的聲傳遞函數。

當在實際中考慮內部空間座椅等問題時，通常認為聲空間中包含座椅體積所占據的空間，聲學模型中包含了座椅體積所占據的空間，這里是通過定義座椅空間部分為不同的流體材料來進行模擬計算的[4]。

1.3無網格邊界元模型

對于外部聲場，參考經典的Kirchhoff-Helmholtz方程

(13)

將問題邊界采用無網格的方式進行離散，遍歷所有節點(上文的P)，整理各離散公式后，可得到系統的線性方程組[12]：

Cp=GG·p-GH·p+F

(14)

式中，C為立體角矩陣，GG和GH為各積分點對不同節點的影響矩陣，F為表征聲源影響的載荷矩陣，p為聲壓列向量，pI表示所有積分點對節點I的聲壓貢獻(I=1,2,3,…,n)。它們的表示形式為：

(15)

(16)

(17)

(18)

(19)

F(ω)=jρωq(r0)G(P,r0)

(20)

(21)

式(17)、(19)為影響矩陣各元素計算公式，i，j分別表示節點序號和計算點序號，當遍歷所有節點和計算點時，可得影響矩陣，再通過求解式(14)即可得到聲場邊界曲面上任意節點處的聲壓值，再通過法向阻抗與聲壓振速的關系即可得到邊界節點上的振速及其他聲學邊界條件。

1.4無網格FEM-BEM算法

為更清楚地說明本文方法的實現過程，圖4給出了本文方法的計算流程圖。

圖4　混合模型計算流程圖Fig.4 Hybrid computing flowchart

對于Meshless BEM，在獲得離散無網格模型后，遍歷所有節點和積分點，計算式(14)中的各影響矩陣，同時設定相應的邊界條件，最后求解式(14)可獲得結構邊界的振速信息。而對于Meshless FEM，經無網格離散化，遍歷所有節點和積分點，計算式(8)中的不同矩陣，同時將施加所得的邊界速度，最后求解式(8)所形成的線性方程組，可得到空間內任意位置的聲場分布。

2數值算例及結果分析

2.1算例1

為了便于與商用軟件對比，首先以一個圓柱殼結構(內部全空)的內部聲場分析為例，如圖5所示。源點位置(3,0,0)，體積速度取1。接收點位置從左到右分別為(0,0,0.9)、(0,0,1.5)、(0,0,2.1)，坐標原點固定在圖中左端板的中心點處。圓柱殼半徑為1 m，軸向長度為3 m。吸聲條件是在所有壁面上具有阻尼邊界條件，內壁聲阻抗為40ρ0c0，外壁聲阻抗取ρ0c0，其中ρ0=1.21 kg/m3，c0=344 m/s。聲源的體積速度取為1，源信號為聲脈沖信號[6]。

圖5　圓柱殼結構示意圖Fig.5 Schematic of the cylindrical shell structure

圖6為無網格模型的節點分布示意圖。取支撐域半徑為1.5 m。根據所建立的無網格理論模型，計算了50 Hz～250 Hz頻段內聲源點到接收點的聲傳遞函數，同時在SYSNOISE中進行相應的數值計算，計算的頻帶范圍保持一致。圖7為相應的各場點結果。

圖6　圓柱形結構內部及截面無網格節點分布Fig.6 Nodal distribution of the meshless model

圖7　各場點聲傳遞函數Fig.7 Sound transfer functions at different field-points

對于場點p1、p2、p3，在特定的頻率上都有極大峰值存在，如58 Hz、123 Hz、209 Hz、240 Hz，三個場點的平均聲壓級分別為79.1 dB、85.4 dB、79.4 dB。為驗證本方法的準確性，將計算結果與商用軟件SYSNOISE進行對比分析。圖8為場點2的聲傳遞函數，“MMFEMBEM”表示無網格FEM-BEM模型計算結果，“SYSNOISE”表示SYSNOISE軟件得到的結果。

圖8　p2場點幅值頻率響應Fig.8 frequency response of the magnitude at p2

由圖8可知，在250 Hz以下，兩種方法獲得的聲壓變化趨勢非常接近。SYSNOISE得到的峰值頻率位置為114 Hz、208 Hz和236 Hz，無網格法得到的峰值頻率為114 Hz、208 Hz和239 Hz，二者得到的峰值頻率保持一致，這些頻率點也正與圓柱聲腔主要的模態相對應。從平均聲壓級來看，無網格FEM-BEM模型的結果為88.9 dB，SYSNOISE的結果為92.6 dB，偏差為3.7 dB，相對于SYSNOISE的相對誤差為4.0%。這些結果初步驗證了無網格FEM-BEM模型的正確性。

2.2算例2

為驗證本方法適用于艙室空間內部存在其他物體的情況，這里以矩形封閉空間內部放置矩形塊為例進行仿真計算。選取的矩形封閉空間長1.0 m，寬1.1 m，高1.2 m。外壁面為剛性壁面，y=0的內壁面為吸聲壁面，阻抗大小取10ρ0c0。聲源為點聲源，體積速度取1，坐標為(3,0,0)。測點坐標為(0.8,0.6,0.9)。矩形塊底面中心與空間底面中心重合，長寬高分別為0.5 m、0.55 m、0.24 m，密度為22 kg/m3，其中聲速為70 m/s，矩形封閉空間形狀、聲源點和測點的分布情況，如圖9所示。

圖9　矩形空間示意圖Fig.9 Schematic of rectangular space

分別采用有網格FEM-BEM方法和無網格的FEM-BEM方法對上述模型進行計算。FEM-BEM的單元數為125，節點數為216，MMFEMBEM使用的節點數為216，積分點數為1 000，計算的頻段均為0～300 Hz。圖10為使用兩種不同方法計算得到的聲傳遞函數。

由圖10可知，有網格FEM-BEM法和無網格FEM-BEM法所得的聲傳遞函數結果很接近。從曲線整體走勢來分析，在低頻段二者的聲傳遞函數變化趨勢是一致的，在250～300 Hz范圍內，兩條曲線的差異有逐漸加大的趨勢。分析原因，算例中所用模型的Schroeder頻率[13]為302 Hz，根據文獻理論可知，在Schroeder頻率附近，是模態由稀疏變密集的過渡頻帶，嚴格地講，在頻率大于302 Hz時，波動方法已不再適用于求解該問題。以局部峰值頻率的位置進行對比，在144 Hz、173 Hz、225 Hz、274 Hz處，兩種方法得到的相應峰值是對應的，這幾個點也正與該矩形空間的前幾階模態對應。總結起來，在低頻率時，兩種算法得到聲傳遞函數基本保持一致，局部峰值頻率點與相應的聲模態階次保持一致。

圖10　聲傳遞函數對比Fig.10 Comparison of transfer functions using different methods

3實驗驗證

為了進一步驗證無網格有限元-邊界元模型的正確性，在半消聲室中對一個實際艙段進行了測量實驗。測量系統如圖11所示，圖12為艙段模型。圖13為傳聲器位置分布示意圖。

圖11　實驗測量系統Fig.11 Experimental measuring principle

圖12　艙段模型Fig.12 Cabin model

圖13　傳聲器位置截面分布示意圖Fig.13 Distribution of microphones and source

以測點2為例。圖14為數值計算和實驗測試的1/3倍頻帶聲壓級結果的對比圖。圖15為測點結果的誤差曲線。“fp2”表示測點2。

由圖可知，在0～200 Hz時，數值計算結果和實驗測量結果偏差在5 dB以內，而在200～500 Hz時，實驗結果和數值仿真結果的誤差大于5 dB，250 Hz時，聲壓級的誤差最大。根據整體的聲壓級變化趨勢，測點2的數值計算結果和實驗測量結果基本保持一致。

圖14　測點2聲壓級結果對比Fig.14 SPL comparison at fp2

圖15　測點2聲壓級誤差曲線Fig.15 SPL error at fp2

為計算各測點的平均聲壓級，將實驗測得的1/3倍頻帶聲壓級結果進行對數平均，同時將數值仿真計算得到的各測點的聲傳遞函數結果也作對數平均。表1為各測點的實驗測量和數值計算總聲壓級結果和相對誤差。

表1　數值計算與測量結果對比

根據表1中數據可求得7個測點的平均相對誤差為5.26%，說明無網格FEM-BEM法與實際測試結果的總體偏差很小，這表明無網格FEM-BEM計算模型是正確的，可以用于外部聲源作用下的封閉結構內部的聲場計算問題。

4結論

本文提出一種無網格化的有限元和邊界元相結合的方法來計算聲傳遞函數，可用于復雜條件下的艙室內部聲場分布預測。仿真算例和實測結果都表明，提出的方法在低頻范圍內，可以取得較為準確的預測效果，因而具有比單一的有限元或邊界元法更寬的適用范圍。在后面的研究工作，還將進一步分析本方法提升精度和計算效率方面的工作。

參 考 文 獻

[1] Marburg S,Nolte B,Bernhard R J，et al. Computational acoustics of noise propagation in fluids:finite and boundary element methods[M]. Springer,2008.

[2] Ihlenburg F. Finite element analysis of acoustic scattering[M]. Springer,1998.

[3] Otsuru T,Tomiku R. Basic characteristics and accuracy of acoustic element using spline function in finite element sound field analysis[J]. Acoustical Science and Technology,2000,21(2):87-95.

[4] 李增剛,詹福良. Virtual. Lab Acoustics 聲學仿真計算高級應用實例[M]. 北京：國防工業出版社,2010.

[5] Citarella R,Federico L,Cicatiello A. Modal acoustic transfer vector approach in a FEM-BEM vibro-acoustic analysis[J]. Engineering Analysis with Boundary Elements,2007,31(3):248-258.

[6] Wang H T,Zeng X Y,Chen L. Calculation of sound fields in small enclosures using a meshless model[J]. Applied Acoustics,2013,74(3):459-466.

[7] 王海濤,曾向陽,陳玲. 小尺度封閉空間無網格Galerkin聲場數值計算方法[J]. 西北工業大學學報,2012,30(1):102-107.

WANG Hai-tao,ZENG Xiang-yang,CHEN Ling.An effective galerkin meshless model for calculating sound fields in arbitrarily shaped enclosures[J].Journal of Northwestern Polytechnical University,2012 30(1):102-107.

[8] Craggs A. Acoustic modeling:finite element method[J]. Handbook of Acoustics,MJ Crocker,Ed. New York:Wiley,1998:149-156.

[9] 王海濤,曾向陽. 周期結構聲散射系數的無網格數值計算方法[J]. 計算物理,2013,30(2):229-236.

WANG Hai-tao,ZENG Xiang-yang. Meshless models for numerical calculation of scattering coefficients of periodic structures[J].Chinese Journal of Computational Physics,2013,30(2):229-236.

[10] 張雄,劉巖,馬上. 無網格法的理論及應用[J]. 力學進展,2009,39(1):1-36.

ZHANG Xiong,LIU Yan,MA Shang.Meshfree methods and their applications[J].Advances in Mechanics,2009,39(1):1-36.

[11] 張雄,劉巖. 無網格法[M]. 北京:清華大學出版社,2004.

[12] Gunda R. Boundary element acoustics and the fast multipole method (FMM)[J]. Sound and Vibration,2008,42(3):12-16.

[13] Schroeder M R. The “Schroeder frequency” revisited[J]. The Journal of the Acoustical Society of America,1996,99(5):3240-3241.

Meshless FEM-BEM method for computing sound transfer function in enclosed spaces

LIU Yan-shan, ZENG Xiang-yang, WANG Hai-tao

(School of Marine Science and Technology, Northwestern Polytechnic University, Xi’an 710072, China)

Abstract:The finite element method(FEM) or boundary element method(BEM) is usually employed in sound field calculation of a cabin space at lower frequency, but a single method can not satisfy computing requirements of complex structures. These methods have some disadvantages, such as, lack of pre-processing and smoothing processing. Aiming at these problems, here, FEM and BEM were combined into a hybrid algorithm, and then transform it was converted into a meshless algorithm. It was a new approach for a sound prediction problem in enclosed environment. At first, the combined meshless FEM-BEM model was derived, and then the specific details including discretization of differential equations, construction of shape functions, nodal disposal, and were combined integral computing scheme were described. Finally, numerical simulations and tests were conducted to verify this method. The simulation results showed that the results using the proposed method agree well with those using SYSNOISE. The test results showed that the average relative error of each position’s sound pressure of a cabin space is less than 5.26%. These results showed that the proposed method not only is applicable for complex problems, but also has a good computation accuracy.

Key words:finite element method (FEM); boundary element method (BEM); meshless method; sound transfer function

基金項目：國家自然科學基金(11374241)；航空科學基金(20151553021)

收稿日期：2015-05-05修改稿收到日期：2015-11-09

通信作者曾向陽 男，博士，教授，1974年4月生

中圖分類號:O422

文獻標志碼：A

DOI:10.13465/j.cnki.jvs.2016.09.002

第一作者 劉延善 男，博士生，1987年4月生