， ,b，，

(華中科技大學(xué) a.船舶與海洋工程學(xué)院；b.船舶和海洋水動(dòng)力湖北省重點(diǎn)實(shí)驗(yàn)室， 武漢 430074)

由于光波和無線電波在水中傳播時(shí)能量衰減非?？欤远叨疾槐阕鳛樗潞叫畜w探測(cè)和導(dǎo)航的信息傳遞工具。相比之下，聲波在水中具有良好的傳播性能，為人類進(jìn)行各種海洋活動(dòng)提供了很好的平臺(tái)。因此，研究水下物體的聲輻射和聲散射性質(zhì)具有重要的現(xiàn)實(shí)意義。目前求解聲學(xué)問題的數(shù)值方法有有限元法(FEM)、邊界元法(BEM)、邊界積分法、無網(wǎng)格技術(shù)和T矩陣法[1]等，其中FEM和BEM最為常用。FEM的缺點(diǎn)在于其數(shù)值離散誤差會(huì)隨計(jì)算頻率的增加而增加，為了縮小高頻段的計(jì)算誤差，雖然可以細(xì)化網(wǎng)格，但是極大地增加了計(jì)算時(shí)間，并且沒有從根本上降低誤差。

為了從根本上降低FEM的數(shù)值離散誤差，文獻(xiàn)[2]結(jié)合FEM和無網(wǎng)格技術(shù)中的應(yīng)力光滑技術(shù)提出了邊光滑有限元法(ES-FEM)。該方法認(rèn)為FEM數(shù)值誤差的來源是模型剛度太大，而ES-FEM能通過光滑技術(shù)建立剛度更接近真實(shí)情況的模型，因此能有效減小數(shù)值離散誤差。文獻(xiàn)[3]用ES-FEM研究了固體力學(xué)中的自由振動(dòng)和受迫振動(dòng)問題，發(fā)現(xiàn)ES-FEM能得到更準(zhǔn)確的固有頻率，因此位移和應(yīng)力結(jié)果也更準(zhǔn)確。文獻(xiàn)[4]用ES-FEM研究了聲固耦合問題，發(fā)現(xiàn)ES-FEM能消除剪切自鎖現(xiàn)象。

目前，國內(nèi)對(duì)ES-FEM研究較少且主要集中在固體力學(xué)領(lǐng)域，其在聲學(xué)領(lǐng)域的研究幾乎是空白。本文采用ES-FEM計(jì)算外域聲學(xué)問題，研究其在聲學(xué)領(lǐng)域的應(yīng)用價(jià)值，重點(diǎn)探討其計(jì)算效率和準(zhǔn)確性。

1 光滑有限元理論

1.1 有限元控制方程

有限元法不能直接計(jì)算無界域聲學(xué)問題，需引入一個(gè)人工邊界將無界域截?cái)嗖⒃谌斯み吔缟霞虞d一定的邊界條件，以消除人工邊界上的雜亂反射。那么，在聲學(xué)域Ω內(nèi)的控制方程有

(1)

狄利克雷邊界條件p=pD

(2)

(3)

(4)

(5)

式中：ρ——流體介質(zhì)密度;

vn——質(zhì)點(diǎn)法向振動(dòng)速度;

ω——角頻率;

An——阻尼系數(shù);

M——DtN算子，認(rèn)為人工邊界上的聲壓與其法向偏導(dǎo)數(shù)存在某種映射關(guān)系，這種關(guān)系可以用DtN算子加以描述。

對(duì)于二維問題，Givoli[5]推導(dǎo)并計(jì)算了當(dāng)人工邊界是半徑為R的圓時(shí)DtN算子的表達(dá)式

(θ-?)

(6)

采用Galerkin理論對(duì)式(1)進(jìn)行離散化處理，并帶入邊界條件得外聲場(chǎng)的總控制方程

(7)

其中剛度矩陣：

(8)

質(zhì)量矩陣：

(9)

阻尼矩陣：

(10)

人工邊界矩陣：

(11)

R——人工邊界的半徑；

NA和NB——人工邊界上的任意兩個(gè)節(jié)點(diǎn)的型函數(shù)。

力矩陣:

(12)

1.2 光滑剛度矩陣

光滑有限元法是對(duì)有限元法的一種改進(jìn)，二者求解同一問題的步驟基本相同，二者的不同之處在于總剛度矩陣的建立過程不同。把式(8)中形函數(shù)的梯度項(xiàng)換成梯度光滑項(xiàng)可得光滑有限元?jiǎng)偠染仃嚤磉_(dá)式為

(13)

(14)

聲學(xué)問題中光滑技術(shù)的本質(zhì)是對(duì)質(zhì)點(diǎn)振動(dòng)速度進(jìn)行光滑處理，光滑域Ωk內(nèi)的速度表示為

(15)

式中：W——光滑函數(shù)，其表達(dá)式為

(16)

把式(16)代入式(15)，并把速度在光滑域的積分變成聲壓在邊界上的線積分

(17)

聲壓的光滑梯度可以簡寫成

(18)

式中：Nk——光滑域的節(jié)點(diǎn)數(shù)。

(19)

式中：bId——下標(biāo)I表示節(jié)點(diǎn)號(hào)，d=1,2分別表示梯度沿x方向和y方向的值；

Γk——第k個(gè)光滑域的邊界；

Nb——光滑邊周圍單元的個(gè)數(shù)，當(dāng)光滑邊處于邊界上時(shí)取1，當(dāng)光滑變處于內(nèi)部時(shí)取2；

把式(18)和(19)帶入式(14)得到總光滑剛度陣的表達(dá)式：

(20)

2 數(shù)值算例

2.1 ES-FEM的計(jì)算效率

FEM與ES-FEM的求解過程是一致的，二者具有相同的質(zhì)量矩陣、阻尼矩陣、邊界矩陣和力矩陣，二者的惟一區(qū)別在于剛度矩陣。因此比較FEM和ES-FEM的計(jì)算效率，就是比較二者單元?jiǎng)偠染仃嚨挠?jì)算時(shí)間和總剛度矩陣的裝配時(shí)間。調(diào)整網(wǎng)格尺寸，分別采用三角形單元(T3)和四邊形單元(Q4)劃分計(jì)算域，各得到總節(jié)點(diǎn)數(shù)為576、1 196、1 620、2 176和2 660的5種模型，然后分別計(jì)算FEM和ES-FEM求解總剛度矩陣的時(shí)間，見表1。

從表1中可以觀察到：隨總結(jié)點(diǎn)數(shù)的增加，ES-FEM和FEM計(jì)算總剛度矩陣時(shí)間呈線性增加趨勢(shì)；對(duì)于T3單元和Q4單元，ES-FEM計(jì)算總剛度矩陣的時(shí)間都明顯少與FEM；ES-FEM計(jì)算總剛度矩陣時(shí)間受單元類型小于FEM。ES-FEM計(jì)算效率高于FEM的原因是：無論是T3單元還是Q4單元，ES-FEM將形函數(shù)梯度轉(zhuǎn)換為光滑域邊界積分，繼而用高斯公式將積分運(yùn)算變?yōu)榧臃ㄟ\(yùn)算，因此ES-FEM無需進(jìn)行數(shù)值微分運(yùn)算；根據(jù)式(19)得到的光滑梯度是常數(shù)值，因而ES-FEM也無需進(jìn)行數(shù)值積分運(yùn)算。

2.2 ES-FEM計(jì)算的準(zhǔn)確性

以無界域中二維舵的對(duì)平面波散射場(chǎng)計(jì)算為例，分析ES-FEM計(jì)算結(jié)果的準(zhǔn)確性。設(shè)平面波沿x軸正向傳播，其表達(dá)式為

p=exp(-jkx)

(21)

計(jì)算域是由舵的外輪廓和人工邊界所圍成的區(qū)域，用T3單元?jiǎng)澐志W(wǎng)格，調(diào)整網(wǎng)格尺寸得到34 285節(jié)點(diǎn)的精細(xì)網(wǎng)格和3 783節(jié)點(diǎn)的粗糙網(wǎng)格。由于舵表面曲率復(fù)雜不能推導(dǎo)解析解，把精細(xì)網(wǎng)格下的FEM計(jì)算結(jié)果作為參考解，見圖1a)、2a)。分別用FEM和ES-FEM計(jì)算波數(shù)k=3π(低頻)和波數(shù)k=15π(高頻)時(shí)的散射聲場(chǎng)，結(jié)果見圖1和圖2。

圖1 波數(shù)k=3π時(shí)的散射聲場(chǎng)

圖2 波數(shù)k=15π時(shí)的散射聲場(chǎng)

由圖1可見，當(dāng)入射波頻率較低時(shí)，粗糙網(wǎng)格下FEM和ES-FEM計(jì)算結(jié)果均與精細(xì)網(wǎng)格下FEM計(jì)算結(jié)果十分吻合；而由圖2中可以看出，當(dāng)入射波頻率較高時(shí)，ES-FEM解的云圖中干涉條紋和舵首尾處的散射場(chǎng)仍比較接近參考解，而同網(wǎng)格質(zhì)量下FEM解的云圖中干涉條紋不清晰，舵首尾處的散射場(chǎng)與參考解相差較大。

圖3是人工邊界上反向散射點(diǎn)的散射聲壓值隨入射波波數(shù)的變化曲線，圖中曲線是振蕩變化的。根據(jù)文獻(xiàn)[1]的解釋，對(duì)于剛性體的散射，這種振蕩現(xiàn)象是由于鏡反射波和繞射波發(fā)生干涉引起的。從圖中可以看到，當(dāng)入射波波數(shù)較小時(shí)，ES-FEM解和FEM解十分接近，并且二者曲線僅在波峰波谷位置與參考解有差異。當(dāng)入射波波數(shù)較大時(shí)，ES-FEM解與FEM解都和參考解有較大差異，但是ES-FEM解中波峰波谷處所對(duì)應(yīng)的波數(shù)值比FEM解更加接近參考解，就這一方面而言，ES-FEM解比FEM解更準(zhǔn)確。

圖3 反向散射聲壓隨波數(shù)變化曲線

3 結(jié)論

1)ES-FEM用光滑技術(shù)把形函數(shù)梯度的域積分變?yōu)橛蜻吔绲木€積分，減少了FEM數(shù)值積分和數(shù)值微分的計(jì)算時(shí)間，極大地提高了計(jì)算效率。

2)用ES-FEM計(jì)算無界域中二維舵對(duì)平面波的散射聲場(chǎng)，并與精細(xì)網(wǎng)格下FEM的計(jì)算結(jié)果進(jìn)行比較，發(fā)現(xiàn)ES-FEM在低頻段的計(jì)算結(jié)果十分準(zhǔn)確，而在高頻段的計(jì)算結(jié)果比同網(wǎng)格質(zhì)量下FEM的計(jì)算結(jié)果更準(zhǔn)確。

3)水下聲學(xué)是海洋工程的重要分支，目前對(duì)于大型水下航行體如潛艇的回波特性分析需借助工程軟件，ES-FEM對(duì)于優(yōu)化基于FEM理論的聲學(xué)工程軟件提供了新的思路。

[1] 李 威，趙 耀，張 濤，等.水下剛硬體聲散射場(chǎng)計(jì)算及分析[J].聲學(xué)技術(shù)，2007,26(5):844-849.

[2] LIU G R, NGUYEN T T,LAM K Y.An edge-based smoothed finite element method (ES-FEM) for static free, and forced vibration analysis [J]. Journal of Sound and Vibration,2009,320:1100-1130.

[3] DAI K Y,LIU G R.Free and forced vibration analysis using the smoothed finite element method (SFEM) [J]. Journal of Sound and Vibration,2007,301:803-820.

[4] HE Z C,LIU G R,ZHONG Z H,et al.Coupled analysis of 3D structural-acoustic problems using the edge-based smoothed finite element method/finite element method [J]. Finite Elements in Analysis and Design,2010,46:1114-1121.

[5] KELLER B,GIVOLI D.Exact non-reflecting boundary conditions [J]. Journal of Computational. Physics,1988,82:172-192.