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

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

基于面光滑有限元法的三維聲傳播問題研究

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

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

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

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

0 引 言

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

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

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

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

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

(1)

式中：▽2為Laplace算子；p為聲場中任意一點的聲壓；c為聲速；t為時間。

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

▽p+jρωv=0，

(2)

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

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

p=pD，

(3)

▽p·n=-jρωvn，

(4)

▽p·n=-jρωAnp。

(5)

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

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

(6)

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

(7)

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

運用廣義梯度光滑技術(shù)對聲場當中的光滑域進行分片光滑操作，可得如下光滑的Galerkin弱積分形式：

(8)

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

(9)

式中：

(10)

(11)

為光滑剛度矩陣；

(12)

(13)

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

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

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

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

在運用廣義梯度光滑技術(shù)對原有的梯度場進行光滑操作時實際上是對每一個光滑域上的速度梯度場進行光滑，光滑后的速度由式(14)得到：

(14)

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

(15)

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

對于三維問題

(16)

(17)

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

進一步可得每個光滑域Ωk的單元剛度矩陣：

(18)

和標準有限元類似，將每個基于面的光滑域上的單元剛度矩陣總體組裝，就可得到系統(tǒng)的總體剛度矩陣。

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

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

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

3.1 模態(tài)分析

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

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

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

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

3.2 聲傳播問題

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

(19)

(20)

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

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

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

(16)

3.3.1 單元尺寸對全局誤差的影響

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

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

3.3.2 計算頻率對全局誤差的影響

在用數(shù)值方法求解聲學(xué)問題時，計算頻率的影響也是至關(guān)重要的，通常在低頻區(qū)域，很多數(shù)值方法都可得到較為理想的數(shù)值結(jié)果，隨著頻率的提高，數(shù)值方法的誤差就會變大，有時甚至不能得到可靠的結(jié)果。本小節(jié)主要討論計算頻率對全局誤差及收斂性的影響。網(wǎng)格計算模型如圖4所示，單元平均尺寸為0.025m，為滿足聲學(xué)網(wǎng)格中每個波長至少要包含6個單元的要求，這里的分析頻率為0～9 000Hz，計算結(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的全局誤差都逐漸增大，但在整個計算頻率范圍內(nèi)FS-FEM的全局誤差要小于FEM，這進一步說明了FS-FEM的計算精度更高，收斂性更好。另外從圖中可看出，全局誤差在系統(tǒng)的各階固有頻率附近會急劇增大，說明在系統(tǒng)的固有頻率附近，F(xiàn)EM和FS-FEM的計算誤差都較大，這與先前其他學(xué)者所得的結(jié)果相吻合[15-16]。出現(xiàn)這種現(xiàn)象的原因是：在求解式(9)時，當聲波的傳播頻率和系統(tǒng)的固有頻率相近時，導(dǎo)致式(9)的系數(shù)矩陣特性接近于一個奇異陣，這樣在計算過程當中任何微小的數(shù)值誤差都會導(dǎo)致求解聲壓時出現(xiàn)較大的誤差。

4 結(jié) 語

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

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

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

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

[1] 劉敬喜，李天勻，劉土光，等.彈性介質(zhì)中充液管道的波衰減特性[J].華中科技大學(xué)學(xué)報(自然科學(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] 何智成，李光耀，成艾國，等.光滑有限元的聲學(xué)研究：時域和頻域分析[J].振動與沖擊,2012,31(16):122-127.

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

[15] 夏百戰(zhàn)，于德介，姚凌云.光滑有限元的聲學(xué)研究：二維多流體域耦合聲場的光滑有限元解法[J].聲學(xué)學(xué)報,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-)，男，高級工程師，研究方向為艦船總體設(shè)計。

O42

A

1672-7649(2015)12-0104-06

10.3404/j.issn.1672-7649.2015.12.021