一種基于流固耦合問題的低頻散射聲場預報方法

唐永壯,周其斗,謝志勇,劉文璽,趙 鵬

(海軍工程大學 艦船與海洋學院,武漢 430033)

水下聲學問題在軍事領域備受關注,在潛艇聲學隱身研究方面表現(xiàn)為潛艇輻射噪聲問題和聲散射問題。隨著輻射噪聲降低和聲吶朝著低頻、大功率方向發(fā)展,低頻聲散射研究愈發(fā)重要。Junger等[1]對剛性散射和彈性散射問題給出了較多理論推導和論述,認為結(jié)構(gòu)聲學的關鍵在于弄清流體和結(jié)構(gòu)的相互作用。Chung等[2]在圓柱殼水下共振散射試驗中檢測到彈性回波,較早發(fā)現(xiàn)彈性作用不可忽略。為研究低頻彈性散射特性和調(diào)控機制,需要理清聲散射中流固耦合機理并實現(xiàn)散射聲場準確計算。Morse等[3]給出了球形等規(guī)則形狀物體的散射理論解,但不適用于復雜結(jié)構(gòu)。湯渭霖等[4]在國內(nèi)較早研究聲散射,相繼提出亮點模型、修正亮點模型、板塊元等方法用于工程近似計算。這些方法在聲吶頻率高、不考慮彈性作用時比較有效,然而不能很好解決彈性散射問題。Zhen等[5-6]基于薄殼理論研究了柱殼結(jié)構(gòu)的彈性散射解析解,并用彈性波理論來描述散射問題,但不能完全反映水下復雜結(jié)構(gòu)的散射特性。

數(shù)值方法是研究復雜結(jié)構(gòu)聲散射特征的有力手段,Everstine[7]對多種聲學數(shù)值方法做了較詳細的論述,包括有限元法、邊界元法、無限元法等。有限元和無限元法都是不同程度近似表達波的遠場傳播條件。邊界元法可看作積分方程法的數(shù)值離散,具有理論完備性,而且可以自動滿足遠場傳播條件。文獻[8]和[9]分別用兩種不同的解耦思路開發(fā)了結(jié)構(gòu)有限元耦合流體邊界元算法應用于聲輻射計算。雖然散射與輻射本質(zhì)上都是流固耦合問題,但在機理解釋與程序算法上不能完全等同。文獻[10]用阻抗的概念闡釋了聲散射中的耦合作用,但只推導了無限長圓柱殼的情況。文獻[11]通過ANSYS和Sysnoise,實現(xiàn)了用有限元結(jié)合邊界元法計算填充球殼的散射問題。文獻[12]針對顆粒的聲散射問題提出了一種改進IB-LBM的直接模擬方法,結(jié)果表明顆粒彈性對散射聲場有較大影響。文獻[13]用有限元軟件COMSOL的完美匹配層(perfect matched layer,PML)方法計算球體的結(jié)構(gòu)表面聲壓與位移,然后通過Helmholtz積分得到遠場散射聲壓,這種方法計算精度或取決于PML能否準確模擬聲波遠場傳播條件。文獻[14]分別用ANSYS和COMSOL計算了同一簡化二維潛艇模型的聲散射,二者結(jié)果在整體上基本吻合,但部分結(jié)果相差不小。商業(yè)軟件往往對計算規(guī)模和軍事應用領域進行限制,并缺少與散射機理的對應,同時也未給出剛性散射與彈性散射的關系。

本文從邊界積分方程出發(fā)推導流體與結(jié)構(gòu)耦合方程并實現(xiàn)解耦,對散射機理做出了物理解釋,并將彈性散射分解為剛性散射項與二次輻射項。在理論推導的基礎上利用Fortran和DMAP語言編寫數(shù)值計算程序,結(jié)合Nastran的有限元功能實現(xiàn)了流固耦合作用下的散射聲場計算。最后對加肋圓柱殼的彈性散射影響因素和散射特征進行了初步研究。

1 公式推導

1.1 聲散射問題的數(shù)學表述

聲場聲壓滿足波動方程、聲壓有界條件、遠場輻射條件和邊界條件,聲散射的定解問題可表示為

(1a)

(1b)

(1c)

(1d)

式中:p為聲場中總聲壓;Γ為結(jié)構(gòu)與其外部流體相接觸的表面(下文統(tǒng)稱為結(jié)構(gòu)濕表面);r為場點到源點的距離;k為波數(shù);ω為角頻率;C0為聲速;ρ為流體密度;νn為結(jié)構(gòu)濕表面法向振速,聲場總壓滿足邊界積分方程

(2)

式中:pI為入射聲壓;Ω1為Γ外部空間;Ω0為Γ內(nèi)部空間;x為場點;y為源點;G(x,y)為自由空間格林函數(shù),r=x-y。

G(r)=e-ikr/4πr,

(3)

(4)

式中,α為源點到場點的矢量r與法向n的夾角。結(jié)構(gòu)濕表面法向振速和法向位移存在關系vn=iωun,若把式(1d)代入式(2)得到:

(5)

p(x)-pIxinΩ1

(6)

式(5)和(6)描述了散射問題中聲壓與結(jié)構(gòu)響應的耦合關系。

1.2 邊界積分方程離散和流固耦合問題的解耦

為進一步推導和用于數(shù)值計算,需對式(5)和式(6)進行離散。將結(jié)構(gòu)濕表面離散為Ne個三角形單元,Ng個節(jié)點。整個曲面上的積分可表示為對所有單元上的積分求和的形式,式(5)寫為

(7)

當單元尺度相比于聲波長足夠小,可用單元平均壓力來表示單元上任一點壓力,單元平均法向位移表示單元上任意一點的法向位移。把p(y)、un(y)放到積分符號外面,單元j表面的壓力表示為

(8)

式(8)說明單元j上的壓力由結(jié)構(gòu)濕表面上所有單元的壓力和法向位移決定,式(8)寫為矩陣形式

(9)

E1{p}=C1{un}+{pI}

(10)

(11a)

(11b)

式(6)也可離散為式(12b)矩陣的形式,E2和C2的形式與式(11)類似。矩陣E1、E2、C1、C2僅與結(jié)構(gòu)濕表面的特征和頻率有關,聲散射問題中的流固耦合關系以及場點聲壓可用式(12)描述。式(12a)建立了表面聲壓與濕表面法向位移的關系,式(12b)說明聲場中任一點總聲壓{pt}可由表面聲壓和法向位移求解。

(12a)

{pt}-{pI}=E2{p}+C2{un}

(12b)

結(jié)構(gòu)振動問題同樣可以用矩陣形式表示,聯(lián)立式(12a)后得到散射問題的耦合方程

(13)

式中:{F}為外載荷;Ms、Bs、Ks分別為結(jié)構(gòu)質(zhì)量矩陣、阻尼矩陣和剛度矩陣;G為表面壓力{p}與節(jié)點力{Fp}等效轉(zhuǎn)換矩陣,滿足G{p}={Fp}。

引入節(jié)點位移與單元法向位移的轉(zhuǎn)化矩陣L,即{un}=L{u}。矩陣G與L具體形式見附錄A。式(13)中存在兩個獨立方程和兩個未知量,有先求解結(jié)構(gòu)位移{u}和先求解表面聲壓{p}兩種解耦途徑。這里采取先求解{u}的思路,消去變量{p}得到下式

[-ω2Ms+iωBs+Ks]{u}=

(14)

[-ω2(Ms+Ma)+iω(Bs+Ba)+Ks]{u}=

(15)

1.3 散射聲場的計算

根據(jù)式(15)求得結(jié)構(gòu)表面位移{u},進而用L求得法向位移{un}。將式(12a)代入(12b)得到:

(16)

總聲壓減去入射聲壓即為散射聲壓,{ps}={pt}-{pI}。

1.4 剛性散射與彈性散射的關系

一般而言,復雜的公式可退化為特殊情況下的簡單形式。1.1～1.3節(jié)討論的是彈性結(jié)構(gòu)聲散射問題。在剛性散射中不存在結(jié)構(gòu)響應,式(13)中{u}為0,退化為

(17)

式中,{prig}Γ為純剛性散射時,結(jié)構(gòu)表面的聲壓。遠場散射聲壓計算式(16)退化為

(18)

進一步將式(17)代入式(15)以及式(18)代入式(16)發(fā)現(xiàn)

[-ω2(Ms+Ma)+iω(Bs+Ba)+Ks]{u}=

{F}-G{prig}Γ

(19)

(20)

2 程序流程

(1) 建立結(jié)構(gòu)有限元模型并提取結(jié)構(gòu)濕表面為邊界元網(wǎng)格,同時繪制聲場點;

(2) Fortran程序A:讀取邊界元網(wǎng)格數(shù)據(jù)、聲場點數(shù)據(jù),計算轉(zhuǎn)換矩陣G和L,計算系數(shù)矩陣E1、C1、E2、C2;

(3) Fortran程序B:讀取邊界元網(wǎng)格數(shù)據(jù),通過編寫的平面波、球面波程序生成入射聲壓矩陣{pI};

(4) DMAP程序C:計算剛性散射表面聲壓矩陣{prig}Γ和剛性散射場點聲壓矩陣{prig};

(5) DMAP程序D:計算附加質(zhì)量矩陣Ma、附加阻尼矩陣Ba、附加壓力矩陣{Fa};

(6) Nastran 108求解序列:利用DMAP語言將Ma、Ba和{Fa}插入結(jié)構(gòu)頻響求解序列,計算完成后得到節(jié)點位移{u};

(7) DMAP程序E:根據(jù)式(20)計算二次輻射項和彈性散射聲場;

(8) 可視化:結(jié)構(gòu)響應結(jié)果已在步驟(6)生成,導入后處理軟件即可實現(xiàn)可視化,聲場結(jié)果修改為后處理軟件接受的格式后實現(xiàn)可視化。

當僅需要計算剛性散射時,執(zhí)行步驟(1)～(4)。

3 算例驗證:水中彈性球殼的散射聲場

通過MATLAB編程,得到球殼在平面波作用下的剛性散射與彈性散射理論解。采用本文方法,計算了1 Pa平面簡諧波作用下的剛性散射與彈性散射聲場。球殼內(nèi)部為真空、外部為水。其中水的密度為1 000 kg/m3,水中聲速取為1 450 m/s。球殼參數(shù)如表1,a為球殼半徑,h為球殼厚度,E為彈性模量,ν是泊松比,ρs是結(jié)構(gòu)密度。球殼被劃分為3 200個三角形單,元、邊界元網(wǎng)格與有限元網(wǎng)格一致,不考慮結(jié)構(gòu)阻尼,結(jié)構(gòu)有限元模型如圖1所示。為滿足遠場要求,在距離球心100 m的圓周上均勻設置360個聲場點。根據(jù)散射強度定義式(21)得到散射強度計算值。

圖1 球殼有限元模型Fig.1 The finite element model of spherical shell

表2 圓柱殼初始參數(shù)Tab.2 Initial parameters of cylindrical shell

(21)

首先設置球殼厚度為3 mm,計算了不同頻率下的剛性散射和彈性散射聲場。然后計算ka=2時,球殼厚度為3 mm和50 mm兩個不同厚度情況下的彈性散射聲場。如圖2所示,不同頻率下剛性散射數(shù)值結(jié)果與理論值的曲線完全吻合,證明了算法計算剛性散射的正確性。在不同頻率和球殼厚度下本文算法都能預報準確的散射強度和指向性,如圖3所示。背向散射理論值與計算值對比表明,當邊界元網(wǎng)格數(shù)目一定時,本文數(shù)值方法在較寬的頻段內(nèi)能得到較準確的結(jié)果,如圖4所示。

(a) ka=1

(a) h=3 mm,不同頻率

圖4 剛性球散射理論值與計算值(180°)Fig.4 Theoretical and numerical results of rigid sphere scattering (180°)

4 帶端蓋加肋圓柱殼散射聲場預報與特性分析

針對圓柱殼這一典型結(jié)構(gòu)的散射問題,Junger等[15-16]用解析方法研究了無限長、有限長以及加肋圓柱殼的散射問題。Zhen等和潘安等先后研究了雙層柱殼和加肋板雙層柱殼聲散射解析解。以上工作主要在解析基礎上用彈性波理論對圓柱殼散射進行分析,很少通過結(jié)構(gòu)響應說明散射特征,也未給出彈性散射影響因素的一般結(jié)論。復雜結(jié)構(gòu)難以推導解析解,數(shù)值計算基礎上的結(jié)構(gòu)響應分析對于散射特性研究和結(jié)構(gòu)聲學設計或是一種直接有效的方法。帶端蓋的有限長加肋圓柱殼更接近某些水下結(jié)構(gòu),殼體厚度、肋骨有無、自由液面位置等都可能是影響散射特征的因素。為進行相關研究,下面選取水下典型聲學試驗[17]中的圓柱殼數(shù)據(jù)作為初始參數(shù)進行數(shù)值仿真研究,如圖5所示。

圖5 加肋圓柱殼模型Fig.5 The model of cylindrical shell with stiffeners

4.1 結(jié)構(gòu)彈性的影響

在距離圓柱軸心100 m的圓周上,均勻布置360個聲場點,剖面如圖5所示。入射聲壓方位為180°,結(jié)果用無因次公式|pr/p0a|表示。其中,r為聲場點到圓柱殼軸心距離,p0為入射聲壓,p是場點聲壓。設定t/a=1%不變,無肋骨和自由液面,計算了ka=0.1、0.5、1、2、5、10共6個工況。

圖6表明,結(jié)構(gòu)的彈性作用對散射指向性和強度都不可忽略。ka=0.1時,剛性散射幾乎為零,說明在很低頻率且僅考慮剛性情況時聲波可以完全繞過圓柱,但彈性情況下二次輻射項的存在使散射回波強度明顯增大。在ka=0.5時,彈性作用使背向散射減小,回波強度幾乎為零。ka=1時,彈性作用對正向散射和背向散射都有增強作用。ka=2時,彈性作用使回波增強,正向散射減弱。圖6(e)、(f)表明,ka=5,10時二次輻射作用使正向散射增強,剛性散射對背向散射強度的貢獻較大。文獻[18]中認為彈性作用使?jié)撏Щ夭◤姸葴p弱的結(jié)論或有些片面,彈性作用不是單一使散射強度降低或者增強,具體與頻率有關,不能簡單認為結(jié)構(gòu)彈性吸收入射能量使回波強度減小。由于不同頻率下附加壓力激發(fā)的二次輻射項具有不同特征,與剛性項疊加后的彈性散射可能增強也可能減弱。

4.2 柱殼厚度及加肋骨對彈性散射影響

如表3所示,在ka=1情況下共計算了6個工況,并繪制圖7所示的散射圖樣。結(jié)果表明,隨著t/a增大,彈性散射與剛性散射圖樣仍有明顯不同,在較大的t/a范圍內(nèi)彈性作用仍不可忽略。工況4和6的結(jié)果表明,肋骨對散射強度和指向性影響較小。即使在工況5中肋骨厚度增加至2b,僅強度值略有減小。文獻[19]計算復合材料圓柱殼散射時也得到類似結(jié)果,但沒有給出肋骨作用不顯著的原因。圖8中結(jié)構(gòu)響應云圖顯示肋骨的變形主要是隨著柱殼平動,這與附加壓力在該頻率下激發(fā)的柱殼振動形式有關。此時柱殼響應以1階彎曲模態(tài)為主,主要發(fā)生軸向的彎曲變形。由于環(huán)肋骨對軸向彎曲剛度貢獻小而加厚殼體可以明顯提高彎曲剛度,增加殼體厚度比肋骨對散射的影響顯著。由于軸向的縱肋骨能夠增加抗彎剛度,可以推斷在該頻率下軸向的縱肋骨會對柱殼彈性散射影響較大。因此,肋骨對柱殼散射的影響與對應頻率下的振動形式和肋骨布置形式有關。

圖7 圓柱殼厚度和加肋骨對散射影響Fig.7 The influence of cylindrical shell thickness and stiffening to scattering

(a) 圓柱殼位移云圖

表3 圓柱殼彈性散射計算工況Tab.3 Calculation conditions for elastic scattering of cylindrical shells

4.3 自由液面的影響

在邊界元法中,可以通過修改格林函數(shù)形式來考慮自由液面的影響。如圖5,定義自由液面到圓柱殼上端蓋的距離為d。圖9結(jié)果表明,自由液面的存在對彈性散射和剛性散射的效果一致,僅強度值發(fā)生改變而散射指向性無明顯變化。圖10顯示了背向散射強度與d的關系。隨著結(jié)構(gòu)在水下位置變深,背向散射強度幅值在無自由液面的散射強度值附近起伏,變化幅度逐漸減小。當d大于150 m,|pr/p0a|變化幅度約在10%以內(nèi)。

圖9 自由液面對圓柱殼散射指向性的影響Fig.9 The influence of free liquid surface on the directivity of scattering from a cylindrical shell

圖10 自由液面位置對背向散射影響Fig.10 The influence of the position of the free liquid surface on the back-scattering

5 結(jié) 論

相比已有的聲散射計算方法,本文從積分方程法出發(fā)推導的散射問題解耦方程從機理上解釋了彈性散射聲場的組成和來源。編寫通用計算程序進行計算,與理論值對比證明了算法的正確性和把彈性散射劃分為剛性項和二次輻射項的合理性。這種分解厘清了彈性散射的物理機制,可用于進一步研究散射特性和調(diào)控機制。針對圓柱殼的散射聲場計算得到了一些有益的結(jié)果:在較大的殼體厚徑比范圍內(nèi),彈性作用均不能忽略。結(jié)構(gòu)彈性作用使圓柱殼散射強度與指向性發(fā)生較大改變,散射強度受彈性影響的改變量不是一致地增強或減弱,具體與頻率有關。肋骨有無對彈性散射影響與附加壓力激發(fā)的振動形式有關,在以1階模態(tài)為主的振動形式下,圓柱殼環(huán)形肋骨的有無對彈性散射指向性和強度影響很小,改變柱殼厚度使散射特性變化顯著。自由液面對指向性影響小,使散射強度值在無自由液面的散射強度附近隨深度振蕩且逐漸逼近。

附錄A

根據(jù)虛功原理,節(jié)點位移與單元位移存在如下關系

(A1)

un=un1+un2+un3

(A2)

(A3)

式中:uni為單元上節(jié)點i對單元法向位移的貢獻量;(xi,yi,zi)為節(jié)點的局部坐標。對于三角形單元,三個節(jié)點貢獻求和即為單元法向位移。一個節(jié)點共有三個平動自由度(u,v,w)和三個轉(zhuǎn)動自由度(θx,θy,θz),單元內(nèi)的平動和繞法向軸的轉(zhuǎn)動不引起法向位移。矩陣L為Ne行Ng列矩陣,每行由相關節(jié)點的[D]組裝而來。

利用做功等效原理,節(jié)點力和單元壓力存在關系

(A4)

(A5)