尤小健，李 威，柴應(yīng)彬

(1.武漢第二船舶設(shè)計(jì)研究所,湖北 武漢 430064； 2.華中科技大學(xué) 船舶與海洋工程學(xué)院，湖北 武漢 430074； 3.船舶和海洋水動(dòng)力湖北省重點(diǎn)實(shí)驗(yàn)室，湖北 武漢 430074；4.高新船舶與深海開(kāi)發(fā)裝備協(xié)同創(chuàng)新中心，上海 200240)

基于面光滑有限元法的三維聲傳播問(wèn)題研究

尤小健1，李 威2,3,4，柴應(yīng)彬2

(1.武漢第二船舶設(shè)計(jì)研究所,湖北 武漢 430064； 2.華中科技大學(xué) 船舶與海洋工程學(xué)院，湖北 武漢 430074； 3.船舶和海洋水動(dòng)力湖北省重點(diǎn)實(shí)驗(yàn)室，湖北 武漢 430074；4.高新船舶與深海開(kāi)發(fā)裝備協(xié)同創(chuàng)新中心，上海 200240)

為解決標(biāo)準(zhǔn)有限元方法“過(guò)剛”的特性，采用一種基于面的光滑有限元法(FS-FEM)來(lái)求解三維聲學(xué)問(wèn)題，在FS-FEM模型當(dāng)中，應(yīng)用四節(jié)點(diǎn)的四面體單元來(lái)離散三維聲場(chǎng)域，通過(guò)與四面體單元面相關(guān)的光滑域并運(yùn)用伽遼金弱形式來(lái)建立線性離散系統(tǒng)的控制方程，在面光滑域上進(jìn)行的梯度光滑操作能夠“軟化”原來(lái)“過(guò)剛”的FEM模型，因此能夠減少數(shù)值離散誤差。數(shù)值仿真結(jié)果表明：運(yùn)用相同的網(wǎng)格密度，F(xiàn)S-FEM能夠得到比FEM更加精確的計(jì)算結(jié)果。

面光滑有限元法；梯度光滑技術(shù)；聲學(xué)；數(shù)值方法

0 引 言

充液管道在機(jī)械、船舶和石油化工等許多工程領(lǐng)域有著極為廣泛的應(yīng)用背景，因此對(duì)充液管道內(nèi)聲傳播及聲輻射特性問(wèn)題的研究具有重要的學(xué)術(shù)價(jià)值和工程意義，國(guó)內(nèi)外許多學(xué)者在這個(gè)研究領(lǐng)域做了大量的工作[1-2]。

對(duì)于一些簡(jiǎn)單的結(jié)構(gòu)，其聲學(xué)問(wèn)題的解析解是可以得到的。但是對(duì)于現(xiàn)實(shí)當(dāng)中大多數(shù)的復(fù)雜工程結(jié)構(gòu)，往往沒(méi)有解析解。此時(shí)通常用各種各樣的數(shù)值方法來(lái)得到其近似解，其中有限元方法和邊界元方法是求解聲學(xué)問(wèn)題當(dāng)中發(fā)展的最成熟，應(yīng)用最廣泛的數(shù)值方法[3-5]。在運(yùn)用有限元方法進(jìn)行聲場(chǎng)數(shù)值模擬時(shí)，由于有限元模型總體剛度與系統(tǒng)的真實(shí)剛度相比總是表現(xiàn)的“偏剛”，這就不可避免地就會(huì)出現(xiàn)數(shù)值計(jì)算誤差[6]，而且這種數(shù)值計(jì)算誤差會(huì)隨著計(jì)算頻率的提高而不斷變大。為了盡量減少這種數(shù)值誤差，提高聲學(xué)數(shù)值仿真的計(jì)算精度，通常的做法就是不斷加大有限元網(wǎng)格的密度，從而使有限元模型的剛度更加接近模型的真實(shí)剛度。但是隨著網(wǎng)格密度的增加，計(jì)算成本也大大增加，尤其對(duì)于復(fù)雜的三維聲學(xué)問(wèn)題，這個(gè)問(wèn)題就顯得更加突出。

為了解決上述有限元模型表現(xiàn)的“過(guò)剛”的問(wèn)題，國(guó)內(nèi)外許多學(xué)者在這方面做了大量的工作。近年來(lái)，Liu等[7-10]提出的光滑有限元方法(SFEM)是將由Chen等[11-12]提出的無(wú)網(wǎng)格分區(qū)應(yīng)力光滑技術(shù)與傳統(tǒng)標(biāo)準(zhǔn)的有限元結(jié)合在一起的一種新型數(shù)值方法。該方法通過(guò)對(duì)分區(qū)光滑域的光滑操作，將原來(lái)“過(guò)剛”的有限元計(jì)算模型進(jìn)行適當(dāng)“軟化”，使其能夠更加接近計(jì)算模型的真實(shí)剛度，因而能夠得到更加精確的數(shù)值計(jì)算結(jié)果。本文應(yīng)用一種基于面的光滑有限元方法(FS-FEM)計(jì)算三維充液管道內(nèi)部的聲傳播問(wèn)題，計(jì)算結(jié)果表明，與標(biāo)準(zhǔn)有限元法相比，F(xiàn)S-FEM具有更高的計(jì)算精度和更好的收斂性，具有良好的工程應(yīng)用前景。

1 聲學(xué)光滑有限元計(jì)算模型

在理想流體介質(zhì)中，忽略流體粘性的影響，簡(jiǎn)諧小振幅振動(dòng)聲波滿足以下的波動(dòng)方程：

(1)

式中：▽2為L(zhǎng)aplace算子；p為聲場(chǎng)中任意一點(diǎn)的聲壓；c為聲速；t為時(shí)間。

在聲場(chǎng)當(dāng)中任意一點(diǎn)的質(zhì)點(diǎn)振動(dòng)速度和該點(diǎn)的聲壓有以下關(guān)系：

▽p+jρωv=0，

(2)

式中：ρ為流體介質(zhì)的密度；ω為圓頻率；v為質(zhì)點(diǎn)的振動(dòng)速度。

在聲場(chǎng)計(jì)算中，考慮3種聲學(xué)邊界條件：Dirichlet邊界條件，Neumann邊界條件和Robin邊界條件，這3類控制方程的表述如下[9,13]：

p=pD，

(3)

▽p·n=-jρωvn，

(4)

▽p·n=-jρωAnp。

(5)

式中vn和An分別為質(zhì)點(diǎn)振動(dòng)的法向速度和邊界上的導(dǎo)納系數(shù)。

采用加權(quán)值法并應(yīng)用上述邊界條件，可得式(1)的Galerkin弱積分形式[9]:

(6)

式中Γ1和Γ2分別為式(4)給定的速度邊界條件和式(5)給定的聲阻抗邊界條件。聲壓變量由式(7)得到：

(7)

式中：Ni為形函數(shù)；pi為節(jié)點(diǎn)聲壓。

運(yùn)用廣義梯度光滑技術(shù)對(duì)聲場(chǎng)當(dāng)中的光滑域進(jìn)行分片光滑操作，可得如下光滑的Galerkin弱積分形式：

(8)

由式(8)可得到如下的聲學(xué)有限元計(jì)算模型[9]：

(9)

式中：

(10)

(11)

為光滑剛度矩陣；

(12)

(13)

上述光滑剛度矩陣不同于標(biāo)準(zhǔn)有限元方法中的剛度矩陣，它通過(guò)在每個(gè)光滑域上而不是在每個(gè)單元內(nèi)部進(jìn)行數(shù)值積分得到。這種操作是對(duì)原來(lái)“過(guò)剛”的有限元模型的一種修正。

2 基于面的光滑域構(gòu)造

對(duì)任意一個(gè)三維問(wèn)題，每個(gè)光滑域的構(gòu)造是基于每個(gè)單元的面，構(gòu)造過(guò)程是：首先將三維問(wèn)題域離散為4節(jié)點(diǎn)的四面體單元，然后將相鄰四面體的中點(diǎn)及每個(gè)面上的3個(gè)頂點(diǎn)順序相連，這就形成了基于面的光滑域。因此，整個(gè)問(wèn)題域當(dāng)中有多少個(gè)面，就會(huì)形成多少個(gè)基于面的光滑域[9,14]，具體如圖1所示。

圖1 基于四面體網(wǎng)格的光滑域Fig.1 The smoothing domain based on tetrahedral mesh

在運(yùn)用廣義梯度光滑技術(shù)對(duì)原有的梯度場(chǎng)進(jìn)行光滑操作時(shí)實(shí)際上是對(duì)每一個(gè)光滑域上的速度梯度場(chǎng)進(jìn)行光滑，光滑后的速度由式(14)得到：

(14)

式中Vk為基于面k的光滑域的體積。將式(7)代入式(14)可得：

(15)

式中n為與每個(gè)光滑域相關(guān)的節(jié)點(diǎn)數(shù)。

對(duì)于三維問(wèn)題

(16)

(17)

式中：Vk為每個(gè)光滑域的體積；Ns為每個(gè)光滑域邊界被分成的份數(shù)；Ng為每份積分域當(dāng)中高斯點(diǎn)的個(gè)數(shù)；w為與之相應(yīng)的權(quán)系數(shù)。

進(jìn)一步可得每個(gè)光滑域Ωk的單元?jiǎng)偠染仃嚕?/p>

(18)

和標(biāo)準(zhǔn)有限元類似，將每個(gè)基于面的光滑域上的單元?jiǎng)偠染仃嚳傮w組裝，就可得到系統(tǒng)的總體剛度矩陣。

3 數(shù)值計(jì)算結(jié)果

三維充液管道計(jì)算模型如圖3所示，管道的長(zhǎng)度L=1m，半徑R=0.05m，管道內(nèi)部區(qū)域Ω充滿密度ρ=1 000kg/m3, 聲速c=1 500m/s的流體介質(zhì)。管道的邊界條件為：在左端面施加一個(gè)vn=10sinωtm/s的速度邊界條件，管道的右端和管壁為剛性邊界條件。

圖2 三維充液管道模型Fig.2 The three dimensional tube model filled with liquid

3.1 模態(tài)分析

為了驗(yàn)證本文所提出計(jì)算方法的有效性與可行性，首先對(duì)圖3所示的計(jì)算模型進(jìn)行模態(tài)分析，計(jì)算過(guò)程中不考慮阻尼的影響。計(jì)算網(wǎng)格如圖3所示，由于對(duì)稱性，圖3只給出了1/4的網(wǎng)格模型，該網(wǎng)格模型包含894個(gè)節(jié)點(diǎn)和3 141個(gè)單元，單元的平均尺寸為h=0.025m。根據(jù)聲學(xué)網(wǎng)格單元尺寸所要滿足的“拇指法則”，即每個(gè)計(jì)算頻率所對(duì)應(yīng)的波長(zhǎng)至少應(yīng)包含6個(gè)單元[9,17]，該網(wǎng)格模型下可分析的上限頻率為10 000Hz。表1給出了FEM和FS-FEM在相同網(wǎng)格模型下前10階模態(tài)的計(jì)算結(jié)果，以及解析解結(jié)果。

圖3 三維充液管道1/4網(wǎng)格模型Fig.3 A quarter mesh of the three dimensional tube model filled with liquid

從表1可看出：在低階模態(tài)時(shí)，F(xiàn)EM和FS-FEM都可得到較為理想的計(jì)算結(jié)果，但與FEM相比，F(xiàn)S-FEM的計(jì)算精度更高。隨著模態(tài)階數(shù)的增高，F(xiàn)EM和FS-FEM的計(jì)算誤差都逐漸增大，但是，F(xiàn)S-FEM的計(jì)算結(jié)果仍然要優(yōu)于FEM。這說(shuō)明通過(guò)前面所述的“梯度光滑技術(shù)”得到的FS-FEM計(jì)算模型有效地“軟化”了原來(lái)“過(guò)剛”的FEM計(jì)算模型，因而能夠得到更加精確的計(jì)算結(jié)果。

表1 三維充液管道模態(tài)計(jì)算結(jié)果Tab.1 The results of the three dimensional tube eigenfrequencies

3.2 聲傳播問(wèn)題

分別運(yùn)用FEM和FS-FEM在前述的邊界條件下對(duì)圖2所示的管道模型的聲傳播問(wèn)題進(jìn)行計(jì)算和分析，計(jì)算網(wǎng)格仍為圖3所示的網(wǎng)格模型。為了考察聲傳播頻率對(duì)計(jì)算結(jié)果的影響，這里計(jì)算了f=2 000Hz，f=4 000Hz和f=7 000Hz等3種不同頻率下管道內(nèi)的聲壓幅值分布情況，結(jié)果如圖4～圖6所示，另外該聲場(chǎng)問(wèn)題的解析解為[9]：

(19)

(20)

為了便于比較，圖中也給出了解析解的計(jì)算結(jié)果。

從圖4～圖6可看出，在低頻區(qū)域(1 000Hz)，F(xiàn)EM和FS-FEM的計(jì)算結(jié)果都與解析解吻合得很好；在高頻區(qū)域，F(xiàn)EM和FS-FEM結(jié)果的誤差都增加，但FS-FEM的結(jié)果要優(yōu)于FEM，這就進(jìn)一步驗(yàn)證了FS-FEM的精確性和有效性。

圖4 聲傳播頻率f=2 000 Hz聲壓幅值分布圖Fig.4 Acoustic pressure distribution along the tube with frequency f=2 000 Hz

圖5 聲傳播頻率f=4 000 Hz聲壓幅值分布圖Fig.5 Acoustic pressure distribution along the tube with frequency f=4 000 Hz

圖6 聲傳播頻率f=7 000 Hz聲壓幅值分布圖Fig.6 Acoustic pressure distribution along the tube with frequency f=7 000 Hz

3.3 收斂性分析

眾所周知，運(yùn)用有限元等數(shù)值方法求解聲學(xué)問(wèn)題時(shí)不可避免地會(huì)出現(xiàn)數(shù)值誤差，數(shù)值誤差越小，數(shù)值方法的收斂性就越好。為了驗(yàn)證FS-FEM的收斂性，本文給出了聲學(xué)數(shù)值計(jì)算結(jié)果的全局部誤差因子[15-16]：

(16)

3.3.1 單元尺寸對(duì)全局誤差的影響

網(wǎng)格劃分的質(zhì)量對(duì)數(shù)值方法的計(jì)算結(jié)果有著重要的影響，一般情況下，網(wǎng)格單元的平均尺寸越小，網(wǎng)格越精細(xì)，計(jì)算結(jié)果的誤差就越小。這里給出了4種不同的網(wǎng)格，單元平均尺寸分別為h=0.025m，h=0.02m，h=0.015m和h=0.012m。計(jì)算頻率f=4 000Hz，該計(jì)算頻率滿足上述4種聲學(xué)網(wǎng)格所決定的上限頻率，圖7給出了FEM和FS-FEM在不同單元尺寸網(wǎng)格下的全局誤差。從圖7可看出：隨著單元平均尺寸的減小，網(wǎng)格密度的增大，F(xiàn)EM和FS-FEM的全局誤差因子都迅速減小，可見(jiàn)這2種方法都是收斂的，但在同一單元尺寸下，F(xiàn)S-FEM的全局誤差要小于FEM，這說(shuō)明FS-FEM的收斂性要優(yōu)于FEM。

圖7 全局誤差與單元平均尺寸之間的關(guān)系Fig.7 The global error against the mesh size

3.3.2 計(jì)算頻率對(duì)全局誤差的影響

在用數(shù)值方法求解聲學(xué)問(wèn)題時(shí)，計(jì)算頻率的影響也是至關(guān)重要的，通常在低頻區(qū)域，很多數(shù)值方法都可得到較為理想的數(shù)值結(jié)果，隨著頻率的提高，數(shù)值方法的誤差就會(huì)變大，有時(shí)甚至不能得到可靠的結(jié)果。本小節(jié)主要討論計(jì)算頻率對(duì)全局誤差及收斂性的影響。網(wǎng)格計(jì)算模型如圖4所示，單元平均尺寸為0.025m，為滿足聲學(xué)網(wǎng)格中每個(gè)波長(zhǎng)至少要包含6個(gè)單元的要求，這里的分析頻率為0～9 000Hz，計(jì)算結(jié)果如圖8所示。

圖8 全局誤差與聲傳播頻率之間的關(guān)系Fig.8 The global error against the frequency

從圖8可看出，在低頻區(qū)域，F(xiàn)EM和FS-FEM的全局誤差都較小，隨著頻率的提高，F(xiàn)EM和FS-FEM的全局誤差都逐漸增大，但在整個(gè)計(jì)算頻率范圍內(nèi)FS-FEM的全局誤差要小于FEM，這進(jìn)一步說(shuō)明了FS-FEM的計(jì)算精度更高，收斂性更好。另外從圖中可看出，全局誤差在系統(tǒng)的各階固有頻率附近會(huì)急劇增大，說(shuō)明在系統(tǒng)的固有頻率附近，F(xiàn)EM和FS-FEM的計(jì)算誤差都較大，這與先前其他學(xué)者所得的結(jié)果相吻合[15-16]。出現(xiàn)這種現(xiàn)象的原因是：在求解式(9)時(shí)，當(dāng)聲波的傳播頻率和系統(tǒng)的固有頻率相近時(shí)，導(dǎo)致式(9)的系數(shù)矩陣特性接近于一個(gè)奇異陣，這樣在計(jì)算過(guò)程當(dāng)中任何微小的數(shù)值誤差都會(huì)導(dǎo)致求解聲壓時(shí)出現(xiàn)較大的誤差。

4 結(jié) 語(yǔ)

本文采用基于面的光滑有限元法(FS-FEM)對(duì)三維充液管道內(nèi)的聲傳播問(wèn)題進(jìn)行研究，通過(guò)對(duì)計(jì)算結(jié)果的分析和比較可以得到如下的結(jié)論：

1)在三維充液管道的模態(tài)預(yù)報(bào)中，無(wú)論是低階模態(tài)還是高階模態(tài)，F(xiàn)S-FEM都能得到比FEM更加精確的結(jié)果。

2)在聲場(chǎng)計(jì)算中，當(dāng)聲傳播頻率較低時(shí)，F(xiàn)EM和FS-FEM都可以得到較為理想的結(jié)果，但隨著聲傳播頻率的提高，F(xiàn)S-FEM的計(jì)算精度要明顯高于FEM。

3)隨著網(wǎng)格單元平均尺寸的減小，F(xiàn)EM和FS-FEM的全局誤差都迅速減小，但FS-FEM的全局誤差始終要小于FEM,表明FS-FEM在收斂性方面要優(yōu)于FEM。

[1] 劉敬喜，李天勻，劉土光，等.彈性介質(zhì)中充液管道的波衰減特性[J].華中科技大學(xué)學(xué)報(bào)(自然科學(xué)版),2003,31(10):90-92.

[2]PER-ANDERSH,G?RANS.Dynamicfiniteelementanalysisoffluid-filledpipes[J].ComputerMethodsinAppliedMechanicsandEngineering,2001,190:3111-3120.

[3]CHENGL,ROBERTD,WHITEKG.Three-dimensionalviscousfiniteelementformulationforacousticfluid-structureinteraction[J].ComputerMethodsinAppliedMechanicsandEngineering,2008,197:60-72.

[4]CHENGH,CHENJ,ZHANGYB,etal.Amulti-domainboundaryelementformulationforacousticfrequencysensitivityanalysis[J].EngineeringAnalysiswithBoundaryElements,2009,33:815-821.

[5]WUHJ,LIUYJ,JIANGWK.Alow-frequencyfastmultipoleboundaryelementmethodbasedonanalyticalintegrationofthehypersingularintegralfor3Dacousticproblems[J].EngineeringAnalysiswithBoundaryElements,2013,37:309-318.

[6]BOUILLARDPH,SULEAUS.Element-freegarlekinsolutionsforhelmholtzproblems:formulationandnumericalassessmentofthepollutioneffect[J].ComputerMethodsinAppliedMechanicsandEngineering,1998,162:117-135.

[7]LIUGR,NGUYENTT,LAMKY.Anedge-basedsmoothedfiniteelementmethod(ES-FEM)forstaticfree,andforcedvibrationanalysis[J].JournalofSoundandVibration,2009,320:1100-1130.

[8]CUIXY,LIUGR,LIGY,etal.Analysisofplatesandshellsusingedge-basedsmoothedfiniteelementmethod[J].ComputationalMechanics,2010,45:141-156.

[9]HEZC,LIUGR,ZHONGZH,etal.Anedge-basedsmoothedfiniteelementmethod(ES-FEM)foranalyzingthree-dimensionalacousticproblems[J].ComputerMethodsinAppliedMechanicsandEngineering，2009,199:20-33.

[10]HEZC,LIUGR,ZHONGZH,etal.AcoupledES-FEM/BEMmethodforfluid-structureinteractionproblems[J].EngineeringAnalysiswithBoundaryElements,2011,35:140-147.

[11]CHENJS,WUCT,YOONS,etal.Astabilizedconformingnodalintegrationorgalerkinmeshfreemethods[J].InternationalJournalforNumericalMethodsinEngineering,2001,50:435-466.

[12]YOOJW,MORANB,CHENJS.Stabilizedconformingnodalintegrationinthenatural-elementmethod[J].InternationalJournalforNumericalMethodsinEngineering,2004,60:861-890.

[13] 何智成，李光耀，成艾國(guó)，等.光滑有限元的聲學(xué)研究：時(shí)域和頻域分析[J].振動(dòng)與沖擊,2012,31(16):122-127.

[14]HEZC,LIUGR,ZHONGZH,etal.Acouplededge-/face-basedsmoothedfiniteelementmethodforstructural-acousticproblems[J].AppliedAcoustics,2010,71:1955-1964.

[15] 夏百戰(zhàn)，于德介，姚凌云.光滑有限元的聲學(xué)研究：二維多流體域耦合聲場(chǎng)的光滑有限元解法[J].聲學(xué)學(xué)報(bào),2012,37(6):601-609.

[16]HEZC,LIUGR,ZHONGZH,etal.Dispersionfreeanalysisofacousticproblemsusingthealphafiniteelementmethod[J].ComputationalMechanics，2010,46:867-881.

[17]HEZC,LIUGR,ZHONGZH,etal.Coupledanalysisof3Dstructural-acousticproblemsusingtheedge-basedsmoothedfiniteelementmethod/finiteelementmethod[J].FiniteElementsinAnalysisandDesign，2010，46：1114-1121.

Analysis of sound of propagation in three-dimensional fluid-filled pipes based on face-based on smoothed finite element method (FS-FEM)

YOU Xiao-jian1，LI Wei2,3,4，CHAI Ying-bin2

(1.Wuhan Second Ship Design and Research Institute,Wuhan 430064,China； 2.School of Naval Architecture and Ocean Engineering, Huazhong University of Science and Technology,Wuhan 430074,China； 3.Hubei Key Laboratory of Naval Architecture & Ocean Engineering Hydrodynamics (HUST) ,Wuhan 430074,China； 4.Collaborative Innovation Center for Advanced Ship and Deep-Sea Exploration (CISSE) ,Shanghai 200240,China)

To solve the “overly-stiff” property of the standard FEM, an face-based smoothed finite element method (FS-FEM) was adopted for analyzing three-dimensional acoustic problems.In the FS-FEM model, four-node tetrahedron elements are used to discretize the three-dimensional acoustic domain, the discretized linear system equations are established using the smoothed Galergin weak form with smoothing domains associated the surface of the tetrahedrons.The gradient smoothing operations over the faced-based smoothing domains provides a proper softening effect to the “overly-stiff” FEM model and hence can reduce the numerical dispersion error.The numerical simulation results show that the FS-FEM model can achieve more accurate results than the FEM model with the same mesh owing to the softened stiffness to the acoustic model.

face-based smoothed finite element method(FS-FEM);gradient smoothing technique;acoustics;numerical method

2015-01-04；

2015-02-02

尤小健(1975-)，男，高級(jí)工程師，研究方向?yàn)榕灤傮w設(shè)計(jì)。

O42

A

1672-7649(2015)12-0104-06

10.3404/j.issn.1672-7649.2015.12.021