內波作用下水下航行體應急上浮的運動及水動力數值模擬研究

黃苗苗，張 楠，李迎華，王文濤

（中國船舶科學研究中心水動力學重點實驗室，江蘇 無錫 214082）

0 引 言

海洋內波是發生在密度穩定層化的海水內部的一種波動，是世界各大洋尤其是邊緣海中一種普遍的海洋現象，其最大振幅出現在海水內部[1]。相對于海水與空氣接觸的自由面上下密度差來說，內波躍層的密度差很小，相當于將分層介質置于微重力場中，所以恢復力也更小，僅為水面波的0.1%量級[2]。因此在能量相同的條件下內波波幅可以是水面波的20～30倍。海洋內波特別是強流速、大波幅的內孤立波對海洋結構物具有極大的破壞性。

在海洋內波對水下航行體的水動力影響方面，國內開展了一定研究[3-7]。試驗及計算結果均表明，遭遇海洋內波會引起水動力突變，由此影響水下航行體的安全航行。對于水下航行體來講，在危急情況比如舵卡、損失浮力和危險縱傾時，可采用高壓氣吹除主壓載水艙的方式實現應急上浮，這也是水下航行體抗沉的重要措施[8]。水下航行體遭遇內波，必要時亦可采用這種方式脫險，但是這方面的研究較少。

目前，水下航行體運動模型形式較多，其中美國的葛特勒方程最具權威性。由于水動力系數項較多且求解困難，目前研究應急上浮運動大多數使用應急上浮簡化運動模型，其中使用較多的是垂直面運動方程，其他附加力的計算模型主要采用經驗公式來求解[9-10]。隨著近幾年CFD 技術及計算機硬件的快速發展，越來越多的學者開始采用數值計算方法來研究水動力與運動問題。美國的Carrica等[11-12]基于自主研發的流體軟件CFDShip-Iowa，實現了不同工況下航行體六自由度運動的數值模擬計算，并對其操縱性進行了預報；Bettle[13-14]等和Watt等[15]用六自由度求解器數值模擬并分析了水下航行體上浮過程的運動特性；國內周廣禮等[16]和閆朋[17]采用粘流方法開展了靜水中的應急上浮運動研究，為解決復雜環境下的應急上浮問題打下了良好基礎。

在前期開展的內波對水下航行體水動力、運動特征的影響研究[6，18]的基礎上，本文考慮海洋內波的非定常影響，開展內波環境下的應急上浮運動數值模擬研究，針對無航速水下航行體模擬實施高壓氣吹除主壓載水艙的上浮運動，并對上浮過程中的運動特征開展詳細分析。本研究可為后續內波下有航速的水下航行體主動操控研究打下基礎。

1 數值計算方法

1.1 控制方程

本文數值計算采用的控制方程包括連續性方程及動量方程，湍流模式選取Realizablek-ε模型。

連續性方程：

動量方程：

式中，ρ為密度，P為壓力，τxx、τxy和τxz等是因分子粘性而產生的作用在微元體表面的粘性應力τ的分量，Fx、Fy和Fz是微元體上的體積力分量。

依據牛頓第二定律，應用質心運動定理和相對于質心的動量矩定理，水下航行體的六自由度運動方程表達為

式中，F與M分別為航行體所受外力及外力矩，B=（mu,mv,mw）為動量，K=（Ixp,Iyq,Izr）動量矩，U=（u,v,w）為速度，Ω=（p,q,r）為角速度。

內波對于水下航行體運動的影響主要表現在垂直面上。因此，本文數值模擬中采用了三自由度運動控制方程，具體為

1.2 數值造波理論

本文采用造波邊界條件進行內孤立波的數值造波，通過波的解析解或數值解在邊界上給定一個波動速度作為邊界條件，以實現造波。根據內孤立波的定義可知，海洋內波產生的根源是海水密度分層，因此不管是物理試驗還是數值模擬，造波必然是在分層流中進行的。根據海洋學中著名的“剛蓋假設”[2]，將海水表面視為剛性平面，內波簡化為兩層流問題。本文根據KdV 理論[19]開展數值造波，內孤立波的波面表達式為

式中，a為波幅，λ為內孤立波的特征波長，c為內孤立波的傳播速度，x為坐標值，t為時間。

1.3 計算模型及網格劃分

本文研究對象為國際通用的SUBOFF標模，主尺度見表1。

表1 SUBOFF標模主要參數Tab.1 Main parameters of SUBOFF model

開展內波對于水下航行體影響研究的前提是建立內孤立波的數值造波水池。如圖1 所示，計算選用的數值造波水池總長為1000 m，上層水深為40 m，下層水深為175 m，上下兩層流體的密度比ρ1/ρ2=0.993，波幅a=5.0 m。計算中大地坐標系為E-XYZ，X軸正方向指向內孤立波前進方向，Z軸正方向垂直指向上方，Y軸正方向滿足右手法則。另外，為了在數值模擬中計算水下航行體的自身運動，建立隨水下航行體運動的動坐標系G-xyz。該坐標系滿足右手法則，原點取在水下航行體的質心處，x軸正方向指向航行體首部，y軸正方向指向右舷。

圖1 坐標系的定義Fig.1 Definition of the coordinate systems

采用重疊網格方法處理運動計算過程中的網格問題。重疊網格方法將復雜的流動區域分成幾何邊界比較簡單的子區域，流場信息通過插值在重疊區邊界進行匹配和耦合。本文計算域包括兩套網格：背景網格和水下航行體網格。背景網格計算域的邊界條件設置如圖2 所示，具體包括：速度入口造波；壓力出口；水池頂部根據海水自由表面剛性假定設為壁面；水池底部也設置為壁面。

圖2 計算域和邊界條件Fig.2 Computational model and boundary conditions

在計算域中生成六面體網格，對航行體模型表面附近網格加密，其中第一層網格間距根據y+確定（y+范圍為40～100）。為更好地捕捉波面，網格劃分時對內孤立波通過的區域進行了網格加密。此外，為準確地模擬水下航行體的水動力及流場特征，對模型主附體周圍的網格進行細化處理，艇體表面及近壁區域的網格劃分如圖3 所示。本文計算模型網格總數為805 萬，其中數值水池背景網格為405萬，水下航行體區域網格為400萬。

圖3 SUBOFF模型附近的網格劃分Fig.3 Meshes around SUBOFF model

2 數值計算方法驗證

為了驗證應急上浮運動的數值計算方法，本文開展了水下航行體應急上浮運動數值模擬，并與相同工況的水池模型試驗對比。該試驗在中國船舶科學研究中心的露天水池開展，通過調整模型的配重，使模型處于靜均衡狀態。初始時刻水下航行體位于水深7 m 的位置，初始航速為0。試驗狀態為靜水中拋載自由上浮，拋載量為排水量的1.48%。

水下航行體的潛深變化實時監測對比如圖4所示，其中實線為水池試驗結果，虛線為數值計算結果。拋載之后航行體的潛深在初期變化較小，后期隨著速度的增加，潛深快速減小，直至接近水面。這符合物體從靜止狀態在外力作用下加速度運動的物理現象。圖5是水下航行體上浮過程中縱傾角變化的對比曲線，上浮過程中抬首角度先是增大，之后變小，并快速轉變為一定程度的埋首。對比可見，應急上浮運動的數值模擬結果與試驗結果一致，也說明本文建立的水下航行體應急上浮運動數值計算方法是可靠的。

圖4 潛深變化對比曲線Fig.4 Comparison of depth change curves

圖5 俯仰角變化對比曲線Fig.5 Comparison of pitch angle history curves

3 數值結果分析

3.1 自由運動

水下航行體初始位置位于分層流的下方水體，相對于波幅0.4 倍深度的位置，初始時刻水下航行體自身重力與浮力相等，航行體水平放置首部迎流、0 航速。在波形穩定了之后水下航行體開始自由運動。

內波作用下，水下航行體的自由運動模擬結果如下述所示。圖6 是航行體運動過程中水平位移（Δx）及垂向位移（Δz）的時歷曲線。從中可以看出，遭遇內孤立波后，水下航行體迅速下沉并被拖向內波波谷方向，自由運動75 s 之后垂向位移已接近一倍航行體長度。說明水下航行體遭遇大幅內波，若不采取主動操控，潛深會大幅度變化，并很可能到達危險深度。

圖6 模型位移變化時歷曲線Fig.6 Displacement histories of the model

初始運動時刻水下航行體周圍的速度場如圖7 所示。整體來看，航行體周圍流場的水平速度指向左方，航行體垂直速度指向下方。在這種速度場的作用下，水下航行體自由運動時自然朝左下方移動，也就是下沉并被卷拖進內波波谷中心。這里需要說明的是：在圖7（b）局部放大圖中，水下航行體近體下方的垂向速度場與其他區域不同，這主要是由于航行體受到來自上方的強烈速度場作用，在其近體下方形成了典型的背風區特征。可以說內波的出現改變了水下航行體周圍的速度場，速度場的改變引起了水下航行體表面壓力的變化，從而導致其運動的改變。

圖7 初始運動時刻的速度場云圖Fig.7 Velocity distributions at start moving moment

圖8 是水下航行體運動過程中自身俯仰角度的變化曲線。圖9 是運動過程中水下航行體受到內波水動力作用產生的俯仰力矩時歷曲線，力矩作用中心為航行體質心位置，其中無量綱俯仰力矩系數為，My為俯仰力矩。My正值代表抬首力矩，負值為埋首力矩。整體來看，在185 s左右其遭受的埋首力矩明顯大于抬首力矩，所以這時俯仰角度快速增大（埋首）；之后俯仰力矩逐漸變為正值，205 s附近抬首力矩占主導，因此航行體埋首角度快速減小；之后俯仰力矩又變為負值，所以在220 s開始又出現了埋首加劇現象，并在之后出現輕度反復。可見水下航行體的姿態與俯仰力矩密切相連。水下航行體除了受內波流場水動力的作用外，還受到自身恢復力矩的作用，在兩者同時作用下，水下航行體不斷運動從而尋求平衡，這是一個典型的水動力與運動相互耦合的作用過程。

圖8 模型俯仰角時歷曲線Fig.8 Pitch angle history of the model

圖9 俯仰力矩系數時歷曲線Fig.9 History curve of pitch moment coefficient

另外值得注意的是，圖9 中水下航行體在運動過程中俯仰力矩的時歷振蕩十分顯著。為了研究其振蕩特性，對該時域數據開展傅立葉變換，獲得功率譜密度（power spectral density，PSD）分布，結果如圖10 所示。功率譜密度較大的幾個點對應的頻率分別是0.034 Hz、0.026 Hz 和1.05 Hz。這說明水下航行體在遭遇內波后屬于低頻運動變化，是一個非穩態的過程，圖9力矩曲線的振蕩特征是正常的。

圖10 俯仰力矩系數的功率譜密度Fig.10 Power spectral density of pitch moment coefficient

3.2 應急上浮

由上一節的模擬結果可知，在內波速度場的卷帶作用下，水下航行體向著內波波谷中心快速下沉。如圖11 所示，圖中曲線為內波的波面，此時水下航行體垂向位移Δz/L＞0.8。針對該情況，實施水下航行體的應急上浮。數值計算中采用施加可變外力的方式來模擬高壓氣吹除主壓載水艙，吹除時間按照斯特羅哈爾數相似原則計算給定。

圖11 上浮時刻水下航行體周圍的流場Fig.11 Flow field around submarine at the moment to float

圖12～14 中虛線對應的時間t=225 s，即開始實施應急上浮的時刻。圖12 是實施應急上浮措施前后的水下航行體位移變化曲線。可以看出，在實施應急上浮之后，由于慣性作用，水下航行體繼續下沉了一段距離，但很快垂向位移與水平位移曲線都出現了拐點變化，垂直方向由下沉轉變為快速上浮，水平方向表現為水下航行體的前移速度趨緩。總體來看，實施應急上浮措施之后，水下航行體朝著偏離內波波谷的方向發展。

圖12 模型應急上浮垂向及水平位移時歷曲線Fig.12 Vertical and horizontal displacement curves of floating submarine

圖13 是應急上浮過程的水下航行體俯仰角度的時歷曲線。可以看到自由運動階段一直保持埋首姿態下沉的水下航行體，在實施應急上浮之后，俯仰角快速轉為抬首，在t=238 s 時抬首角度達到最大，之后變小并轉為埋首。

圖13 應急上浮俯仰角變化曲線Fig.13 Pitch angle history of SUBOFF emergency floating

圖14 是航行體垂向速度監測曲線，下沉為負，上升為正。從中可以看到，在實施應急上浮措施之后，垂向速度由負值逐漸拉升為正值。值得注意的是，垂向速度在t=238 s時達到最大值，之后逐漸變小。

圖14 應急上浮垂向運動速度時歷曲線Fig.14 Time history of vertical velocity

圖15（a）～（d）是上浮過程中水下航行體周圍的速度場及內波的波形變化。圖15（a）是實施應急上浮的初始時刻，水下航行體周圍速度場主要受內波的影響，波面上方速度較大，波面下方速度較小；從圖15（b）～（c）可以明顯看出在實施應急上浮操縱之后，水下航行體開始抬首、上浮，航行體下方附近的海水由于粘帶作用速度也越來越大；圖15（d）表示水下航行體穿越波面之后，抬首角度變小，之后轉為埋首。

圖15 應急上浮過程中的速度場云圖Fig.15 Velocity distributions around floating submarine

圖16 是航行體上浮過程中的波形變化云圖，波面上方是低密度海水，波面下方是高密度海水為平均密度，ρ'=。從中可以看到航行體與內波波面的變化，在這個過程中水下航行體逐步上浮、接近內波波面并最終穿越波面。總的來講，水下航行體上浮穿越內波波面的運動特征與水下航行體的應急上浮出水現象相似。

圖16 應急上浮過程中的水下航行體和波形變化Fig.16 Model and internal waves in the process of floating up

為了探究水下航行體在內波中應急上浮運動特征的形成原因，本文進一步開展了理論分析。圖17是航行體穿越內波波面前后的受力分析示意圖，其中G為重力、F為浮力。當航行體以抬首姿態沖出波面時，在波面上下海水密度差的影響下，浮力變小并且作用點后移，而重力的大小及位置不變，于是形成了一對較強的逆時針作用力矩。在這對力矩的作用下，水下航行體沖出內波波面之后抬首角度迅速變小并開始轉向埋首姿態。這就是前面數值模擬中航行體俯仰角及上浮速度時歷曲線在238 s時出現拐點變化的原因。

圖17 受力分析示意圖Fig.17 Forces diagram of submarine

4 結 論

本文采用RANS 方法結合運動控制方程，開展了海洋內波作用下的水下航行體運動特征數值計算研究，其中靜水中的應急上浮運動與模型試驗吻合，證明了計算方法切實可行，主要結論如下：

（1）數值計算表明，處于分層流下方的水下航行體在遭遇內波之后，即使未穿越波面，也會由于速度場的改變導致航行體表面壓力的改變，從而出現埋首下沉現象。這個過程中，水下航行體除了受內波流場水動力的作用外，還受到自身恢復力矩的作用，在兩者同時作用下水下航行體不斷運動從而尋求平衡，這是一個典型的水動力與運動相互耦合的作用過程。

（2）在此研究基礎上，文中采用數值計算方法模擬了高壓氣吹除主壓載水艙的作用，針對水下航行體遭遇內波迅速下沉的情況，開展了應急上浮運動的數值模擬。結果表明，實施緊急操控之后，水下航行體能夠快速改變運動軌跡、由下沉轉變為上浮并脫離內波影響；上浮后期，航行體以抬首姿態快速沖出波面之后轉為埋首。文中對上浮運動特征開展了詳細的分析，數值模擬結果與理論分析一致。