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

（中國(guó)船舶科學(xué)研究中心水動(dòng)力學(xué)重點(diǎn)實(shí)驗(yàn)室，江蘇 無(wú)錫 214082）

0 引 言

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

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

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

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

1 數(shù)值計(jì)算方法

1.1 控制方程

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

連續(xù)性方程：

動(dòng)量方程：

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

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

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

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

1.2 數(shù)值造波理論

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

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

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

本文研究對(duì)象為國(guó)際通用的SUBOFF標(biāo)模，主尺度見(jiàn)表1。

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

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

圖1 坐標(biāo)系的定義Fig.1 Definition of the coordinate systems

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

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

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

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

2 數(shù)值計(jì)算方法驗(yàn)證

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

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

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

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

3 數(shù)值結(jié)果分析

3.1 自由運(yùn)動(dòng)

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

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

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

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

圖7 初始運(yùn)動(dòng)時(shí)刻的速度場(chǎng)云圖Fig.7 Velocity distributions at start moving moment

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

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

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

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

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

3.2 應(yīng)急上浮

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

圖11 上浮時(shí)刻水下航行體周?chē)牧鲌?chǎng)Fig.11 Flow field around submarine at the moment to float

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

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

圖13 是應(yīng)急上浮過(guò)程的水下航行體俯仰角度的時(shí)歷曲線?？梢钥吹阶杂蛇\(yùn)動(dòng)階段一直保持埋首姿態(tài)下沉的水下航行體，在實(shí)施應(yīng)急上浮之后，俯仰角快速轉(zhuǎn)為抬首，在t=238 s 時(shí)抬首角度達(dá)到最大，之后變小并轉(zhuǎn)為埋首。

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

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

圖14 應(yīng)急上浮垂向運(yùn)動(dòng)速度時(shí)歷曲線Fig.14 Time history of vertical velocity

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

圖15 應(yīng)急上浮過(guò)程中的速度場(chǎng)云圖Fig.15 Velocity distributions around floating submarine

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

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

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

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

4 結(jié) 論

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

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

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