內(nèi)波作用下水下航行體應急上浮的運動及水動力數(shù)值模擬研究

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

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

0 引 言

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

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

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

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

1 數(shù)值計算方法

1.1 控制方程

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

連續(xù)性方程：

動量方程：

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

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

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

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

1.2 數(shù)值造波理論

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

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

1.3 計算模型及網(wǎng)格劃分

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

表1 SUBOFF標模主要參數(shù)Tab.1 Main parameters of SUBOFF model

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

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

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

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

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

圖3 SUBOFF模型附近的網(wǎng)格劃分Fig.3 Meshes around SUBOFF model

2 數(shù)值計算方法驗證

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

水下航行體的潛深變化實時監(jiān)測對比如圖4所示，其中實線為水池試驗結果，虛線為數(shù)值計算結果。拋載之后航行體的潛深在初期變化較小，后期隨著速度的增加，潛深快速減小，直至接近水面。這符合物體從靜止狀態(tài)在外力作用下加速度運動的物理現(xiàn)象。圖5是水下航行體上浮過程中縱傾角變化的對比曲線，上浮過程中抬首角度先是增大，之后變小，并快速轉變?yōu)橐欢ǔ潭鹊穆袷住Ρ瓤梢?，應急上浮運動的數(shù)值模擬結果與試驗結果一致，也說明本文建立的水下航行體應急上浮運動數(shù)值計算方法是可靠的。

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

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

3 數(shù)值結果分析

3.1 自由運動

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

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

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

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

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

圖8 是水下航行體運動過程中自身俯仰角度的變化曲線。圖9 是運動過程中水下航行體受到內(nèi)波水動力作用產(chǎn)生的俯仰力矩時歷曲線，力矩作用中心為航行體質心位置，其中無量綱俯仰力矩系數(shù)為，My為俯仰力矩。My正值代表抬首力矩，負值為埋首力矩。整體來看，在185 s左右其遭受的埋首力矩明顯大于抬首力矩，所以這時俯仰角度快速增大（埋首）；之后俯仰力矩逐漸變?yōu)檎担?05 s附近抬首力矩占主導，因此航行體埋首角度快速減??；之后俯仰力矩又變?yōu)樨撝担栽?20 s開始又出現(xiàn)了埋首加劇現(xiàn)象，并在之后出現(xiàn)輕度反復?？梢娝潞叫畜w的姿態(tài)與俯仰力矩密切相連。水下航行體除了受內(nèi)波流場水動力的作用外，還受到自身恢復力矩的作用，在兩者同時作用下，水下航行體不斷運動從而尋求平衡，這是一個典型的水動力與運動相互耦合的作用過程。

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

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

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

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

3.2 應急上浮

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

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

圖12～14 中虛線對應的時間t=225 s，即開始實施應急上浮的時刻。圖12 是實施應急上浮措施前后的水下航行體位移變化曲線。可以看出，在實施應急上浮之后，由于慣性作用，水下航行體繼續(xù)下沉了一段距離，但很快垂向位移與水平位移曲線都出現(xiàn)了拐點變化，垂直方向由下沉轉變?yōu)榭焖偕细。椒较虮憩F(xiàn)為水下航行體的前移速度趨緩?？傮w來看，實施應急上浮措施之后，水下航行體朝著偏離內(nèi)波波谷的方向發(fā)展。

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

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

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

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

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

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

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

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

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

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

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

4 結 論

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

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

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