波浪作用下水下航行體出水成功概率預報*

孫龍泉，任澤宇，李志鵬，馬貴輝，趙建洲

(哈爾濱工程大學 船舶工程學院， 黑龍江 哈爾濱 150001)

潛射導彈發射技術中，對潛射導彈在不同的海浪環境下出水姿態參數的預報至關重要。潛射導彈點火技術主要分為兩類：①導彈發射一段時間后在水中點火；②導彈出筒后依靠慣性在水中航行，離開水面時點火。對于水面點火的水下發射方式，如果出水姿態參數不夠理想，導彈出水后的二次推進就無法調整彈體的軌跡，甚至會影響到潛射導彈的成功出水，因此預報在波浪作用下導彈出水的姿態參數對導彈的成功發射有著至關重要的作用。而實際中波浪的運動是非常復雜的，不同浪向、頻率的波浪相互疊加造成了波浪屬性的隨機性，所以實現實測波浪下水下航行體出水成功的預報更有實際的意義。

水下航行體出水過程中受到波浪的擾動，二者之間的作用屬于流體與固體之前的相互耦合作用，流固耦合數學模型主要包括理想勢流模型和NS模型。由于NS模型具有耗時以及大幅度運動過程中網格易發生嚴重變形等缺點，本文的計算模型采用理想勢流模型。倪寶玉[1]、Longuet-higgins等[2]就理想勢流模型的可行性做出了詳細的論述和驗證，使在忽略黏性等因素后的數學模型得到了大大簡化，提高了計算效率。對于理想勢流模型的求解，目前的主要求解方法包括邊界元方法和有限元方法。李帥[3]、Li[4-5]、Kashiwagi[6]等使用邊界元方法求解，Eatock[7]、Ma[8]等采用有限元方法求解。

然而，在水下航行體出水計算過程中，結構濕表面網格的重新布置以及自由面非線性邊界條件的處理都顯得十分復雜。為了解決這一系列復雜問題，Ma[8]和Yan[9]等使用QALE-FEM方法，在整個出水計算過程中只生成一次復雜的結構化網格，在下一時間步直接計算到邊界位置，大量節省了網格重置的時間。倪寶玉[1]提出在水下航行體出水過程中，由于其高速運動，流場可以近似認為是準靜態流場，可以繞射理論為基礎，采用勢流理論計算方法解決水下航行體出水問題。Xue[10]、Miloh[11]等就自由面非線性邊界條件的問題，進行了計算模型和邊界條件的簡化，并與試驗對比進行了有效性驗證。

在水下航行體彈道計算方面，李代金[12]、李體方[13]等對航行體發射過程中彈道的數學模型進行了建立和分析，這些數學模型可用于潛射導彈等水下航行體的彈道計算，也為其他彈道模型的建立提供了參考依據。劉曜[14]、姜濤[15]通過解析法計算了水下航行體受到波浪作用后所受載荷大小，并且通過軌跡方程得到了水下航行體運動姿態變化。王紅萍等[16]以自建的高精度數值波浪水池為基礎，探究了波浪擾動力以及運動姿態參數的變化規律，為后續的彈道設計打下基礎。尤天慶等[17]針對航行體尾空泡對流體阻尼力影響的問題，通過求解雷諾平均的斯托克斯方程組，進行阻尼力計算研究，驗證了彈道設計中考慮尾空泡的重要性。后續研究發現在航行體表面進行通氣可以改善航行體表面的流體動力特征，從而提高水下航行體的出水成功概率。例如：劉濤濤[18]、王克林[19]等對氣孔排氣的氣液兩相流特性進行了研究，為排氣孔的設計優化提供了參考；馬貴輝[20]、張耐民[21]就通氣空泡技術改善航行體水下彈道及出水姿態參數的有效性和可行性做了詳細論述。

綜上所述，雖然國內外的科研人員對水下航行體出水問題已經做了很多研究，但是大多需要借助計算流體力學(computational fluid dynamics, CFD)相關的商業軟件，計算時間長、工作量大，沒有辦法形成規模計算，對海浪和水下航行體相互作用的定性分析也存在困難,同時，缺少在實時海況下的水下航行體出水成功概率的快速預報，在此方面無法為潛艇作戰提供幫助。本文以邊界元方法，在波浪和水下航行體的流固耦合分析上做了一定程度的改進，大大縮短了計算時間，并且將計算結果納入規則波浪數據中，以此為基礎，能夠在一定程度上實現實時海況下的航行體出水成功概率預報。

1 理論與方法

航行體在水下運動時，所受到的重力、浮力、流體力會影響水下航行體的原始軌跡并發生偏離。出水姿態參數是影響航行體出水成功與否的重要因素之一，如果出水彈道設計不好，潛射導彈出水后具有的角速度或偏轉角可能會超出允許的控制范圍，無法通過水面點火等方式調節平衡，從而直接導致發射失敗[22]。本文主要以俯仰角為例，研究波浪參數對尾部觸水俯仰角的影響，并且以尾部觸水俯仰角為判定依據進行水下航行體出水成功概率的預報。水下航行體運動姿態的示意如圖1所示，俯仰角φ為航行體軸線和水平方向的夾角。

圖1 水下航行體運動軌跡Fig.1 Trajectories of underwater vehicles

如圖1所示，小尺度航行體在水中高速運動時，雷諾系數很大，慣性力的影響遠大于黏性力的影響[1]，所以對水下航行體的數學模型做簡化，能夠實現規模計算。彈體在水下運動的數學模型簡圖如圖2所示，航行體發射時有垂直向上的初始速度，同時加入了水平艇速的影響，上升過程中主要受到重力、流場慣性力和波浪力的作用。

圖2 彈體出水坐標系Fig.2 Projectile body outlet coordinate system

如圖3所示，以半球頭回轉體作為水下航行體，其中圓柱體長L0，直徑D0，距頭部頂端0.56L0處為質心。在航行體表面總速度勢中計入波浪入射勢，以考慮波浪對水下運動航行體的影響，這樣可簡化自由液面處網格的劃分，將實際物理模型(圖1)簡化為邊界元計算模型(如圖3(a)所示)。模型中，航行體發射水深為30 m，垂向速度為25 m/s，具有1.5 m/s的向彈體坐標系x軸負方向(艇速方向)的初始艇速。因為模型中不涉及出筒階段，為考慮航行體水下發射實際情況，在數值模型中設置水下航行體初始的俯仰角速度為ω0，初始俯仰角為φ0。

(a) 航行體網格(a) Grid of the vehicle

如圖3(a)所示，水下航行體結構表面用三角形網格進行劃分，計算過程中將其視為剛體。模型中，設水下航行體在水中運動時為完全浸濕狀態，忽略空泡等強非線性條件影響。

1.1 控制方程

在無旋、無黏假設下的理想勢流理論中，航行體總速度勢為：

Φ=φ+φi

(1)

式中，φ為結構表面的擾動勢，φi為波浪入射勢。

航行體總速度勢滿足拉普拉斯方程：

?2Φ=0

(2)

結構表面滿足不可穿透條件：

(3)

式中，n為物面法向向量，V為速度向量。

流場無窮遠處滿足：

?Φ→?φi

(4)

流場域內任意一點的速度勢都可以用其邊界上的速度勢及其法向導數來確定[23]，所以需要采用邊界積分方法[24]求解流場擾動勢φ：

(5)

式中，S為航行體表面邊界面，p和q分別為場點和源點。自由空間的格林函數G(p,q)=1/|R-r|，R和r分別是場點和源點的位置矢量。

將求解得出的速度勢代入非定常的伯努利方程中，得到式(6)，將由式(5)求得的擾動勢φ代入式(6)中求得水下航行體運動過程中所受壓力，同時計及流場靜壓和大氣壓力。

(6)

式中：P0為大氣壓，水下航行體出水時其表面將只承受大氣壓力；h為距離水面的高度；ρ為海水密度。

水下航行體的結構動力學方程如式(7)所示：

(7)

式中，Ms為結構質量矩陣，Ks為剛度矩陣，fs為結構外力，Isf為聯系結構單元與流體邊界單元的轉換矩陣，Af為單元面積矩陣。

將水下航行體表面新的單元、節點位置信息輸入到流場動壓力計算程序之中，計算在新時間步內結構表面各單元所受壓力大小，然后更新節點、單元信息，循環計算到要求時間步長為止。

1.2 波浪入射勢

本文計算工況中，波高相對于水深為小量，所以可以采用小振幅波理論對水下航行體在規則波浪中的運動進行分析。小振幅波模型如圖4所示，圖中，λ為波長，H為波高，A為波幅。

圖4 小振幅波波面Fig.4 Small amplitude wave surface

小幅波面方程如式(8)所示：

η(x,t)=Acos(kx-ωt+ε)

(8)

式中，ω代表波浪頻率，k為波數，ε為初相位。

小振幅波理論中，水質點運動緩慢，故自由表面非線性運動和動力邊界條件可以簡化為線性邊界條件。因此，波浪向前傳播時，波浪的入射勢為：

(9)

若波浪的傳播方向和坐標系之間存在夾角γ，則可以通過坐標轉換得到新的波浪數學模型：

(10)

當艇速水平運動方向和波浪傳播方向相同時為順浪，此時γ=180°；當方向相反時為逆浪，此時γ=0°。

1.3 有效性驗證

為驗證算法的有效性，將在開放水池所測得的試驗結果與數值模擬結果對比。試驗工況為：初始艇速3.78ω0L0，垂直速度63.01ω0L0，波面狀態為靜水。

(11)

將計算所得的水下航行體頭部觸水俯仰角和尾部觸水俯仰角，同試驗結果進行對比，來驗證算法的有效性。對比結果如表1所示，頭部觸水俯仰角和尾部觸水俯仰角的相對誤差較小，可以認為在姿態計算上算法有效。

表1 水下航行體出水俯仰角對比

為驗證波浪力的有效性，采用海洋工程中常用的經驗公式進行比較。對于與波長相比尺度較小的細長彈體(彈體D/λ<0.2)的波浪力計算，在工程設計中廣泛采用莫里森公式[25]。該公式假定，柱體的存在對波浪運動無顯著影響，認為波浪對柱體的作用主要是黏滯效應和附加質量效應。

(12)

由于莫里森公式一般應用于靜態或者準靜態的結構物波浪力計算當中，所以本文在校核勢流方法計算結果有效性時，假定柱體靜止于水面附近，然后分別用莫里森公式和勢流方法計算柱體單位長度所受波浪力大小，計算結果如表2所示。在水面處，應用勢流方法比應用莫里森公式計算相比偏小，3級浪和5級浪誤差分別為6.71%和9.56%，可認為計算結果符合波浪力有效性要求。結合上述對彈體出水俯仰角的驗證結果，本文所建立的基于勢流理論的航行體出水姿態計算模型符合有效性的要求，可以應用到本文的計算當中。

表2 莫里森公式和勢流方法計算的單位長度波浪力對比

、

2 實時波浪預報方法

2.1 實際海況分析

實際海浪的波面運動是一種十分復雜的自然現象，受到自身重力、海面附近風速、風向、水面以下地形地貌特點，甚至是地球自轉偏向力、潮汐作用等諸多復雜因素的影響，從而使其波動規律表現出極大的隨機性，增加了對水面附近艦船和水下航行體的運動進行準確預測的困難程度。圖5所示為渤海灣某一海域在臺風下的30 min波面時歷曲線,本文以該波浪作為實測海況進行水下航行體的出水姿態成功概率預報。

圖5 波面30 min的時歷曲線Fig.5 30-minute time-history curve of the wave surface

圖6 實測波浪的波高統計直方圖Fig.6 Statistical histogram of wave height of measured waves

為驗證實測波高的分布規律，將實測波高進行累計概率統計，得到如圖7所示的實測波高的概率分布曲線。同時為得到理論上的波高概率分布曲線，參考式(13)得到瑞利分布[25]下的波高概率分布，并將所得結果與實測波高概率分布進行對比，對比結果如圖7所示，發現統計值和理論值的概率分布曲線擬合良好，認為波高的分布符合瑞利分布。

(13)

式中，f(H)為波高概率密度，P為對應的累計概率分布，Hmax為符合邊界要求的最大波高。

圖7 理論值和統計值對比Fig.7 Comparison of simulated values and statistical values

2.2 波浪參數對出水俯仰角的影響

波浪力可分為兩部分[16, 22]：①由波浪壓力場引起的壓差力；②由于波浪質點的軌跡速度引起的附加慣性力。波浪參數通過影響波浪力的大小從而影響水下航行體的出水姿態參數。波高主要影響波浪力的幅值，對航行體出水姿態參數有著直接影響；波浪的周期T、初相位ε和浪向通過改變水下航行體出水相位來影響波浪質點的軌跡速度，改變附加慣性力，從而影響出水姿態參數。

將不同工況下的尾部觸水俯仰角和靜水下尾部觸水俯仰角作差比得到俯仰角偏差值β：

(14)

圖8 不同波高下的俯仰角偏差值Fig.8 Deviation of pitch angle under different wave height

圖9 不同出水相位下的俯仰角偏差值Fig.9 Deviation of pitch angle under different water-exit phase

通過圖9可以發現，在該計算工況下，俯仰角偏差隨著出水相位呈余弦變化規律，因為隨著出水相位的變化，流體質點的軌跡速度是不同的，流體質點和波浪的傳播存在相對運動，受到波浪相位的影響，流體質點運動速度的方向近似呈余弦變化。浪流的相對運動如圖10所示，在波峰位置，波浪傳播方向和流體質點的運動方向相同；在波谷位置，波浪傳播方向和流體質點的運動方向相反。波峰位置和波谷位置浪流相對運動的不同，造成兩個相位下航行體尾部觸水俯仰角相反的情況：波峰位置出水，流體質點始終沿著x軸正方向運動[15]，此時彈體所受的附加慣性力偏向于航行體俯仰方向，波浪對尾部觸水俯仰角的負影響達到最大；波谷位置出水，流體質點會沿x軸負方向運動[15]，此時彈體所受的附加慣性力對航行體俯仰運動呈抑制作用，波浪對尾部觸水俯仰角的正影響達到最大。

圖10 浪流相對運動Fig.10 Relative motion of waves and fluid particles

2.3 水下航行體出水成功概率預報

在潛艇實際海上作戰中，海況是復雜多變的，在進行潛射導彈發射之前，需要對該海況下的發射成功概率進行快速預報。本小節針對此工程背景提出一種水下航行體出水成功概率的快速預報方法。

1)規則波浪數據的建立：為了提供基本的數據參考，將各個工況所算得的尾部觸水俯仰角進行集成，形成規則波浪數據庫。通過2.2節分析可知，相同波高下，俯仰角的變化受到出水相位的影響，在某一出水相位會出現最惡劣情況(即尾部觸水俯仰角最小)。所以，在生成規則波浪數據庫時，計算水下航行體在不同波高下最惡劣的工況即尾部觸水時的最小俯仰角納入數據庫中。波高的計算范圍參考北半球大洋海浪年平均統計資料，浪向取為逆浪，波高取變化范圍[0,14]m，波高間隔0.5 m，每個工況計算所得的尾部觸水俯仰角納入規則波浪數據庫中。通過圖8可以發現隨著波高的增加，俯仰角偏差值近似呈線性增加，所以，對于未曾計算到的波高(比如H=0.3 m)，可采用分段線性插值得到出水俯仰角。該數據庫用于后續水下航行體出水成功的判定。

3)概率預報模型的建立：波高H的分布服從瑞利分布[25]，所以參考式(13)建立波高分布概率預報模型。

該實測海況下概率分布如圖11所示，Hmax為滿足要求的最大波高，0～Hmax為所滿足的波浪范圍，通過概率分布函數即可得到水下航行體的出水成功概率。水下航行體出水成功概率預報的流程如圖12所示。

圖11 波高的概率分布Fig.11 Probability distribution of wave height

圖12 水下航行體出水成功概率預報流程Fig.12 Successful probability prediction process of underwater vehicles coming out of water

表3 水下航行體出水成功概率預報結果

3 結論

本文通過邊界元計算方法實現了對水下航行體水下發射的數值模擬，對水下航行體在規則波浪下的出水姿態參數形成了規模計算，從而建立了規則波浪數據庫。計算結果表明，在考慮最惡劣工況下，出水俯仰角偏差隨著波高的增加而減小；相同波高下，出水俯仰角偏差隨著出水相位呈余弦變化，俯仰角極值發生在波峰和波谷位置。

以尾部觸水時的俯仰角作為水下航行體能否成功出水的判定，在認為波浪的波高符合瑞利分布的前提下，實現了在實測海浪下對水下航行體出水成功概率的預報。基于規則波浪數據庫下的水下航行體出水成功概率預報方法具有如下特性：

1)快速性：該預報方法是從已經計算完成的數據庫中抓取邊界值，利用波浪分布的概率模型對出水成功概率進行計算，用時短暫，對于潛艇作戰判定發射時機具有一定參考意義。

2)靈活性：該方法是參考北半球大洋海浪年平均統計資料進行工況劃分，在實際應用時，可以參考不同海域的統計資料進行計算以形成數據庫，進行預報的彈體模型也可以根據實際進行更換。

3)實用性：由于計算結果的線性變化，該方法只需要一定數量的原始數據，數據獲取容易而且精度較高，具有一定的實用性。