基于體積力模型的潛艇應(yīng)急上浮運(yùn)動(dòng)數(shù)值模擬分析

魏可可，高霄鵬

海軍工程大學(xué) 艦船與海洋學(xué)院，湖北 武漢 430033

0 引 言

潛艇的應(yīng)急上浮是指當(dāng)發(fā)生艙室破損進(jìn)水、卡舵等事故時(shí)，采取排除部分或者全部主壓載水艙內(nèi)的水，從而使?jié)撏а杆俑〕鏊娴囊环N緊急操縱措施。為保證潛艇應(yīng)急上浮的安全性和穩(wěn)定性，準(zhǔn)確預(yù)報(bào)潛艇在應(yīng)急上浮過程中的運(yùn)動(dòng)特性十分有必要。

和正常的上浮相比，潛艇應(yīng)急上浮中艇的排水方式不同且排水速率快，艇體姿態(tài)變化幅度大。常見的潛艇應(yīng)急上浮運(yùn)動(dòng)預(yù)報(bào)方法是通過拘束模試驗(yàn)[1-4]或參數(shù)辨識(shí)[5]等方法求解水動(dòng)力系數(shù)，然后再基于六自由度運(yùn)動(dòng)模型進(jìn)行預(yù)報(bào)。戴余良等[6-7]應(yīng)用非線性系統(tǒng)運(yùn)動(dòng)穩(wěn)定性與分叉理論，以及同倫延拓?cái)?shù)值計(jì)算方法，分析了潛艇應(yīng)急上浮運(yùn)動(dòng)的穩(wěn)態(tài)響應(yīng)及運(yùn)動(dòng)穩(wěn)定性。王曉玢等[8]通過對(duì)潛艇垂直面操縱運(yùn)動(dòng)進(jìn)行非線性建模，利用分叉與突變理論方法對(duì)潛艇上浮運(yùn)動(dòng)中由靜態(tài)分叉和動(dòng)態(tài)分叉引發(fā)的狀態(tài)突變進(jìn)行了分析。潛艇水動(dòng)力系數(shù)的求解較復(fù)雜且有一定的難度，而上述預(yù)報(bào)方法主要是針對(duì)潛艇某自由度上的運(yùn)動(dòng)，并不能真實(shí)反映潛艇在六自由度完全釋放情況下的應(yīng)急上浮運(yùn)動(dòng)特性，且預(yù)報(bào)精度還有待提高。

隨著計(jì)算機(jī)技術(shù)的快速發(fā)展，計(jì)算流體力學(xué)在船舶運(yùn)動(dòng)研究中的應(yīng)用越來越廣泛。Carrica等[9-10]基于自主研發(fā)的流體軟件CFDShip-Iowa V4，實(shí)現(xiàn)了不同工況下潛艇六自由度運(yùn)動(dòng)的數(shù)值模擬計(jì)算，并對(duì)其操縱性進(jìn)行了較高精度的預(yù)報(bào)。Chase 等[11]采用重疊網(wǎng)格技術(shù)，對(duì)Suboff 潛艇模型的直航拖曳和自由自航等運(yùn)動(dòng)進(jìn)行了數(shù)值模擬計(jì)算。周廣禮等[12-13]基于RANS 方程及流體體積模型，采用整體動(dòng)網(wǎng)格技術(shù)，對(duì)Suboff 模型的應(yīng)急上浮六自由度運(yùn)動(dòng)及其出水運(yùn)動(dòng)予以了研究。廖歡歡等[14]采用流體體積法（volume of fluid，VOF）六自由度求解器和重疊網(wǎng)格技術(shù)，研究了潛艇在一定潛深下的應(yīng)急上浮六自由度運(yùn)動(dòng)，并對(duì)上浮過程中的姿態(tài)和水動(dòng)力進(jìn)行了分析。

以上大多學(xué)者對(duì)潛艇應(yīng)急上浮運(yùn)動(dòng)的研究都沒有考慮螺旋槳以及初始航速等問題的嵌入，為此，本文將基于STAR-CCM+軟件平臺(tái)，采用基于體積力模型的虛擬槳盤，針對(duì)潛艇在不同航速下的應(yīng)急上浮運(yùn)動(dòng)進(jìn)行數(shù)值模擬，并對(duì)比分析潛艇在應(yīng)急上浮運(yùn)動(dòng)中的姿態(tài)和上浮運(yùn)動(dòng)時(shí)間等，用以為潛艇在實(shí)際中的上浮運(yùn)動(dòng)操控提供一定的借鑒。

1 計(jì)算對(duì)象

本文將以美國泰勒研究中心提供的Suboff 潛艇試驗(yàn)?zāi)Ｐ蜑檠芯繉?duì)象。國內(nèi)外針對(duì)該模型開展了大量的拖曳及流場測量試驗(yàn)，其三維模型圖如圖1 所示，模型主尺度參數(shù)如表1 所示。

圖1 Suboff 潛艇三維模型Fig. 1 Three dimensional model of Suboff submarine

表1 模型主要參數(shù)Table 1 Main parameters of the model

2 數(shù)值模擬

2.1 控制方程及湍流模型

應(yīng)用黏流CFD 方法開展數(shù)值計(jì)算的基礎(chǔ)是納維?斯托克斯（Navia-Stokes，N-S）方程，然而受限于計(jì)算能力，當(dāng)前還不能通過直接求解N-S 方程（direct numerical simulation，DNS）來獲得潛艇周圍復(fù)雜的黏性流場。將N-S 方程中的湍流脈動(dòng)項(xiàng)進(jìn)行時(shí)均化處理可以極大地減小計(jì)算量，所得雷諾平均方程（RANS）的具體表達(dá)如式（1）所示。然而，N-S 方程在時(shí)均化的過程中增加了含未知量時(shí)均的雷諾應(yīng)力項(xiàng)，通過引入湍流模型，可使方程封閉。目前，該方法已在艦艇流體力學(xué)理論研究及工程實(shí)踐中得到廣泛應(yīng)用。

開展計(jì)算時(shí)，控制方程離散均采用二階迎風(fēng)格式，并應(yīng)用分離式解法（segregated flow）對(duì)離散后的方程組進(jìn)行求解。壓力速度耦合迭代采用SIMPLE 算法，涉及非穩(wěn)態(tài)計(jì)算時(shí)，設(shè)定計(jì)算過程中的時(shí)間步長為 ?t且固定不變，時(shí)間迭代采用二階形式，即

采用RANS 方法開展?jié)撏畡?dòng)力計(jì)算時(shí)，湍流模型的選取是確保計(jì)算精度的關(guān)鍵要素，眾多學(xué)者已針對(duì)該問題開展了大量研究，但所得結(jié)論存在較大的不同，且湍流模型的適用性也與網(wǎng)格布局及操作者的數(shù)值計(jì)算經(jīng)驗(yàn)相關(guān)。綜合各類文獻(xiàn)資料，本文選取Wilcoxk?ω 與SSTk?ω湍流模型進(jìn)行優(yōu)選分析。

Wilcoxk?ω湍 流 模型 最 早由Wilcox 于1988年提出，后經(jīng)過不斷的修正，模擬精度及其適用范圍均得到了提升，STAR-CCM+軟件平臺(tái)中就提供有最新的修正模型，具體表達(dá)形式如下：

式中：k為湍動(dòng)能； ω為特殊湍動(dòng)能耗散；其余參數(shù)的具體定義及參數(shù)值參見文獻(xiàn)[15]。

Menter 等[16]結(jié) 合 標(biāo) 準(zhǔn)k?ε 與k?ω湍 流 模 型，提出了剪切應(yīng)力模型SSTk?ω，式（4）給出了該模型的最新修正形式，具體推導(dǎo)過程及參數(shù)定義詳見文獻(xiàn)[15]。

2.2 計(jì)算域及網(wǎng)格劃分

計(jì)算域選取為方形域，除域尾端邊界設(shè)置為壓力出口外，其余邊界條件均設(shè)置為速度入口，域大小及計(jì)算域與自由液面的初始相對(duì)位置如圖2 所示。

進(jìn)行網(wǎng)格劃分時(shí)，為能更好地捕捉兩相流自由液面，對(duì)自由液面可能通過的區(qū)域進(jìn)行了網(wǎng)格加密。此外，為準(zhǔn)確模擬潛艇出水后艇體周圍的流場，對(duì)模型首部及指揮室圍殼周圍網(wǎng)格進(jìn)行了細(xì)化處理。同時(shí)，為滿足網(wǎng)格貼體的要求，對(duì)近壁面區(qū)域進(jìn)行了局部加密，以確保艇體壁面y+在30～300 之間[12]。當(dāng)y+在30～200 之間時(shí)，表征第1 層網(wǎng)格布置于對(duì)數(shù)層，可通過假定的壁面函數(shù)對(duì)近壁面流動(dòng)進(jìn)行模擬。該方法不直接求解黏性底層與過渡層內(nèi)的流動(dòng)，適用性較強(qiáng)，可極大地減少近壁面網(wǎng)格數(shù)量。艇體表面、域剖面及近壁面域的網(wǎng)格如圖3 所示，網(wǎng)格總數(shù)為700 萬。

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

圖3 艇體整體網(wǎng)格及域剖面網(wǎng)格Fig. 3 Overall grid and domain profile grid of hull

2.3 體積力模型

通過對(duì)螺旋槳入流面流速進(jìn)行加權(quán)積分，獲得槳的進(jìn)流速度，然后代入到已知螺旋槳水動(dòng)力曲線中，預(yù)估得到推力及扭矩，接著，通過添加源項(xiàng)的方式，對(duì)設(shè)定的體積力域內(nèi)的流體施加等效的力及力矩，該方法即為體積力模型。體積力模型不僅能夠?qū)⒘傲厥┘佑跐撏?，而且能更為?zhǔn)確地模擬由螺旋槳引起的潛艇尾流場變化，相較于直接力模型，其計(jì)算量大，但計(jì)算效率遠(yuǎn)優(yōu)于開展的真實(shí)螺旋槳計(jì)算。

為了重點(diǎn)分析潛艇的運(yùn)動(dòng)特征，提高計(jì)算效率，螺旋槳用基于體積力模型的虛擬槳代替。該虛擬槳的參數(shù)采用DTMB 4 383 標(biāo)準(zhǔn)槳的，模型主參數(shù)如表2 所示。

表2 DTMB 4 383 螺旋槳主要參數(shù)Table 2 Main parameters of DTMB 4 383 propeller

2.4 數(shù)值計(jì)算的驗(yàn)證

針對(duì)漂角 β=2?～14?、 間隔 2?的Suboff 模型的無因次縱向水動(dòng)力X′進(jìn)行數(shù)值計(jì)算，計(jì)算結(jié)果與試驗(yàn)值[17]的對(duì)比如圖4 所示。

圖4 不同漂角下艇體縱向水動(dòng)力的試驗(yàn)值與數(shù)值計(jì)算值的對(duì)比Fig. 4 Comparison between experimental and numerical values of longitudinal hydrodynamics of the hull at different drift angles

由圖4 可看出，不同漂角下艇體縱向水動(dòng)力數(shù)值計(jì)算結(jié)果與試驗(yàn)值吻合良好，證明該方法可行。因此，選取該方法進(jìn)行后續(xù)的數(shù)值仿真計(jì)算。

2.5 計(jì)算工況

為探討不同航速對(duì)潛艇上浮運(yùn)動(dòng)的影響規(guī)律，選取上浮力為總排水量的5%、螺旋槳轉(zhuǎn)速分別為300，400，500，600 r/min 的工況進(jìn)行數(shù)值模擬計(jì)算，具體計(jì)算工況如表3 所示。

表3 計(jì)算工況Table 3 Calculation conditions

3 不同航速下數(shù)值模擬結(jié)果分析

3.1 螺旋槳轉(zhuǎn)速與航速的匹配

假定潛艇上浮之前在水下做定深直航運(yùn)動(dòng)，為獲取模型水下穩(wěn)定直航狀態(tài)，設(shè)定初始直航航速U=0，在螺旋槳推力的作用下，潛艇沿x方向做直航運(yùn)動(dòng)。為避免在初始計(jì)算階段出現(xiàn)計(jì)算結(jié)果不穩(wěn)定的現(xiàn)象，在t=0～4 s 時(shí)間段內(nèi)將艇體按靜態(tài)處理，待螺旋槳推力及流場計(jì)算結(jié)果穩(wěn)定后，即在t1時(shí)刻，將x方向運(yùn)動(dòng)完全釋放，潛艇開始做直航運(yùn)動(dòng)。

不同轉(zhuǎn)速下螺旋槳推力及艇體阻力的時(shí)歷曲線如圖5 所示。由圖可知，艇體阻力及螺旋槳推力隨時(shí)間逐漸向自航值逼近，且二者隨時(shí)間的絕對(duì)變化率較為相近，航速隨時(shí)間逐步增大并在t2時(shí)刻后趨于穩(wěn)定；隨著轉(zhuǎn)速的增大，螺旋槳獲得的推力以及艇體的阻力值越大，其運(yùn)動(dòng)穩(wěn)定的時(shí)間就越短，最終穩(wěn)定的航速值也越大，這從圖6 即可看出。

表4 給出了不同轉(zhuǎn)速下潛艇做定深直航運(yùn)動(dòng)時(shí)的計(jì)算結(jié)果。由表4 可知，4 個(gè)轉(zhuǎn)速下t1值間的差別較小，在t2時(shí)刻是隨轉(zhuǎn)速的增大而減小，轉(zhuǎn)速越大，需要穩(wěn)定的時(shí)間就越短；轉(zhuǎn)速為300～600 r/min 時(shí)，最終穩(wěn)定的定深下的航速分別為0.97，1.35，1.72 和2.1 m/s。

3.2 艇體上浮運(yùn)動(dòng)

圖7 給出了不同轉(zhuǎn)速下艇體垂向位移運(yùn)動(dòng)的時(shí)歷曲線。由圖可知，艇體做應(yīng)急上浮至出水一般會(huì)經(jīng)歷水下定深直航運(yùn)動(dòng)—上浮運(yùn)動(dòng)—水面航行這3 種狀態(tài)的運(yùn)動(dòng)。當(dāng)螺旋槳轉(zhuǎn)速為300 r/min時(shí)，其水下定深直航時(shí)間為101 s，上浮時(shí)間8 s，上浮出水的時(shí)刻為109 s；當(dāng)轉(zhuǎn)速為400 r/min 時(shí)，其水下定深直航時(shí)間為92 s，上浮時(shí)間9 s，上浮出水的時(shí)刻為101 s；當(dāng)轉(zhuǎn)速為500 r/min 時(shí)，其水下定深直航時(shí)間為91 s，上浮時(shí)間9 s，上浮出水的時(shí)刻為100 s；當(dāng)轉(zhuǎn)速為600 r/min 時(shí)，其水下定深直航時(shí)間為71 s，上浮時(shí)間12 s，上浮出水的時(shí)刻為83 s。由此可知，隨著轉(zhuǎn)速的增大，潛艇獲得穩(wěn)定航速的時(shí)間越短，艇體在水下做定深直航運(yùn)動(dòng)的時(shí)間也就越短；不同轉(zhuǎn)速下艇體上浮的時(shí)間雖然有所差別，但其差值量級(jí)較小，上浮運(yùn)動(dòng)的時(shí)間主要還是受上浮力的影響；同時(shí)，潛艇的上浮運(yùn)動(dòng)是一個(gè)非線性的復(fù)雜的六自由度運(yùn)動(dòng)，當(dāng)上浮力一致時(shí)，螺旋槳轉(zhuǎn)速越大，其上浮運(yùn)動(dòng)速度不一定越大。

圖5 不同轉(zhuǎn)速下艇體阻力和螺旋槳推力時(shí)歷曲線Fig. 5 Time histories of hull resistance and propeller thrust at different rotation speeds

圖6 不同轉(zhuǎn)速下艇體航速的時(shí)歷曲線Fig. 6 Time histories of hull velocity at different rotation speeds

表4 潛艇定深直航運(yùn)動(dòng)計(jì)算結(jié)果Table 4 Calculational results of direct sailing motion of submarine at fixed depth

圖7 不同轉(zhuǎn)速下艇體的垂向位移運(yùn)動(dòng)時(shí)歷曲線Fig. 7 Time histories of vertical displacement motion at different rotation speeds

3.3 艇上浮時(shí)的姿態(tài)分析

3.3.1 縱傾角分析

圖8 所示為不同轉(zhuǎn)速下艇體縱傾角的時(shí)歷曲線，表5 給出了縱傾角第1 次產(chǎn)生變化的時(shí)刻、發(fā)生峰值的時(shí)刻、峰值以及來回振蕩波動(dòng)的時(shí)間。由圖8 和表5 可知，當(dāng)轉(zhuǎn)速為300 r/min 時(shí)，縱傾角第1 次產(chǎn)生變化的時(shí)刻是100.2 s，第1 次產(chǎn)生峰值的時(shí)刻是109 s，其峰值為?22°，縱傾角來回振蕩波動(dòng)的時(shí)間是29.8 s；當(dāng)轉(zhuǎn)速為400 r/min時(shí)，縱傾角第1 次產(chǎn)生變化的時(shí)刻是92 s，第1 次產(chǎn)生峰值的時(shí)刻是99 s，其峰值為?23°，縱傾角來回振蕩波動(dòng)的時(shí)間是47 s；轉(zhuǎn)速為500 r/min 時(shí)，縱傾角第1 次產(chǎn)生變化的時(shí)刻是91 s，第1 次產(chǎn)生峰值的時(shí)刻是96 s，其峰值為?22°，縱傾角來回振蕩波動(dòng)的時(shí)間是47 s；當(dāng)轉(zhuǎn)速為600 r/min 時(shí)，縱傾角第1 次產(chǎn)生變化的時(shí)刻是71 s，第1 次產(chǎn)生峰值的時(shí)刻是75 s，其峰值為?27°，縱傾角來回振蕩波動(dòng)的時(shí)間是54 s。結(jié)果顯示，隨著轉(zhuǎn)速的增大，縱傾角第1 次產(chǎn)生變化及峰值的時(shí)刻呈減小趨勢，而縱傾角來回振蕩波動(dòng)的時(shí)間則越來越大，即轉(zhuǎn)速越大，艇體出水時(shí)受自由液面擾動(dòng)越小?？v傾角的峰值變化主要發(fā)生在潛艇上浮運(yùn)動(dòng)過程中，來回振蕩波動(dòng)主要發(fā)生在潛艇出水時(shí)刻。對(duì)比3.2 節(jié)有關(guān)上浮的時(shí)刻可知，無論轉(zhuǎn)速多大，潛艇第1 次傾角的峰值都是發(fā)生在潛艇上浮運(yùn)動(dòng)過程中，由對(duì)比分析可知，隨著轉(zhuǎn)速的增大，縱傾角的峰值并未呈現(xiàn)出明顯的變化規(guī)律。

圖8 不同轉(zhuǎn)速下艇體縱傾角時(shí)歷曲線Fig. 8 Time histories of trim angle at different rotation speeds

表5 不同轉(zhuǎn)速下縱傾角的分析Table 5 Analysis of trim angle at different rotation speeds

3.3.2 橫傾角分析

圖9 不同轉(zhuǎn)速下艇體橫傾角時(shí)歷曲線Fig. 9 Time histories of heeling angle at different rotation speeds

表6 不同轉(zhuǎn)速下橫傾角的分析Table 6 Analysis of heeling angle at different rotation speeds

圖9 所示為不同轉(zhuǎn)速下艇體橫傾角的時(shí)歷曲線，表6 給出了橫傾角第1 次產(chǎn)生變化的時(shí)刻、發(fā)生峰值時(shí)刻、峰值大小以及來回振蕩波動(dòng)的時(shí)間。由圖9 和表6 可知，當(dāng)轉(zhuǎn)速為300 r/min 時(shí)，橫傾角第1 次產(chǎn)生變化的時(shí)刻是101 s，第1 次產(chǎn)生峰值的時(shí)刻是108 s，其峰值為28.5°，橫傾角來回振蕩波動(dòng)的時(shí)間是87 s；當(dāng)轉(zhuǎn)速為400 r/min時(shí)，橫傾角第1 次產(chǎn)生變化的時(shí)刻是92 s，第1 次產(chǎn)生峰值的時(shí)刻是99 s，其峰值為22.3°，橫傾角來回振蕩波動(dòng)的時(shí)間是79 s；當(dāng)轉(zhuǎn)速為500 r/min時(shí)，橫傾角第1 次產(chǎn)生變化的時(shí)刻是88 s，第1 次產(chǎn)生峰值的時(shí)刻是98 s，其峰值為18.1°，橫傾角來回振蕩波動(dòng)的時(shí)間是73 s；當(dāng)轉(zhuǎn)速為600 r/min時(shí)，橫傾角第1 次產(chǎn)生變化的時(shí)刻是70 s，第1 次產(chǎn)生峰值的時(shí)刻是78 s，其峰值為22.2°，橫傾角來回振蕩波動(dòng)的時(shí)間是69 s。結(jié)果顯示，隨著轉(zhuǎn)速的增大，橫傾角第1 次產(chǎn)生變化時(shí)刻和發(fā)生峰值的時(shí)刻以及來回振蕩波動(dòng)的時(shí)間越來越小。橫傾角的峰值變化主要發(fā)生在潛艇上浮運(yùn)動(dòng)過程中，來回振蕩波動(dòng)主要發(fā)生在艇出水時(shí)刻。對(duì)比3.2 節(jié)有關(guān)上浮的時(shí)刻可知，無論轉(zhuǎn)速是多少，潛艇第1 次橫傾角的峰值均發(fā)生在潛艇上浮運(yùn)動(dòng)過程中，隨著轉(zhuǎn)速的增大，橫傾角的峰值并未呈現(xiàn)出明顯的變化規(guī)律。

4 結(jié) 語

本文基于體積力模型的虛擬槳盤模型，對(duì)不同轉(zhuǎn)速下潛艇的水下定深直航、六自由度上浮、出水以及水面航行進(jìn)行了數(shù)值模擬，并對(duì)比分析了不同轉(zhuǎn)速下潛艇螺旋槳與轉(zhuǎn)速的匹配及其姿態(tài)等，該方法可為潛艇的應(yīng)急上浮運(yùn)動(dòng)提供一定的參考。

通過對(duì)比分析不同轉(zhuǎn)速下艇體上浮運(yùn)動(dòng)過程中的姿態(tài)，可得出如下結(jié)論：隨著轉(zhuǎn)速的增大，潛艇獲得穩(wěn)定航速的時(shí)間越短，艇體在水下做定深直航運(yùn)動(dòng)的時(shí)間也越短；不同轉(zhuǎn)速下其上浮運(yùn)動(dòng)的時(shí)間相差不大；轉(zhuǎn)速越大，無論是縱傾角還是橫傾角，其傾角第1 次產(chǎn)生變化的時(shí)刻和發(fā)生峰值的時(shí)刻會(huì)越來越?。煌瑫r(shí)隨著轉(zhuǎn)速的增大，縱傾角的來回震蕩波動(dòng)時(shí)間會(huì)越來越長，而橫傾角的來回震蕩波動(dòng)時(shí)間則越來越短，且傾角的峰值隨轉(zhuǎn)速的變化并未呈現(xiàn)出明顯的變化規(guī)律。