趙欣陽, 梅志遠, 祝 熠, 杜 度
(1.中國人民解放軍海軍工程大學 艦船與海洋學院,武漢 430000;2.中國人民解放軍92578部隊,北京 100000)
潛艇聲隱身技術一直備受關注,數十年減振降噪工程技術的發(fā)展使得各國在潛艇聲輻射控制領域取得了長足進步,嚴重削弱了傳統(tǒng)被動聲吶的探測性能。而主動聲吶具有探測低噪聲甚至“無聲”目標的優(yōu)勢,因而得到重視。20世紀70年代初,聲學覆蓋層的引入使得潛艇在聲目標強度控制技術方面有了顯著提高,但也帶來了一系列問題,此結構不似傳統(tǒng)金屬結構,其內部細觀特征復雜,材料性質多樣,以往的用于傳統(tǒng)金屬材料的聲目標強度預報方法已不再完全適用。此外,目前對于潛艇聲目標強度的研究大多集中在外殼方面,當然對于潛艇大部分結構來講,一般不需要考慮聲波透射進入內部的情況,但對于潛艇非水密結構卻必須區(qū)別對待。此結構整體被水所填充,僅有一部分升降通路為耐壓結構可視為剛性殼體。這意味著聲波入射到結構上時,外殼不可被視為剛性,必定會有大量聲波透射進入內部。由于潛艇內部存在擋板、肋骨以及各種管道,聲波便會在多種結構之間產生復雜的聲散射,散射聲場不可控,極有可能出現一些未知“亮點”, 所以對于這樣非水密結構的目標強度預報也顯得尤為重要。
聲目標強度預報的實質便是求解水下目標的聲散射問題,即求解目標在聲激勵作用下的散射聲場,目前該問題的解決方法主要分為理論與數值計算兩部分。理論研究方法主要有積分方程法[1]與簡正級數法[2]兩種,數值計算方法可分為以Helmholtz表面積分方程為基礎發(fā)展起來的邊界元法[3]、有限元結合邊界元法[4]、無限空間中的有限元法[5]、T矩陣方法[6]與波疊加方法[7]等等。這些方法在解決簡單結構方面各具優(yōu)勢,針對性也不盡相同,但均在處理含復雜內部結構的聲學覆蓋層時存在較大阻力。在處理此結構時,有限元法雖能對其進行建模,但在離散整個結構過程中會造成網格數量巨大,計算代價極高。為研究聲學覆蓋層的聲學特性,參數反演是一種很好的思路。Schneider等[8]將超聲波法用在測定固體彈性常數上,利用板中的縱波與橫波測量了幾種金屬的彈性常數。陶猛等[9]運用模擬退火算法,聯立多層平板結構反射系數計算傳遞矩陣,反演得到了消聲瓦結構的等效反射系數。孫鐵林等[10]基于遺傳算法,利用水下覆蓋粘彈性板的回波實驗,反演得到了復合板的相速度、衰減系數與彈性模量。金國梁等[11]利用遺傳算法對復雜聲學覆蓋層進行了參數反演以得到其等效參數,在利用等效參數進行圓柱殼體聲散射計算時發(fā)現等效后結構與原有復雜結構的聲散射形態(tài)函數基本吻合。周世根[12]分別采用等效剛度均勻化方法與參數反演法提取了內部含圓柱型空腔與圓錐形空腔的等效材料參數,研究結果表明,在利用等效材料進行聲目標強度計算時,殼體內部介質為空氣時利用此方法更好。
目前參數反演方法在對聲管試驗亦或是無限大平板仿真反演中基本能達到較好的效果,但在運用到圓柱殼亦或是其他結構上卻精度不高,這是由兩方面原因所造成的,一是所采用的優(yōu)化算法較為單一,存在一定的局限性;二是考慮不足,未考慮斜入射以及聲波透過結構后的多次散射問題。本文針對以上問題,結合現行不同算法的優(yōu)勢,在考慮斜入射、多次散射、吻合效應以及橫波波速閾值的情況下利用分層介質理論構建了雙層彈性介質層反演模型,此模型可替代原有敷設聲學覆蓋層的復雜模型,為非水密結構聲目標強度預報提供了一種有效的方法。
潛艇在結構設計上可分為耐壓結構與非耐壓結構,耐壓結構主要包括潛艇的主體結構,如耐壓船體、耐壓指揮室等等;非耐壓結構則可分為非耐壓水密結構與非耐壓非水密結構,前者主要包括主壓載水艙、燃油壓載水艙等,后者主要由上層建筑和指揮室圍殼等結構構成,本文重點討論的對象是非耐壓非水密結構。對潛艇來講,聲波從外部入射時,以空氣為背襯的耐壓殼結構使得外殼結構內外的特性阻抗失配,聲波幾乎無法透射到結構內部,更不會對內部結構產生影響。但對于圍殼與上建這樣的非水密結構,外殼內外特性阻抗一致,聲波透過外殼直接進入內部,在內部結構間形成復雜的聲散射,從而對整體的回波產生貢獻。
為真實模擬表面敷設聲學覆蓋層的潛艇非水密結構實際特征,同時提高仿真計算的效率,本文構建了如圖1所示的無限長圍殼模型。截面為橢圓型的圍殼浸沒在無限大水域中,a和b分別代表圍殼的長軸半徑與短軸半徑,圍殼內部由幾個圓柱形耐壓容器構成,外殼與耐壓容器間充滿水,在圖1所示的坐標系下,內部耐壓裝置的位置與尺寸信息如表1所示。同時,圍殼外殼材料為鋼,表面敷設了聲學覆蓋層,局部放大結構如圖1右側所示,h1代表鋼外殼厚度,h2代表聲學覆蓋層厚度,h3代表空腔底部到鋼外殼距離,h4代表空腔高度,d代表空腔直徑,圍殼及聲學覆蓋層的幾何參數如表2所示。
表1 內部耐壓裝置位置與尺寸信息Tab.1 Position and size information of internal voltage resistance equipment
表2 圍殼及聲學覆蓋層幾何參數Tab.2 Geometric parameters of submarine sails and acoustic coating
圖1 無限長圍殼模型(橫截面)Fig.1 Infinitely long submarine sails model (cross section view)
在本文所構建的模型基礎上,想要通過反演手段得到與原結構目標強度一致的等效結構就必須認識到非水密結構的特殊性,即聲波從聲學覆蓋層一側入射與從鋼板一側入射時兩者反射系數與透射系數存在差異。這時如果采用單層均質材料等效整個復雜結構則無法達到外側入射與內側入射反射系數與透射系數不一致的效果,這里引入考慮損耗的分層介質波理論解決此問題,通過構建雙層不同材料的彈性介質層達到聲波從兩側分別入射時反射系數與透射系數不一致的效果。
固態(tài)分層介質中的彈性波實際上就是一套遞推公式的運用,這其中考慮了縱波與橫波在固體中的折射與反射,以雙層固態(tài)分層介質為例進行理論推導。如圖2所示,當平面波從1側入射時, 利用場在分界面的連續(xù)性條件可以得到1、2界面與2、3界面、2、3界面與3、4界面的速度關系,如式(1)所示。
圖2 聲波在雙層彈性介質層中的反射與透射Fig.2 Reflection and transmission of sound waves in double elastic dielectric layers
(1)
式中,a11,…,a44由固體介質中的縱波波速c、縱波折射角θ、縱波在厚度上的相移P、橫波波速b、橫波折射角γ、橫波在厚度上的相移Q以及密度ρ組成。c與b可由楊氏模量E、密度ρ與泊松比σ決定,如式(2)所示。
(2)
當考慮材料損耗η后,楊氏模量就變?yōu)榱耸?/p>
E′=E(1+iη)
(3)
將層與層之間的速度關系建立聯系,可直接得到整個層系的速度關系
(4)
當介質1與4是液體時,切向應力消失,借助連續(xù)性條件便可得到層系的反射系數V與透射系數D如式(5)和式(6)所示。
(5)
(6)
當從介質4側入射時,只需調換矩陣[a(2)][a(3)]的位置,便可得到對應的反射系數與透射系數。為了模擬無限大雙層彈性介質層,在COMSOL中建立如圖3所示模型,第一層材料為聚氨酯聲學材料,第二層材料為PVC材料,厚度分別為40 mm、10 mm,兩層彈性介質層兩側流體為水,材料參數如表3所示。
表3 材料參數Tab.3 Material parameters
圖3 無限大雙層彈性介質有限元模型Fig.3 Finite element model of infinite double-layer elastic medium
水域兩端設置完美匹配層以模擬無限遠水域,在兩層彈性介質層與水域的上下表面與左右表面分別設置Floquet周期性條件用以模擬無限大雙層彈性介質層與水域。平面波以0°~80°入射到聚氨酯一側或PVC一側,入射頻率為1~5 kHz,步長為1 kHz。
圖4中的點為根據理論得到的結果,藍線表示仿真結果。從圖4可以看出,在考慮損耗后,利用分層介質波推導的聲學理論解與無限大雙層彈性介質層仿真解完全一致,也就意味著用于之后的聲學系數計算是可靠的。
(a) 反射系數理論解與仿真值對比(聚氨酯一側入射)
(b) 透射系數理論解與仿真值對比(聚氨酯一側入射)
(c) 反射系數理論解與仿真值對比(PVC一側入射)
(d) 透射系數理論解與仿真值對比(PVC一側入射)圖4 聲學系數理論解與仿真解的對比Fig.4 Comparison between theoretical solution and simulation solution of acoustic coefficient
目前聲學參數反演多采用遺傳算法尋找最優(yōu)解,遺傳算法的本質是一種并行、高效的全局搜索手段,其通過自適應尋優(yōu)手段得到最優(yōu)解?,F有的僅依靠遺傳算法反演聲學參數的手段基本能夠滿足正入射下簡單周期結構的求解但在個別頻率下仍存在較大誤差,這是因為遺傳算法在求解過程中重點關注全局尋找最優(yōu)解,對于局部重視不夠。為解決這部分存在的問題,考慮在遺傳算法中加入非線性規(guī)劃,將非線性規(guī)劃的局部強搜索能力與遺傳算法的全局強搜索能力結合在一起,從而提升尋找全局最優(yōu)解的能力。
遺傳算法的基礎步驟主要分為六步,依次為種群初始化、適應度函數選擇、選擇運算、交叉運算、變異運算與終止條件判斷,在本問題的尋優(yōu)過程中會存在陷入局部最優(yōu)解、效率不高、約束考慮不足等問題,對應這些問題,本文提出了對應的解決措施。
(1) 初始化種群的選擇,反演聲學參數的種群容易在模量的選擇上出現問題,如果模量的范圍上選取過大,譬如從104~1012,會使得模量在隨機的過程中傾向于設定的上限,從而使得尋優(yōu)過程中陷入次優(yōu)解而達不到目標,一種方法是主動降低模量設定的上限,使得模量在隨機的過程中更加均勻地分布在整個范圍空間,但是這種方法需要人為估計,具有一定的經驗性。另一種思路比較簡單,就是在隨機過程中將模量做一些預處理,將模量取對數后再進行隨機,隨機完成后再將預處理的模量取10次方得到模量的真實值進行求解。
(2) 適應度函數的選擇,本小節(jié)的目的是通過對主要參數的調整使得反演后的等效結構能夠達到復雜結構的聲學特征,表征上就是優(yōu)化后結構的反射系數與透射系數能夠達到與復雜結構的反射系數與透射系數一致的水平,所以可設定式(7)作為適應度函數,那么本節(jié)的目的就轉化成了求函數最小值的問題,函數值越小的個體,適應度值越小,個體越優(yōu)。
abs(abs(Vij)-Rij))
(7)
式中:Dij、Vij為考慮損耗后推導出來的雙層無限大彈性介質聲學系數;Tij、Rij為復雜結構仿真或試驗得到的聲學系數。i從1取到9,分別對應入射角0°~80°時的情況,j從1取到2,分別對應聲學覆蓋層一側入射與鋼板一側入射的情況。
(3) 非線性規(guī)劃,措施(1)與(2)是針對本文問題以及利用遺傳算法反演聲學參數過程中易存在的問題提出的解決措施,基本能夠滿足尋優(yōu)要求,但在個別頻率下仍存在較大誤差,這是因為遺傳算法在求解過程中重點關注全局尋找最優(yōu)解,對于局部重視不夠,所以需非線性規(guī)劃來加強局部最優(yōu)解的求取。要想在遺傳算法中引入非線性規(guī)劃,首先需要對適應度函數做處理,將適應度函數通過數學上的變換使其可以求導,這樣可以加快非線性規(guī)劃尋找極值的速度,那么適應度函數的形式就變?yōu)榱耸?8)所示。
(8)
非線性規(guī)劃的核心就是研究一個多元函數在無約束或者多約束下的極值問題,經典的非線性規(guī)劃算法采用梯度下降的方法求解,局部搜索能力強,且收斂速度加快。
(4) 約束考慮
為了保證反演過程的收斂性,遺傳算法與非線性規(guī)劃相結合的組合算法在種群初始化階段與進化階段均對參數設置了范圍約束,如將材料1和材料2的楊氏模量E范圍設為104~1012,泊松比σ設置為0.01~0.499,損耗因子η設置為0.01~4,密度ρ設置為1 000~10 000 kg/m3,反演介質1的厚度范圍為10~40 mm、反演介質2與反演介質1的總厚度為50 mm,這其中并不考慮等效后各個參數的真實意義。但若想仿真結果與理論結果完全吻合,還需做進一步的約束處理,本文重點考慮橫波臨界波速與吻合效應。
① 橫波臨界波速的判定
本文對種群材料參數的范圍設置較廣,目的是在較大范圍內使得適應度值盡可能小,從而尋得最優(yōu)解,但這樣的處理會使得兩種反演介質的橫波與縱波波速不可控,很可能出現波速極小的情況。以兩種波中速度較小的橫波為例,當其速度小于5 m/s時,在5 000 Hz的情況下,固體介質的網格尺寸就必須控制在0.2 mm以內才能滿足精度要求(聲學計算時網格尺寸需小于五分之一波長)。這樣的網格尺寸會使得仿真計算量激增,達不到大幅提高計算效率的目的,所以需對反演材料每個頻率下的橫波波速設定一個閾值。因此,本文利用式(9)對反演材料進行判定,僅保留符合條件(9)的種群,從而在滿足網格精度要求的基礎上最大化提高仿真計算效率。
(abs(b)/5f)≥5 mm
(9)
式中:b代表橫波波速;f代表頻率。
② 吻合效應臨界頻率的判定
借助分層介質波理論得到無限大平板的透射系數后,會存在一種較為特殊的情況,即如果板周圍液體的阻抗遠小于板的阻抗,那么在入射波沿板的徑跡速度與板中自由波的相速度一致時,則會發(fā)生全透射。利用分層介質波理論去求解無限大雙層彈性介質的全透射條件過于困難,這里從脹縮波與彎曲波的角度去理解這個現象。平面波入射到平板上時,當液體中的聲波在板上的投影與板的彎曲波吻合,聲波便會激發(fā)板的固有振動,使得結構的聲輻射能力大大增強,這種現象叫做吻合效應。當板在一定條件下產生了吻合效應,板將強烈振動,透射系數顯著變化,此時有限元仿真所得到的聲學系數將難以與理論解吻合,從而導致反演結果不準確。所以為了避免反演材料在所研究的頻段內產生吻合效應,本文引入了王鵬偉等[13]所提出的“等效薄板法”用以判定吻合效應的產生與否。
利用經典層合板理論得到如圖5所示的層合板彎曲剛度計算公式為
圖5 層合板彎曲剛度示意圖Fig.5 Schematic diagram of bending stiffness of laminated plate
(10)
程序主要流程圖如圖6所示。首先進行系統(tǒng)初始化,完成MCU,ESP8266,傳感器模塊以及電源模塊的初始化配置。ESP8266設置為STA模式,接入無線AP實現與主控制器的WIFI連接,將單片機置為LMP3模式,關閉MCLK與SMCLK時鐘,僅開啟ACLK作為串口模塊時鐘。當接收到開始采集指令后,觸發(fā)串口中斷,激活CPU,執(zhí)行中斷服務程序,進行數據采集與AD轉換與無線發(fā)送,數據發(fā)送完成后開啟定時器中斷,當120 s內沒有收到采集命令,結束本次采集,進入待機模式等待下一次采集。
(11)
式中:m為薄板的面密度;c為薄板兩側介質聲速。利用式(11)對每個頻率下的反演材料進行判定,僅保留設置頻率小于吻合效應臨界頻率的種群,從而使得反演材料在所研究的頻段無法產生吻合效應。
為體現預處理后加入非線性規(guī)劃的遺傳算法的優(yōu)勢,制定三種方案進行對比,方案一為遺傳算法,方案二為對模量預處理后的遺傳算法,方案三在方案二的基礎上加入了非線性規(guī)劃,此外,所有方案中均對吻合效應臨界頻率與橫波臨界波速進行判定,遺傳算法的參數選擇如表4所示。
表4 各方案中遺傳算法的主要參數設置Tab.4 Main parameter setting of genetic algorithm in each scheme
適應度函數中囊括復雜結構的仿真值,所以需要構建無限大樣品用以提取聲學系數,此部分依托于COMSOL有限元仿真平臺。利用含圓柱空腔的鋼背襯周期單元模擬無限大樣品結構,如圖6所示。
圖6 圓柱形空腔單元模型Fig.6 Cylindrical cavity element model
基體、鋼板與水的材料參數如表5所示,水域兩端設置完美匹配層以模擬無限遠水域,在結構與水域的上下表面與左右表面分別設置Floquet周期性條件用以模擬無限大樣品與水域。平面波以0°~80°入射到覆蓋層一側或鋼一側,入射頻率為1~10 kHz,步長為1 kHz。將有限元軟件中提取的聲學系數與三型方案反演的聲學系數展示圖7中。
(a) 反射系數(覆蓋層一側入射)
(b) 透射系數幅值(覆蓋層一側入射)
(c) 反射系數幅值(鋼板一側入射)
(d) 透射系數幅值(鋼板一側入射)圖7 各種方案反演解與復雜結構仿真解的對比Fig.7 Comparison between inversion solutions of various schemes and simulation solutions of complex structures
表5 材料參數Tab.5 Material parameters
從圖7中可以看出,在加入吻合效應臨界頻率與橫波臨界波速的判定后,三種方案均能在一定程度上反演出復雜結構的聲學特征,但明顯方案三具有顯著優(yōu)勢,反演出來的聲學系數無論是從覆蓋層一側入射來看還是從鋼板一側入射來看都與復雜結構仿真解基本吻合,方案二雖較方案一在反演準確性上有一定的提高,但在部分頻率與部分角度下仍存在較大偏差。此外,加入了非線性規(guī)劃后的遺傳算法可以大幅提高收斂速度,在收斂速度與求解結果上都優(yōu)于基本的遺傳算法。
為驗證有限元方法計算水下結構目標強度的有效性,將基于COMSOL軟件的數值仿真結果與無限長彈性圓柱簡正級數解進行對比。模型如圖8所示,主要參數包括:① 圓柱體半徑為25 mm,彈性模量E=7×1010,密度ρ=7 800 kg/m3,泊松σ=0.33。② 水域外邊界設完美匹配層,用于吸收入射波,模型無限大水域。平面聲波正入射,頻率為1~300 kHz,步長取1 kHz。根據入射波和散射波計算目標強度
圖8 無限長圓柱有限元模型Fig.8 Finite element model of an infinitely long cylindrical cylinder
(12)
式中:ps表示距離目標r處的散射聲壓;pi表示入射聲壓,測量點選在反向散射處。對于無限長彈性圓柱,在正入射情況下,其入射聲壓與散射聲壓為
Pi=P0expi(kx-ωt)
(13)
(14)
圖9 無限長鋁制圓柱體的目標強度Fig.9 The target strength of an infinitely long aluminum cylinder
為驗證本文構建的參數反演方法用于非水密結構聲目標強度預報的可行性,按照無限大圍殼模型構建了大小相同內部結構一致的反演模型,如圖10所示。此外,考慮到圍殼的介質環(huán)境與反演模型的介質環(huán)境需保持一致,計算時采用如圖11所示的三維有限元模型。在圖11所示的三維模型中,在水域與固體的邊界上均施加Floquet周期性條件以模擬無限長圍殼模型與無限大水域,在水域外側設置PML層來形成無反射邊界條件。
圖10 反演模型Fig.10 Inversion model
圖11 敷設聲學覆蓋層的圍殼三維模型Fig.11 Three dimensional model of submarine sails with acoustic coating
由仿真得到敷設聲學覆蓋層鋼板的聲學系數,利用方案三對結構進行等效參數反演,部分結果如圖12所示。從圖中可以看出,兩種反演介質的材料參數隨頻率無明顯的規(guī)律可循,總體來看,反演介質1的楊氏模量與損耗因子均較反演介質2穩(wěn)定。
(a) 楊氏模量(取對數)
(b) 損耗因子
(c) 泊松比圖12 雙層彈性介質層參數反演結果Fig.12 Inversion results of parameters of double-layer elastic medium
為充分驗證本文構建的反演模型的準確性,在真實模型與反演模型的對比階段選取了x方向入射、y方向入射與45°方向入射三種入射角度,每種入射角度拾取反向、30°兩個方向的聲目標強度值用以對比,入射頻率從100 Hz到5 kHz,步長為100 Hz。
如圖13所示,紅線代表原結構的計算結果、黑線代表等效結構的計算結果。從圖中可以看出無論從哪個方向入射,在收發(fā)合置的情況下,原結構與等效結構的各個頻率下的目標強度值基本一致,僅在1 000~2 000 Hz頻段內等效結構的目標強度較原結構有些許浮動,x方向入射、y方向入射與45°方向入射時研究頻段的整體誤差分別為1.25 dB,1.19 dB與1.28 dB。在30°收發(fā)分置的情況下,原結構與等效結構的各個頻率下的目標強度值波峰波谷的趨勢相同,x方向入射、y方向入射與45°方向入射時研究頻段的整體誤差分別為2 dB,1.30 dB與1.75 dB。
(a) 收發(fā)合置(x方向入射)
(b) 30°收發(fā)分置(x方向入射)
(c) 收發(fā)合置(y方向入射)
(d) 30°收發(fā)分置(y方向入射)
(e) 收發(fā)合置(45°方向入射)
(f) 30°收發(fā)分置(45°方向入射)圖13 目標強度計算結果對比Fig.13 Comparison of target strength calculation results
同時,圖14給出了幾個頻率下原結構與等效結構的目標強度指向性對比。從圖14(a)~14(b)可以看出,低頻段等效結構與原結構指向性基本一致,但幅值存在一定差異。從圖14(c)~14(e)可以看出,在1 000~1 600 Hz頻段,等效結構與原結構在-30°~30°指向性一致,但幅值存在差異,尤其是1 200 Hz時,差異達到5 dB左右,這與此頻率下反演結果與原結構吻合程度不高有關。從圖14(f)~14(i)可以看出,隨著頻率的增加,原結構與等效結構目標強度一致性越來越高,且在-30°~30°范圍內,此頻段的目標強度基本一致??傮w來看,-30°~30°范圍內,除個別頻率外,原結構與等效結構目標強度基本吻合。
(a) 100 Hz
(b) 600 Hz
(c) 1 000 Hz
(d) 1 200 Hz
(e) 1 600 Hz
(f) 2 000 Hz
(g) 3 000 Hz
(h) 4 000 Hz
(i) 5 000 Hz圖14 目標強度指向性對比Fig.14 Comparison of target target strength directivity
本文利用參數反演方法,對內部含耐壓圓柱體的非水密結構的聲目標強度進行了預報。一方面,在已有的遺傳算法反演聲學參數的基礎上加入了非線性規(guī)劃,在算法中充分考慮了斜入射、多次散射、吻合效應以及橫波波速閾值的影響,并以此建立了雙層彈性介質層反演模型,該反演模型與原結構在各個角度下均具有基本相同的反射系數與透射系數。另一方面,本文建立了兩種內部含耐壓圓柱體的非水密結構用以對比原結構與等效結構的目標強度,從整個頻段目標強度對比結果來看,收發(fā)合置情況下,原結構與等效結構的各個頻率下的目標強度值基本一致,各個方向入射時的研究頻段整體誤差在1.5 dB以下;30°收發(fā)分置情況下,原結構與等效結構的各個頻率下的目標強度值波峰波谷的趨勢相同,各個方向入射時的研究頻段整體誤差在2 dB以下。從各個頻率下的目標強度指向性來看,中高頻段,原結構與等效結構目標強度吻合程度較好;中低頻段,原結構與等效結構目標強度指向性基本一致,但幅值存在差異;總體來看,-30°~30°范圍內,除個別頻率外,原結構與等效結構目標強度基本吻合。
綜上所述,本文所構建的參數反演方法在解決外殼敷設聲學覆蓋層,內部存在結構物的非水密結構聲目標強度預報方面具有一定優(yōu)勢。