劉延善， 曾向陽， 王海濤

(西北工業(yè)大學(xué) 航海學(xué)院，西安　710072)

無網(wǎng)格FEM-BEM法計(jì)算艙室空間聲傳遞函數(shù)

劉延善， 曾向陽， 王海濤

(西北工業(yè)大學(xué) 航海學(xué)院，西安710072)

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

關(guān)鍵詞：有限元法；邊界元法；無網(wǎng)格法；聲傳遞函數(shù)

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

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

實(shí)際艙室聲場研究中更關(guān)心外部聲源對(duì)于復(fù)雜艙室結(jié)構(gòu)內(nèi)部聲場的作用，因此，本文首先研究基于傳統(tǒng)聲學(xué)有限元法和聲學(xué)邊界元法相結(jié)合的準(zhǔn)則來進(jìn)行聲場數(shù)值計(jì)算，在此基礎(chǔ)上，將混合模型無網(wǎng)格化，提出一種無網(wǎng)格FEM-BEM模型。然后以一個(gè)矩形結(jié)構(gòu)和一個(gè)圓柱結(jié)構(gòu)為例，對(duì)其計(jì)算準(zhǔn)確性進(jìn)行初步檢驗(yàn)。最后對(duì)一個(gè)小型艙段結(jié)構(gòu)進(jìn)行實(shí)驗(yàn)測試，進(jìn)一步驗(yàn)證該模型的有效性。

1無網(wǎng)格FEM-BEM模型原理

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

本文中的問題是聲激勵(lì)處于結(jié)構(gòu)外部，需要求解艙室結(jié)構(gòu)內(nèi)部的聲場，這需要借鑒傳統(tǒng)的FEM-BEM結(jié)合方法。圖1為本文混合方法的原理示意圖。首先，無網(wǎng)格化BEM用于計(jì)算結(jié)構(gòu)外部聲場，得到內(nèi)部聲場的邊界信息，具體理論實(shí)現(xiàn)在后續(xù)第三節(jié)說明。其次，利用前一步所得邊界信息，借助無網(wǎng)格化的FEM法可獲得結(jié)構(gòu)內(nèi)部的任意位置的聲場信息，具體理論實(shí)現(xiàn)在后續(xù)第二節(jié)說明?；谶@種結(jié)合思路，可以求解結(jié)構(gòu)在外部激勵(lì)作用下的速度邊界，進(jìn)而可以得到內(nèi)部任意位置的聲學(xué)響應(yīng)，由此可進(jìn)一步求解其他的聲學(xué)參量。為說明無網(wǎng)格混合方法，首先闡述無網(wǎng)格法的基本原理，其次分別對(duì)無網(wǎng)格FEM和無網(wǎng)格BEM算法實(shí)現(xiàn)進(jìn)行說明。

圖1　無網(wǎng)格FEM-BEM原理示意圖Fig.1 Schematic diagram of Meshless FEM-BEM

1.1無網(wǎng)格法

無網(wǎng)格法的基本思想是將待求解問題的區(qū)域離散為多個(gè)有限數(shù)目的節(jié)點(diǎn)，采用節(jié)點(diǎn)的權(quán)函數(shù)或者核函數(shù)來表征節(jié)點(diǎn)及其鄰域內(nèi)的聲學(xué)參量，即利用權(quán)函數(shù)來近似地表示影響域內(nèi)的聲場函數(shù)，進(jìn)而形成與節(jié)點(diǎn)聲場相關(guān)的系統(tǒng)方程，之后再離散求解。

無網(wǎng)格法目前有十余種，區(qū)別體現(xiàn)在所用的微分方程的等效積分形式和離散近似方案?；谖⒎址匠痰刃Хe分的加權(quán)余量法是求解偏微分方程的一種十分有效的方法，無網(wǎng)格法也可以用加權(quán)余量法或變分法來建立整體求解的系統(tǒng)方程[8]。加權(quán)權(quán)函數(shù)選取的不同可以給出不同的加權(quán)余量方法。常見的權(quán)函數(shù)有配點(diǎn)法、最小二乘法、Galerkin法等。后續(xù)計(jì)算主要采用Galerkin加權(quán)余量法來構(gòu)造系統(tǒng)方程[9]。

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

圖2　節(jié)點(diǎn)的圓形支撐域和矩形支撐域Fig.2 Nodal compacted supported domain

權(quán)函數(shù)是緊支函數(shù)，它只在對(duì)于節(jié)點(diǎn)xI生成的w(x-xI)≠0的局部區(qū)域有定義，而在其他區(qū)域?yàn)榱?。?quán)函數(shù)有定義的局部區(qū)域一般稱為緊支撐域。權(quán)函數(shù)與空間距離有關(guān)系，設(shè)

w(x-xI)=w(d)

(1)

在問題域內(nèi)，某位置的待求聲學(xué)參量p(x)可以由一組離散節(jié)點(diǎn)xI(I=1,2,…,n)與相對(duì)應(yīng)的緊支函數(shù)NI(x)的線性組合近似表示為

(2)

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

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

(3)

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

(4)

對(duì)于任意位置處待求參量的求解，無網(wǎng)格法中需要對(duì)所有節(jié)點(diǎn)的貢獻(xiàn)量進(jìn)行一次計(jì)算，而不像有限元法只是疊加該位置所包圍單元各節(jié)點(diǎn)的聲學(xué)參量的貢獻(xiàn)，從這一點(diǎn)看，無網(wǎng)格比FEM法計(jì)算量略多。但是在求解系統(tǒng)方程不同矩陣時(shí)，無網(wǎng)格不需要進(jìn)行單元矩陣的組裝過程。

本文采用移動(dòng)最小二乘法近似來構(gòu)造形函數(shù)[11]。利用上述原理，可以分別將有限元算法和邊界元算法與無網(wǎng)格法結(jié)合，最后再實(shí)現(xiàn)二者的融合，下面分別敘述相應(yīng)的理論實(shí)現(xiàn)。

1.2無網(wǎng)格有限元模型

如圖3所示，假設(shè)一封閉空間體積為V，內(nèi)壁總面積為S，其中S1為剛性邊界面，S2為吸聲邊界面，S3為振動(dòng)速度邊界面，n為外法向?qū)?shù)。

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

邊界條件分別為：

S1邊界上：

(5)

S2邊界上：

(6)

S3邊界上：

(7)

封閉聲場內(nèi)聲壓的計(jì)算屬于非齊次邊界條件的二階偏微分方程的求解，使用Galerkin法，經(jīng)過推導(dǎo)[6]可得到系統(tǒng)的有限元方程為：

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

(8)

式中,K為剛度矩陣，C為阻尼矩陣，M為質(zhì)量矩陣，F(xiàn)為載荷矩陣。對(duì)模型使用無網(wǎng)格法離散化，根據(jù)上文的理論可實(shí)現(xiàn)形函數(shù)的計(jì)算，代入系統(tǒng)方程后，各個(gè)矩陣計(jì)算形式為：

(9)

(10)

(11)

(12)

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

當(dāng)在實(shí)際中考慮內(nèi)部空間座椅等問題時(shí)，通常認(rèn)為聲空間中包含座椅體積所占據(jù)的空間，聲學(xué)模型中包含了座椅體積所占據(jù)的空間，這里是通過定義座椅空間部分為不同的流體材料來進(jìn)行模擬計(jì)算的[4]。

1.3無網(wǎng)格邊界元模型

對(duì)于外部聲場，參考經(jīng)典的Kirchhoff-Helmholtz方程

(13)

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

Cp=GG·p-GH·p+F

(14)

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

(15)

(16)

(17)

(18)

(19)

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

(20)

(21)

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

1.4無網(wǎng)格FEM-BEM算法

為更清楚地說明本文方法的實(shí)現(xiàn)過程，圖4給出了本文方法的計(jì)算流程圖。

圖4　混合模型計(jì)算流程圖Fig.4 Hybrid computing flowchart

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

2數(shù)值算例及結(jié)果分析

2.1算例1

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

圖5　圓柱殼結(jié)構(gòu)示意圖Fig.5 Schematic of the cylindrical shell structure

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

圖6　圓柱形結(jié)構(gòu)內(nèi)部及截面無網(wǎng)格節(jié)點(diǎn)分布Fig.6 Nodal distribution of the meshless model

圖7　各場點(diǎn)聲傳遞函數(shù)Fig.7 Sound transfer functions at different field-points

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

圖8　p2場點(diǎn)幅值頻率響應(yīng)Fig.8 frequency response of the magnitude at p2

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

2.2算例2

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

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

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

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

圖10　聲傳遞函數(shù)對(duì)比Fig.10 Comparison of transfer functions using different methods

3實(shí)驗(yàn)驗(yàn)證

為了進(jìn)一步驗(yàn)證無網(wǎng)格有限元-邊界元模型的正確性，在半消聲室中對(duì)一個(gè)實(shí)際艙段進(jìn)行了測量實(shí)驗(yàn)。測量系統(tǒng)如圖11所示，圖12為艙段模型。圖13為傳聲器位置分布示意圖。

圖11　實(shí)驗(yàn)測量系統(tǒng)Fig.11 Experimental measuring principle

圖12　艙段模型Fig.12 Cabin model

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

以測點(diǎn)2為例。圖14為數(shù)值計(jì)算和實(shí)驗(yàn)測試的1/3倍頻帶聲壓級(jí)結(jié)果的對(duì)比圖。圖15為測點(diǎn)結(jié)果的誤差曲線。“fp2”表示測點(diǎn)2。

由圖可知，在0～200 Hz時(shí)，數(shù)值計(jì)算結(jié)果和實(shí)驗(yàn)測量結(jié)果偏差在5 dB以內(nèi)，而在200～500 Hz時(shí)，實(shí)驗(yàn)結(jié)果和數(shù)值仿真結(jié)果的誤差大于5 dB，250 Hz時(shí)，聲壓級(jí)的誤差最大。根據(jù)整體的聲壓級(jí)變化趨勢，測點(diǎn)2的數(shù)值計(jì)算結(jié)果和實(shí)驗(yàn)測量結(jié)果基本保持一致。

圖14　測點(diǎn)2聲壓級(jí)結(jié)果對(duì)比Fig.14 SPL comparison at fp2

圖15　測點(diǎn)2聲壓級(jí)誤差曲線Fig.15 SPL error at fp2

為計(jì)算各測點(diǎn)的平均聲壓級(jí)，將實(shí)驗(yàn)測得的1/3倍頻帶聲壓級(jí)結(jié)果進(jìn)行對(duì)數(shù)平均，同時(shí)將數(shù)值仿真計(jì)算得到的各測點(diǎn)的聲傳遞函數(shù)結(jié)果也作對(duì)數(shù)平均。表1為各測點(diǎn)的實(shí)驗(yàn)測量和數(shù)值計(jì)算總聲壓級(jí)結(jié)果和相對(duì)誤差。

表1　數(shù)值計(jì)算與測量結(jié)果對(duì)比

根據(jù)表1中數(shù)據(jù)可求得7個(gè)測點(diǎn)的平均相對(duì)誤差為5.26%，說明無網(wǎng)格FEM-BEM法與實(shí)際測試結(jié)果的總體偏差很小，這表明無網(wǎng)格FEM-BEM計(jì)算模型是正確的，可以用于外部聲源作用下的封閉結(jié)構(gòu)內(nèi)部的聲場計(jì)算問題。

4結(jié)論

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

參 考 文 獻(xiàn)

[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 聲學(xué)仿真計(jì)算高級(jí)應(yīng)用實(shí)例[M]. 北京：國防工業(yè)出版社,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] 王海濤,曾向陽,陳玲. 小尺度封閉空間無網(wǎng)格Galerkin聲場數(shù)值計(jì)算方法[J]. 西北工業(yè)大學(xué)學(xué)報(bào),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] 王海濤,曾向陽. 周期結(jié)構(gòu)聲散射系數(shù)的無網(wǎng)格數(shù)值計(jì)算方法[J]. 計(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] 張雄,劉巖,馬上. 無網(wǎng)格法的理論及應(yīng)用[J]. 力學(xué)進(jìn)展,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] 張雄,劉巖. 無網(wǎng)格法[M]. 北京:清華大學(xué)出版社,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

基金項(xiàng)目：國家自然科學(xué)基金(11374241)；航空科學(xué)基金(20151553021)

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

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

中圖分類號(hào):O422

文獻(xiàn)標(biāo)志碼：A

DOI:10.13465/j.cnki.jvs.2016.09.002

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