賴喜悅,陳美霞,董文凱
(華中科技大學(xué) 船舶與海洋工程學(xué)院,湖北 武漢 430074)
裝有聲吶設(shè)備的艦艇艏部透聲窗是水下反潛作戰(zhàn)的探測窗口[1]。透聲窗內(nèi)部聲吶設(shè)備工作時,基陣發(fā)射的聲波經(jīng)過導(dǎo)流罩透射到外部聲場,由于透聲窗的遮擋和導(dǎo)流罩內(nèi)部結(jié)構(gòu)的透聲損失,波束會發(fā)生畸變,并伴隨著能量損失。因此,透聲窗設(shè)計在滿足結(jié)構(gòu)強(qiáng)度、阻力和耐用性等要求的前提下,還需要著重考慮其聲學(xué)特性[2]。
因透聲窗遮擋作用帶來的聲波能量損失被定義為遮擋損失,是透聲窗聲學(xué)特性設(shè)計關(guān)注的重點之一。然而,由于透聲窗結(jié)構(gòu)比較復(fù)雜且模型尺寸一般較大,加之工程中關(guān)注的頻率范圍并不局限于低頻,透聲窗遮擋損失的準(zhǔn)確計算較為困難。因此,在許多文獻(xiàn)中往往通過對平面波入射下單層或多層介質(zhì)模型的透聲性能進(jìn)行研究計算來粗略預(yù)報透聲窗的聲學(xué)性能。LEE 等[3]利用三層和四層介質(zhì)模型對復(fù)合透聲窗材料在一定頻率下的插入損失進(jìn)行了理論計算,并開展實驗驗證了計算結(jié)果的準(zhǔn)確性。童暉等[4]利用三層介質(zhì)模型對透聲窗材料在高頻下的聲學(xué)特性進(jìn)行了研究,并探討了材料屬性和平面波入射角度對其插入損失的影響。易燕等[5]提出了一種低頻寬帶測試技術(shù),解決了雙層鈦合金導(dǎo)流罩平板模型低頻段和大角度測試難度問題,得到了模型試樣低頻寬帶范圍內(nèi)插入損失的頻率譜和角度譜。李東升等[6]利用夾芯透聲窗的簡化模型研究了夾芯透聲窗參數(shù)對其聲學(xué)性能的影響,并對比了不同頻段下夾芯透聲窗與單層玻璃鋼透聲窗透聲性能的差異。此外,李源源[7]等建造了實際艦艇球鼻艏導(dǎo)流罩模型,實驗得到了球鼻首的全向和定向發(fā)射透聲損失。
綜上,對透聲窗聲學(xué)性能進(jìn)行準(zhǔn)確預(yù)報是對其進(jìn)行聲學(xué)設(shè)計的前提,而簡化的單層或多層介質(zhì)模型無法考慮實際透聲窗的復(fù)雜結(jié)構(gòu)、復(fù)雜邊界和復(fù)雜聲源,只能粗略地估算透聲窗的聲學(xué)性能。實驗方法研究透聲窗模型的聲學(xué)性能又會帶來較大的實驗成本,因此,提出一種準(zhǔn)確預(yù)報透聲窗遮擋損失的計算方法具有切實意義。
本文提出了一種解析–數(shù)值混合方法,用于快速預(yù)報透聲窗的遮擋損失。首先,利用球諧波展開法得到了單活塞聲源的聲場分布,并通過疊加原理得到陣列聲源的聲場分布,進(jìn)而根據(jù)透聲窗有限元模型各網(wǎng)格節(jié)點的坐標(biāo)得到模型表面任意一點由聲源激勵引起的聲壓和法向振速。將之插值回COMSOL 有限元模型之中,調(diào)用 COMSOL 的Kirchhoff-Helmholtz 模塊通過再輻射的方法求解模型的遮擋損失,最后建立簡單二維透聲窗模型,驗證了方法的有效性。具體的計算流程如圖1所示。
圖1 解析–數(shù)值混合方法計算流程圖Fig.1 Flowchart of the calculation process of analysis-numerical hybrid method
要想對透聲窗模型遮擋損失進(jìn)行準(zhǔn)確計算,聲源的建模和求解是第一步應(yīng)該做的內(nèi)容。圓柱形基陣是目前國內(nèi)外艦殼聲吶最普遍的一種基陣形式,受水面艦吃水深度的限制,通常采用收發(fā)合置的聲基陣,使得這種基陣形式對水平 360°范圍具有相同的探測能力。實際研究中可以將這種基陣簡化為圓面活塞陣列,在定向發(fā)射時,可以根據(jù)相控規(guī)則調(diào)控各活塞振元的振動相位,使得各個活塞輻射聲場波陣面重疊,實現(xiàn)波束聚焦都效果[8]。
文獻(xiàn)[9]和[10]采用球諧波展開法求解了圓盤形輻射體在聲場中的聲輻射,對于二維模型來說,陣列的每一個振元都可以簡化為以一定振速在自由空間中簡諧振動的線段,線段在y>0 方向的法向振速為V0,在y<0 方向的法向振速幅值為–V0,活塞半徑為a,r和θ分別為測點半徑及其與y軸正方向的夾角,模型示意圖如圖2所示。
圖2 單活塞聲場分布示意圖Fig.2 Schematic diagram of sound field distribution of single piston
將無限水域劃分為子域Ⅰ(0≤r≤a,0≤θ≤π /2)、Ⅱ(0≤r≤a,π /2≤θ≤π)、Ⅲ(r≥a,0≤θ≤π)后,各子域聲壓滿足時間簡諧形式如式(1)所示:
式中:ω為振動圓頻率;j為虛數(shù)單位;μ=1,2,3。
忽略時間項exp(-jωt),各子域聲壓還應(yīng)滿足如下邊界條件:Z0為介質(zhì)特性阻抗;k為流體中的波數(shù)。
從而,各子域的聲壓函數(shù)可以假設(shè)為式(3)中的形式:
式中:An,Bn和Cn為未知系數(shù);J和H分別表示一類貝塞爾函數(shù)和漢克爾函數(shù),H的下標(biāo)表示其階數(shù),上標(biāo)表示其為第一類漢克爾函數(shù)。
各子域之間需要滿足的聲壓–振速連續(xù)性條件如下:
將各子域假設(shè)的聲壓函數(shù)代入到上述連續(xù)性條件中,可以得到式(5):
式(5)的第2、3 個方程左右兩邊同時乘以cos(mθ),m和n均為自然數(shù),并在兩方程θ的取值區(qū)間內(nèi)積分,可以得到式(6):
其中分部積分可求得:
式中,pm(nm,n)=qmn(m,n)=,m=n,且m,n為奇數(shù)。
從而可知pmn和qmn有下面的性質(zhì),
取m=2s+1,s為自然數(shù),將式(6)中的2 個方程相加,利用上述pmn和qmn的性質(zhì)可以得到式(9):
式(9)等式兩邊同時在r=a處對r求導(dǎo),可以得到式(10):
聯(lián)立式(9)和式(10),取截斷數(shù)s=n=N,即可以得到一個由2*(N+1)個方程組成的線性方程組,求解線性方程組,得到聲壓函數(shù)中的未知系數(shù)An,Bn和Cn,回代到式(3),即可得到單活塞輻射聲場。
多個活塞組成陣列時,陣列半徑為Rr,活塞個數(shù)為2*Ns+1,即以第Ns+1 個活塞為中心,左右兩邊各有Ns個活塞,活塞間隔圓心角為δθ,第i個活塞與中心活塞的夾角為φi,陣列發(fā)射方向與中心活塞的夾角為θr。為了實現(xiàn)定向發(fā)射的功能,不同活塞振動在時域上需要存在一定的時延關(guān)系,這種關(guān)系可以通過調(diào)節(jié)頻域上的相位差來實現(xiàn),即不同活塞的速度振幅可以由如下形式表出:
式中,V0i表示第i個活塞的法向振動速度幅值。
忽略活塞之間的相互作用,陣列輻射聲場可以看作是單活塞聲場疊加而成,從而總聲場可以表示為
式中:TR,Tθ表示測點坐標(biāo);ri,θi表示測點在第i個活塞局部坐標(biāo)系下的坐標(biāo)。如圖3所示,對于平面上任意一點(TR,Tθ),我們都可以由余弦定理將其在第i個活塞局部坐標(biāo)系下的坐標(biāo)(ri,θi)表示成式(13)的形式:
圖3 局部坐標(biāo)系下測點坐標(biāo)示意圖Fig.3 Schematic diagram of measurement point coordinates in local coordinate system
獲取了陣列輻射聲場任意一點的聲壓P(TR,Tθ)后,結(jié)合透聲窗艙室?guī)缀涡畔?,就可以獲取透聲窗及艙室表面任意一點處的聲壓與法向振速。
以陣列中心為全局坐標(biāo)系的原點建立透聲窗艙室輪廓模型,對于透聲窗艙室表面上任意一點,其坐標(biāo)向量可以設(shè)為=[x,y]T,該點處的外法向單位矢量為=[nx,ny]T,從而該點處由陣列聲源激勵起的法向振速vn()可以表示為式(14)的形式,式中dΔ 為一較小的值,可取對應(yīng)頻率下流場聲波波長的1/20。
依據(jù)惠更斯原理,文獻(xiàn)[11]提出了一種再輻射的方法來計算透聲窗內(nèi)聲吶基陣的遠(yuǎn)場輻射指向性,這種再輻射的方法很好的適應(yīng)了實際透聲窗內(nèi)的聲場分布情況,運算比較簡單。結(jié)合這一方法,可以將遮擋損失的計算分為2 步:第一步是計算陣列在自由場中的遠(yuǎn)場聲壓級指向性,第二步是計算透聲窗遮擋后遠(yuǎn)場聲壓級指向性,兩者之差即為遮擋損失。
依據(jù)圖1,獲取表面振速信息后,調(diào)用有限元軟件COMSOL Multiphysics 6.0 的Pressure Acoustics,Kirchhoff-Helmholtz 模塊(壓力聲學(xué),基爾霍夫–亥姆霍茲模塊)進(jìn)行再輻射計算。通過輸入振速,該模塊基于基爾霍夫–亥姆霍茲積分公式對高頻輻射問題進(jìn)行計算,而無需對周圍流體進(jìn)行建模,可以極大地降低有限元模型的復(fù)雜度,使得計算較為迅速。調(diào)用Kirchhoff-Helmholtz 模塊后,在COMSOL中調(diào)用MATLAB 插值函數(shù)將振速輸入到幾何和網(wǎng)格表面,得到再輻射源。計算無遮擋聲場源特性時,需要將整個艙室外表面選中,給每個節(jié)點插值對應(yīng)的法向振速,而計算遮擋后的輻射聲場時,只需要將透聲表面選中并賦予振速計算。最后,只需要在后處理中將2 次計算結(jié)果相減,即可得到模型的遮擋損失。
在COMSOL 中建立二維單活塞模型,活塞半徑為35.5 mm,計算頻率為10 kHz,活塞法向速度幅值為1 m/s。有限元模型如圖4所示,其中內(nèi)圓的水平直徑即為活塞,內(nèi)圓和中間圓為水域,水的密度為1 000 kg/m3,水中聲速為1 500 m/s,外圓為PML 層。理論與有限元計算得到的該頻率下距離圓心r=1 000 m 處一周的聲壓如圖5所示??梢钥吹絻烧叩挠嬎憬Y(jié)果吻合較好,說明球諧波展開法能較準(zhǔn)確地求解單活塞在水中的輻射聲場。
圖4 單活塞有限元模型Fig.4 Finite element model of single piston
圖5 單活塞r=1 000 m 指向性圖Fig.5 Directivity pattern of single piston at r=1 000 m
在COMSOL 中建立活塞陣列模型,NS=10,各活塞的直徑均為35.5 mm,活塞法向速度幅值為1 m/s,陣列半徑為1.2 m,活塞間隔圓心角為3.5°,θr=0°,計算頻率為10 kHz,有限元模型示意圖如圖6所示,理論與解析計算得到的該頻率下距離圓心r=1 000 m 處一周的聲壓如圖7所示。也可以看到理論結(jié)果能與有限元結(jié)果吻合較好,陣列定向發(fā)射時,能量主要集中在發(fā)射方向。
圖6 活塞陣列有限元模型示意圖Fig.6 Schematic diagram of finite element model for piston array
圖7 活塞陣列r=1 000 m 指向性圖Fig.7 Directivity pattern of piston array at r=1 000 m
為驗證本文提出方法的準(zhǔn)確性,下面將建立一個簡單的二維模型,分別采用本文方法與有限元法對模型的遮擋損失進(jìn)行計算,并對結(jié)果進(jìn)行對比驗證。
在計算時,取陣列發(fā)射角度為0°,通過改變陣列發(fā)射方向與透聲窗中軸的夾角θx,可以在掃描扇面內(nèi)進(jìn)行周向掃描,取遮擋前后陣列平面內(nèi)沿著陣列發(fā)射方向距陣列中心r=1 000 m 處的場點聲壓級之差,為模型在θx角度時的遮擋損失,記為SL(θx),即
式中:LP(θx)為遮擋前的場點聲壓級,LP'(θx)為遮擋后的場點聲壓級。以θx為周向坐標(biāo),SL(θx)為徑向坐標(biāo),繪制極坐標(biāo)曲線,就可以表征透聲窗全向的遮擋損失特性。
以一橢圓模型為驗證模型,橢圓長半徑為2.3 m,短半徑為1.8 m,橢圓中心在原點,截取長邊頂點兩邊各5/16 橢圓弧長作為透聲窗,其余弧為吸聲側(cè)壁驗證模型的陣列參數(shù)與3.2 節(jié)一致,全局坐標(biāo)系原點與陣列中心重合。模型示意圖如圖8所示。
圖8 驗證模型示意圖Fig.8 Schematic diagram of verification model
分別采用純有限元和解析–數(shù)值混合方法得到的模型全向遮擋損失計算結(jié)果如圖9所示,可以看到,本文提出的方法準(zhǔn)確性較高,且相較于有限元,該方法無需建立出完美匹配邊界和復(fù)雜的陣列活塞聲源,模型簡單,計算效率較高。
圖9 模型全向遮擋損失結(jié)果對比Fig.9 Comparison of results of model omnidirectional obstruction loss
本文提出了一種用于透聲窗遮擋損失計算的解析–數(shù)值混合方法,并通過二維簡單模型進(jìn)行了方法驗證,得到了以下結(jié)論。
1)球諧波展開法可以有效求解二維活塞聲場,結(jié)果與有限元吻合較好,本文對各子域的聲壓函數(shù)假設(shè)合理。
2)不考慮活塞間的相互作用,運用疊加原理,可以較準(zhǔn)確的計算出活塞陣列的聲場分布,計算結(jié)果與有限元吻合較好,陣列定向發(fā)射時,聲能量主要集中在發(fā)射方向上。
3)本文提出的解析–數(shù)值混合方法可以較準(zhǔn)確的預(yù)報二維透聲窗模型的全向遮擋損失,思路清晰,計算簡單,相較于純有限元方法,不涉及無限域的模擬,也不需要建出活塞陣列(純有限元模型一般在活塞陣列附近網(wǎng)格質(zhì)量要求較高),對硬件要求較低且耗時較短。
后續(xù)可以依據(jù)這一思路,繼續(xù)將該方法拓展至三維模型,方法快速性的優(yōu)點將得到更充分的體現(xiàn)。