基于Abaqus的水下結(jié)構(gòu)聲輻射仿真方法

吳健++李澤成++熊晨熙

摘要： 提出一種在Abaqus中直接完成輻射聲場(chǎng)仿真的方法，該方法可大大減少流場(chǎng)網(wǎng)格層數(shù).利用無(wú)限元的計(jì)算結(jié)果求解任意遠(yuǎn)處球形聲場(chǎng)的聲壓；結(jié)合彈性球殼的解析解算例驗(yàn)證該方法的正確性，并將該方法拓展到細(xì)長(zhǎng)型結(jié)構(gòu)中.結(jié)果表明：基于Abaqus無(wú)限元的水下結(jié)構(gòu)聲輻射的計(jì)算方法簡(jiǎn)單、可行；采用有限元與無(wú)限元相結(jié)合的方法可以大幅度縮小流場(chǎng)區(qū)域.

關(guān)鍵詞： 水下結(jié)構(gòu)； 聲學(xué)無(wú)限元； 聲介質(zhì)； 聲輻射； 球殼

中圖分類(lèi)號(hào)： TB532文獻(xiàn)標(biāo)志碼： B

0引言

隨著現(xiàn)代潛艇隱身技術(shù)的發(fā)展，帶肋圓柱殼體作為潛艇的典型結(jié)構(gòu)形式，其在無(wú)限介質(zhì)中的流固耦合振動(dòng)和聲輻射問(wèn)題日益成為研究的熱點(diǎn).采用仿真分析的方式研究水下結(jié)構(gòu)聲輻射問(wèn)題可節(jié)省大量的試驗(yàn)費(fèi)用，是當(dāng)前的主要研究方法.

一般來(lái)說(shuō)，水下結(jié)構(gòu)的振動(dòng)聲學(xué)分析方法主要有解析法和數(shù)值法2類(lèi).解析法只適用于計(jì)算少數(shù)具有正交曲面的彈性結(jié)構(gòu)（如球面和圓柱面等結(jié)構(gòu)）的輻射場(chǎng).采用解析法求解復(fù)雜結(jié)構(gòu)的聲學(xué)問(wèn)題很困難，因此對(duì)于一般形狀的物體應(yīng)采用數(shù)值法進(jìn)行求解.目前，聲輻射的數(shù)值法很多基于模態(tài)疊加法，由于高頻模態(tài)較為密集且結(jié)構(gòu)模態(tài)阻尼難以準(zhǔn)確測(cè)量，因此直接法求解結(jié)構(gòu)聲輻射有很多優(yōu)勢(shì).商德江[1]提出首先在有限元中計(jì)算模型在激勵(lì)力作用下的流固耦合動(dòng)力響應(yīng)，再利用結(jié)構(gòu)表面的位移作為邊界元的邊界條件，計(jì)算聲學(xué)物理量的計(jì)算方法.姚雄亮等[2]將ANSYS與Sysnoise相結(jié)合對(duì)雙層殼體水下振動(dòng)和聲輻射進(jìn)行仿真分析，認(rèn)為利用2種軟件相結(jié)合的方法對(duì)水下復(fù)雜結(jié)構(gòu)的振動(dòng)和聲輻射進(jìn)行數(shù)值計(jì)算完全可行.繆旭弘等[3]通過(guò)Abaqus和聲固耦合算法對(duì)雙層圓柱殼的水下噪聲輻射特性進(jìn)行計(jì)算和分析，計(jì)算值與實(shí)驗(yàn)數(shù)據(jù)吻合較好.以上方法均采用聲學(xué)無(wú)限元或邊界阻抗技術(shù)的處理方式實(shí)現(xiàn)無(wú)限大流場(chǎng)，一般需要建立足夠大的水體模型來(lái)滿足聲學(xué)邊界條件，水體單元需要涵蓋到聲壓分析位置，或者使用Sysnoise軟件求解任意遠(yuǎn)場(chǎng)點(diǎn)的聲壓.

Abaqus具有較強(qiáng)的建模與計(jì)算分析能力，但在聲學(xué)后處理方面稍顯薄弱.本文提出一種在Abaqus中直接完成輻射聲場(chǎng)仿真的方法，利用彈性球殼的解析解算例驗(yàn)證不同建模方式時(shí)聲輻射結(jié)果的正確性.結(jié)果表明：Abaqus中的無(wú)限元技術(shù)可以采用較小流場(chǎng)進(jìn)行水下結(jié)構(gòu)振動(dòng)分析，求解得到的任意遠(yuǎn)場(chǎng)點(diǎn)的聲輻射結(jié)果精度滿足分析要求.

1理論背景

1.1有限元方程

有限元對(duì)模型進(jìn)行分析是通過(guò)對(duì)模型進(jìn)行離散后求解波動(dòng)方程和運(yùn)動(dòng)方程實(shí)現(xiàn)的.在流固耦合面上，結(jié)構(gòu)振動(dòng)會(huì)產(chǎn)生流體負(fù)載，而聲壓同時(shí)對(duì)結(jié)構(gòu)產(chǎn)生附加力，所以必須同時(shí)計(jì)算結(jié)構(gòu)動(dòng)力方程和流體域的波動(dòng)方程，用統(tǒng)一的矩陣形式表示為Ms

ρRf 0

Mf u¨

p¨+ Cs

0 0

Cf u·

p·+ Ks

0 -Rf

Kfu

p=Fs

0（1）式中：Ms為結(jié)構(gòu)質(zhì)量矩陣；Ks為結(jié)構(gòu)剛度矩陣；Cs為結(jié)構(gòu)阻尼矩陣；Mf為流體質(zhì)量矩陣；Kf為流體剛度矩陣；Cf為聲阻尼矩陣，Rf為流固面上的耦合矩陣；ρ為流體介質(zhì)密度；u和p為節(jié)點(diǎn)的位移和聲壓向量；Fs為結(jié)構(gòu)的載荷向量.用式（1）可以同時(shí)得到流固面上的位移和聲壓值.

1.2無(wú)限流場(chǎng)模擬

聲波在有限流場(chǎng)中傳播時(shí)，由于邊界阻抗的存在，在流場(chǎng)的邊界必然存在聲波的反射.采用聲學(xué)無(wú)限元技術(shù)后，由聲學(xué)介質(zhì)組成的流場(chǎng)能夠滿足在流場(chǎng)無(wú)窮遠(yuǎn)邊界上的Sommerfield輻射條件lim rr→∞（pr+jkp）=0（2）2有限元建模

2.1網(wǎng)格劃分

有限元分析中的網(wǎng)格劃分好壞直接關(guān)系到模型計(jì)算的準(zhǔn)確性.根據(jù)波動(dòng)理論，單元長(zhǎng)度的細(xì)化標(biāo)準(zhǔn)為在一個(gè)彎曲波波長(zhǎng)內(nèi)要有4個(gè)以上單元.板的彎曲波波長(zhǎng)近似計(jì)算公式為r=δ212（1-μ2）（3）

λB=2πrCB/f（4）式中：r為板的慣性半徑；δ為板厚；μ為泊松比；λB為彎曲波波長(zhǎng)；CB為彎曲波波速；f為計(jì)算頻率上限.為適應(yīng)計(jì)算數(shù)據(jù)的分布特點(diǎn)，應(yīng)在結(jié)構(gòu)不同部位采用疏密變化的網(wǎng)格.

流體網(wǎng)格應(yīng)保證每個(gè)聲波波長(zhǎng)內(nèi)網(wǎng)格數(shù)量大于6個(gè)，同時(shí)流體與結(jié)構(gòu)耦合面網(wǎng)格尺寸應(yīng)與結(jié)構(gòu)保持一致，水體單元采用聲學(xué)單元模擬.

2.2聲學(xué)無(wú)限元技術(shù)建模

聲學(xué)無(wú)限元技術(shù)要求在結(jié)構(gòu)中設(shè)置參考點(diǎn)作為聲場(chǎng)分析的原點(diǎn)，一般應(yīng)將無(wú)限元設(shè)置成外凸形狀.在Abaqus中可以利用“skin”在水體外邊界上定義無(wú)限元，并賦予“acoustic infinite”屬性.在網(wǎng)格劃分時(shí)采用相應(yīng)的聲學(xué)無(wú)限元單元，即完成無(wú)限元邊界定義.結(jié)構(gòu)表面與水體的接觸面用關(guān)鍵字“TIE”實(shí)現(xiàn)，使用該關(guān)鍵字可以使水體與結(jié)構(gòu)表面始終保持接觸狀態(tài)而不分離.[4]

2.3場(chǎng)點(diǎn)聲壓計(jì)算

通過(guò)聲場(chǎng)分析插件“acousticVisualization”可將無(wú)限元的聲壓計(jì)算結(jié)果拓展到任意遠(yuǎn)球殼聲場(chǎng)上[5]，目前僅支持基于直接法的穩(wěn)態(tài)聲振耦合計(jì)算結(jié)果文件，且結(jié)果文件中必須包含聲學(xué)無(wú)限單元的法向矢量INFN，基于參考點(diǎn)的無(wú)限單元半徑INFR和無(wú)限單元的壓力系數(shù)PINF的輸出結(jié)果.

3仿真方法驗(yàn)證

采用受法向簡(jiǎn)諧點(diǎn)集中力作用的彈性球殼在無(wú)限大流體中的結(jié)構(gòu)聲輻射作為考核算例，具體流體結(jié)構(gòu)參數(shù)如下：球殼半徑為0.5 m；球殼壁厚為0.001 m；球殼體密度為7 800 kg/m3；彈性模量為210 GPa；泊松比為0.3；損耗因子為0.02；流體密度為1 025 kg/m3；流體聲速為1 500 m/s.

無(wú)限大流體中球殼的聲輻射解析式為

p（r，θ，t）=-ρω2∞n=0Unrh（2）n（kr）kdh（2）n（x）dxx=kR·endprint

Pn（cos θ）eiωt（5）

式中：Unr為球殼法向振動(dòng)位移的n號(hào)簡(jiǎn)正模的振幅；r為球殼半徑；Pn（）為n階勒讓德多項(xiàng)式；h（2）n（）為第二類(lèi)球漢克爾函數(shù).

應(yīng)用Abaqus程序計(jì)算流場(chǎng)輻射聲壓，兩者相比較，考核計(jì)算方法和計(jì)算程序的正確性.

計(jì)算場(chǎng)點(diǎn)處的聲源級(jí)（換算成距球心1 m處的聲壓級(jí)）Lp=20lgp（r，θ）1×10-6+20lg r（6）式中：（r，θ）為場(chǎng)點(diǎn)在以球心為原點(diǎn)的球坐標(biāo)系中的坐標(biāo)，見(jiàn)圖1.

圖 1球殼和坐標(biāo)系

Fig.1Spherical shell and coordinate system

對(duì)該模型進(jìn)行有限元建模，采用結(jié)構(gòu)化有限元網(wǎng)格劃分方法，結(jié)構(gòu)網(wǎng)格和水體網(wǎng)格分別為四邊形和六面體單元.根據(jù)式（2）和（3）求得本模型在400 Hz時(shí)的彎曲波波長(zhǎng)為0.16 m.取一個(gè)波長(zhǎng)內(nèi)8個(gè)節(jié)點(diǎn)，即結(jié)構(gòu)網(wǎng)格尺寸為0.02 m，共9 126個(gè)結(jié)構(gòu)單元.水體厚度為0.02 m，網(wǎng)格尺寸為0.02 m，水體為1層單元，共28 087個(gè)聲學(xué)單元和9 360個(gè)聲學(xué)無(wú)限元單元，見(jiàn)圖2.數(shù)值結(jié)果與解析結(jié)果的聲源級(jí)比較見(jiàn)圖3.各子圖中2條曲線均幾乎完全重合（微小差異是由于網(wǎng)格劃分等數(shù)值離散引起的），說(shuō)明本文所述的有限元加無(wú)限元的計(jì)算方法和計(jì)算程序是正確的.

圖 2結(jié)構(gòu)化有限元網(wǎng)格

Fig.2Structured finite element mesh

a）r=100，θ=π/2b）r=100，θ=0c）r=1，θ=π/2d）r=1，θ=0圖 3聲源級(jí)比較

Fig.3Acoustic source level comparisons

在實(shí)際問(wèn)題仿真過(guò)程中，結(jié)構(gòu)和水體模型往往難以劃分結(jié)構(gòu)化網(wǎng)格，且在聲學(xué)無(wú)限單元中，無(wú)限方向上（參考點(diǎn)與無(wú)限元節(jié)點(diǎn)射線方向）的聲場(chǎng)變化是通過(guò)一組包含10個(gè)9階多項(xiàng)式的函數(shù)給出的[5]，在水體模型網(wǎng)格縮小后無(wú)限單元是影響計(jì)算速度的主要因素.對(duì)2種不同流場(chǎng)網(wǎng)格劃分方式進(jìn)行結(jié)果對(duì)比：網(wǎng)格1直接采用四面體網(wǎng)格，由細(xì)到粗過(guò)渡，內(nèi)層網(wǎng)格尺寸為0.02 m，外層網(wǎng)格尺寸為0.04 m，共5 320個(gè)聲學(xué)無(wú)限元單元；網(wǎng)格2將流體網(wǎng)格分為2部分，與結(jié)構(gòu)接觸部分采用四面體網(wǎng)格，在流體邊界外側(cè)采用結(jié)構(gòu)化網(wǎng)格，且網(wǎng)格由內(nèi)至外逐漸變粗，內(nèi)層網(wǎng)格尺寸為0.02 m，外層網(wǎng)格尺寸為0.04 m，共2 946個(gè)聲學(xué)無(wú)限單元.這2種網(wǎng)格劃分方式均可明顯減少無(wú)限元網(wǎng)格的數(shù)量.有限元模型見(jiàn)圖4，數(shù)值仿真結(jié)果見(jiàn)圖5.由此可知：采用網(wǎng)格1和2的劃分方法，計(jì)算結(jié)果與解析解吻合仍然較好；100 m處聲場(chǎng)計(jì)算結(jié)果中網(wǎng)格2劃分方法精度要優(yōu)于網(wǎng)格1劃分方法，且這2種網(wǎng)格劃分的計(jì)算時(shí)間都遠(yuǎn)低于結(jié)構(gòu)化網(wǎng)格劃分方式.因此，在分析復(fù)雜結(jié)構(gòu)工程問(wèn)題時(shí)建議采用網(wǎng)格2的網(wǎng)格劃分方式.

a）網(wǎng)格1b）網(wǎng)格2圖 4不規(guī)則網(wǎng)格劃分

Fig.4Irregular mesh

a）r=100，θ=π/2b）r=100，θ=0c）r=1，θ=πd）r=1，θ=0圖 5非規(guī)則網(wǎng)格聲源級(jí)比較

Fig.5Comparisons of acoustic source level with irregular mesh

4細(xì)長(zhǎng)結(jié)構(gòu)聲輻射分析

水下潛器如潛艇、魚(yú)雷等一般為細(xì)長(zhǎng)結(jié)構(gòu)，采用有限元法進(jìn)行聲輻射分析時(shí)如果建立球形水體，水體部分會(huì)形成大量的網(wǎng)格，導(dǎo)致有限元法對(duì)硬件要求極高且計(jì)算頻率上限較低.采用前文仿真方法，針對(duì)球殼結(jié)構(gòu)模型建立圓柱水體求解結(jié)構(gòu)聲輻射.

水體邊界與結(jié)構(gòu)最小距離為0.04 m，水體有限元網(wǎng)格尺寸為0.02 m，與結(jié)構(gòu)網(wǎng)格尺寸一致，采用結(jié)構(gòu)化網(wǎng)格自動(dòng)劃分方式，網(wǎng)格見(jiàn)圖6.根據(jù)無(wú)限元分析結(jié)果進(jìn)行后處理，進(jìn)一步拓展到1和100 m處聲場(chǎng)聲壓見(jiàn)圖7，分析結(jié)果見(jiàn)圖8.有限元分析結(jié)果與理論結(jié)果吻合較好，說(shuō)明這種求解方法對(duì)水體形狀無(wú)特殊要求，因此計(jì)算細(xì)長(zhǎng)結(jié)構(gòu)聲輻射時(shí)僅需建立完全包裹結(jié)構(gòu)的圓柱水體單元.圖 6圓柱形水體網(wǎng)格

Fig.6Cylindrical mesh of water圖 7f=100 Hz時(shí)，r=1 m處聲場(chǎng)聲壓云圖，dB

Fig.7Sound field pressure contour at r=1 m while

f=100 Hz， dB

a）r=100，θ=π/2b）r=100，θ=0c）r=1，θ=πd）r=1，θ=0圖 8圓柱水體網(wǎng)格聲源級(jí)比較

Fig.8Comparisons of acoustic source level of cylindrical mesh of water

5結(jié)論

提出一種在Abaqus中有限元加無(wú)限元直接計(jì)算流固耦合結(jié)構(gòu)聲輻射的仿真方法，采用直接法進(jìn)行穩(wěn)態(tài)響應(yīng)計(jì)算，能夠求解任意距離處的聲場(chǎng)分布，通過(guò)算例驗(yàn)證計(jì)算方法和計(jì)算程序的正確性，并給出提高計(jì)算效率的網(wǎng)格劃分策略，得出以下結(jié)論.

1）該法在Abaqus中可完成任意距離的聲輻射計(jì)算，計(jì)算過(guò)程簡(jiǎn)化，結(jié)果精度滿足分析要求.

2）采用有限元加無(wú)限元的方法可以大幅度縮小流場(chǎng)，采用恰當(dāng)?shù)木W(wǎng)格劃分方法可以大大提高計(jì)算效率，特別有利于細(xì)長(zhǎng)型水下結(jié)構(gòu)的計(jì)算，降低硬件資源需求.參考文獻(xiàn)：

[1]商德江. 復(fù)雜彈性殼體水下結(jié)構(gòu)振動(dòng)和聲場(chǎng)特性研究[D]. 哈爾濱： 哈爾濱工程大學(xué)， 1999.

[2]姚熊亮， 楊娜娜， 陶景橋. 雙層殼體水下振動(dòng)和聲輻射的仿真分析[J]. 哈爾濱工程大學(xué)學(xué)報(bào)， 2004， 25（2）： 136140.

YAO Xiongliang， YANG Nana， TAO Jingqiao. Numerical research on vibration and sound radiation of underwater double cylindrical shell[J]. J Harbin Eng Univ， 2004， 25（2）： 136140.（下轉(zhuǎn)第65頁(yè)）第24卷 第6期2015年12月計(jì) 算 機(jī) 輔 助 工 程Computer Aided EngineeringVol.24 No.6Dec. 2015endprint