于祥，胡開業(yè)

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

0 引 言

潛艇潛航狀態(tài)下產(chǎn)生的水動力尾跡主要有伯努利水丘（Benoulli hump）、開爾文尾跡、湍流尾跡、內(nèi)波尾跡和渦尾跡等。在潛艇非聲學探測領域，該類水動力信號可被遙感與雷達影像偵測[1-4]。其中開爾文尾跡與艇體上方伯努利水丘在遙感影像上識別度較高，潛艇通過增加潛深或降低航速等手段可有效抑制該類信號，從而顯著降低被探測的概率。

常規(guī)方法不易控制分層海洋中潛艇航行產(chǎn)生的內(nèi)波，因其存在時間久且傳播范圍廣，即使在惡劣天氣、雨、夜中，仍可通過星載或機載合成孔徑雷達（SAR）全天候廣泛地偵測潛艇位置，并通過散射影像推算出潛艇的運動狀態(tài)、潛深及所處的位置。Voropayev[5-6]早在1999 年就發(fā)現(xiàn)，當潛艇在連續(xù)分層流體中進行操縱機動航行時可產(chǎn)生大尺度的漩渦，該類由潛艇非定常機動航行誘發(fā)的水動力興波與海表波相互作用，在未來的潛艇非聲學探測領域有重要的學術(shù)與工程研究意義。目前關(guān)于潛艇尾跡特性的理論、數(shù)值模擬及試驗研究，幾乎均假設潛艇作勻速定深航行運動，對非定常與潛艇機動運動方面的研究極少。

本文將運用CFD 技術(shù)，對強分層海洋環(huán)境下潛艇非定常機動模式對海表及內(nèi)部流場的影響進行定量分析，研究潛艇運動狀態(tài)及加速度大小對其隱身性能的影響，為潛艇非聲學探測技術(shù)以及潛艇運動狀態(tài)推演提供相關(guān)理論依據(jù)。

1 理論模型

針對不同時刻、不同空間上的各物理量，CFD 數(shù)值模擬是通過計算域內(nèi)離散點變量來表示，以一定原則、方式為基礎，建立并近似求解相關(guān)方程組來獲取各離散點場變量關(guān)系[7-8]，密度分層流中流動過程遵循質(zhì)量守恒、動量守恒定律，基于溫度分層的密度躍層時，還需遵循能量守恒定律。在水動力CFD 數(shù)值模擬領域，對流體黏性的處理通常采用求解雷諾平均納維?斯托克斯（RANS）方程。

1.1 數(shù)值模擬控制方程

通過相關(guān)控制方程，來表示流體遵循的三大守恒定律。本節(jié)主要運用能量守恒、質(zhì)量守恒和動量守恒方程實現(xiàn)潛艇在不同海洋分層環(huán)境中航行尾跡及流場的數(shù)值模擬。

1） 質(zhì)量守恒方程。

任意流體微粒在單位時間內(nèi)，其質(zhì)量不隨時間變化，且流入、流出微元體中的質(zhì)量相同。也可認為相同時間內(nèi)微元體中增加的流體質(zhì)量等于流入該微元體的流體質(zhì)量。其表達式為

式中，u，v，w為速度矢量u在x，y，z這3 個軸方向上的分量。

2） 動量守恒方程。

流體所受合外力等于其單位時間內(nèi)的動能增加量。不可壓縮流體動量守恒方程為：

式中： ρ為流體密度；p為流體微團上的壓力； τxx為流體微團上不同位置不同切向的黏性應力；Fx為流體微團所受到的各方向的體積力。

3） 流體域體積（VOF）多相流模型。

不可壓縮流體，通過連續(xù)性方程?·u=0的VOF 方程為：

式中：u為 流體速度；Cf為單元中流體體積；Cu為單元體積。當計算的網(wǎng)格單元內(nèi)充滿該流體時，C=1；網(wǎng)格內(nèi)無該流體時，C=0；0

1.2 計算模型及分析方法

1.2.1 潛艇加減速模型

潛艇變速機動過程的加速度通常為變加速，此過程中潛艇螺旋槳有效推力、水動力參數(shù)等諸多因素決定加速度瞬態(tài)值，數(shù)值計算中對其求解極為復雜[9-10]，且許多參數(shù)在求解時無法獲取，所以本文潛艇加減速不考慮螺旋槳。按式a=(vj?vi)/tij求解潛艇在該過程中的平均加速度，tij為變速時間，vi為變速前速度，vj為變速后速度。變速過程中速度瞬態(tài)值為aT+v(i?1)，其中，v(i?1)為前一數(shù)據(jù)周期的速度，T為計算時間步長。直到潛艇速度加速或減速到vj時穩(wěn)定，加（減）速過程通過編程對潛艇運動狀態(tài)場函數(shù)進行自定義。

1.2.2 尾流場分析法

SAR 靈敏度由雷達波頻決定，將波頻設置低風速檔時可消除海表陣風對內(nèi)波影像的影響，即可從復雜海邊信號中偵測到潛艇航行生成的內(nèi)波。美國曾通過該方法有效偵測到前蘇聯(lián)潛艇在300 m 潛深航行時產(chǎn)生的內(nèi)波信號。結(jié)合理論研究與實驗觀測，Alpers 等[11-12]得到SAR 可偵測識別潛艇潛航產(chǎn)生的內(nèi)波與海表流場經(jīng)過波流調(diào)制后的輻聚輻散信號強度區(qū)間為10?3～10?4s?1（s 為時間單位秒，下文均采用此單位計量）。潛艇運動產(chǎn)生的內(nèi)波誘發(fā)海表流場的空間分布特征同內(nèi)波水面映波類似，在波峰前與波谷后出現(xiàn)輻聚流，波谷前與波峰后出現(xiàn)輻散流，即水面映波的波節(jié)處出現(xiàn)明顯的分散、聚集流動。

CFD 中流場輻聚輻散強度為對自由面水平速度場V求其散度，其表達式為：

2 算例分析

2.1 數(shù)值模擬驗證

本文首先驗證RANS 方程、SSTk-ω 湍流模型、VOF 法對潛艇潛航狀態(tài)尾流場模擬的可靠性及自由面捕捉的準確性。通過數(shù)值模擬Rankine ovoid model 在水下勻速潛航時的自由面興波特性，將結(jié)果與1966 年David Taylor 水池試驗[13]的興波縱剖面波高數(shù)據(jù)進行對比，驗證數(shù)值模擬潛艇航行的可行性與精度。

根據(jù)Rankine 卵形體幾何參數(shù)[13]，通過CATIA軟件對Rankine 卵形體進行三維建模（圖1）。本文參照泰勒水池潛艇拖曳航行試驗數(shù)據(jù)進行參數(shù)設置，將模型潛深設置為0.457 m，以勻速方式進行拖曳。選取航速3.048 m/s 工況進行相關(guān)數(shù)值模擬，計算域取去流段長度30.48 m，進流段長度1.219 m。

圖1 Rankine 卵形體三維模型Fig. 1 Three-dimensional model of Rankine ovoid

對該卵形體潛航流場運用STAR-CCM+軟件進行數(shù)值模擬方案驗證，計算域網(wǎng)格劃分選取軟件自帶的切割體網(wǎng)格模型，針對黏性底層對卵形體表面劃分邊界層網(wǎng)格，對海表波高、波長方向網(wǎng)格進行精細加密，使單個波長與波高有充足的網(wǎng)格數(shù)，網(wǎng)格劃分如圖2 所示。

圖2 卵形體計算域網(wǎng)格劃分Fig. 2 Meshing of computational domain for ovoid

當潛艇潛深1.5 ft，航速V=10.0 ft/s，流場趨于穩(wěn)定且結(jié)果收斂時，取海表縱剖面波高數(shù)據(jù)與泰勒水池浮標試驗數(shù)據(jù)進行對比。如圖3（右圖波高以英尺為單位）所示，本文數(shù)值模擬方法的海表縱剖面波高、波形、波長數(shù)據(jù)與泰勒水池試驗結(jié)果近乎一致，說明本文所采用的CFD 方法在模擬潛艇運動時對于捕捉海表“興波”特性是有效且精準的，該方法對于實際海洋中潛艇運動尾跡特性的預報可行。

圖3 數(shù)值模擬驗證Fig. 3 Numerical simulation verification

2.2 計算模型及計算工況

以國外某常規(guī)潛艇為計算模型（圖4），該艇長73.8 m，寬9.9 m，艇體為水滴線型，在艇艏導流罩下方設有艏舵，艇艉設有2 片水平艉舵與1 片垂直艉舵。

圖4 潛艇計算模型Fig. 4 Submarine computational model

潛艇定深變速運動數(shù)值模擬采用重疊網(wǎng)格技術(shù)，流場取頂層為15 m 厚度的空氣（海表空氣為海平面條件），中層淡水厚度取40 m，密度為997.56 kg/m3，自由面以下40 m 處為強分層內(nèi)界面。內(nèi)界面下方取厚度110 m、密度1 024.2 kg/m3的鹽水。計算域取來流段450 m 為入口，去流段1 000 m 處為壓力出口，流場寬度為側(cè)方500 m。網(wǎng)格采用切割體網(wǎng)格，對各界面網(wǎng)格進行垂向加密，確保波高范圍內(nèi)垂向網(wǎng)格30 個以上。同時為準確刻畫海表開爾文波波長與波形，對自由面網(wǎng)格進行V 形雙層加密，在潛艇近場波長方向設置網(wǎng)格數(shù)70 以上。網(wǎng)格總數(shù)約為2 500萬，各角度網(wǎng)格劃分視圖如圖5 所示。

圖5 潛艇計算域網(wǎng)格劃分Fig. 5 Meshing of computational domain for submarine

本文海洋密度分層環(huán)境數(shù)值模擬參數(shù)的選取根據(jù)中國近海相關(guān)水文數(shù)據(jù)[14-15]設置，將其合理簡化為強分層海洋模型。設置航速分別為12，16 和20 kn。本文各加減速工況在0～10 s 內(nèi)勻速直航，待流場穩(wěn)定后，在10～50 s 內(nèi)勻變速，50～60 s內(nèi)以第50 s 速度勻速直航，各工況具體航行參數(shù)的選取如表1 所示。

表1 計算工況Table 1 Calculation cases

2.3 計算結(jié)果及分析

針對12 kn 勻速航行潛艇，劃分了3 套粗細程度不同的網(wǎng)格，并針對相關(guān)物理量，用于驗證數(shù)值模擬過程中的網(wǎng)格無關(guān)性。3 套網(wǎng)格數(shù)量分別為885 萬（粗），2 500 萬（中）和4 500 萬（細），計算得到的興波如圖6 所示（下圖所有工況波高與速度場均以m 為單位）。通過分析發(fā)現(xiàn)，中網(wǎng)格和細網(wǎng)格興波差距最大值在5%內(nèi)，綜合考慮精度與計算耗時，最終選取中等網(wǎng)格方案。

圖6 不同網(wǎng)格數(shù)自由面興波Fig. 6 Free surface wave-making with different grid numbers

圖7 為t=60 s 時潛艇在不同航速與非定常航行工況下內(nèi)界面興波波形圖。勻速工況下，內(nèi)界面興波距離隨航速增加而顯著增加；當航速提升至20 kn 時，內(nèi)界面橫波逐漸清晰，波峰呈拖曳式傳播，衰減較慢且幅值隨航速增加顯著降低，波谷隨航速增加略微增加，其尺度同艇體尺度，狀態(tài)穩(wěn)定。加速工況B1 和B3 在t=60 s 時興波波谷幅值已略高于工況A3 時的情況，波峰幅值也高于工況A2。相比勻速工況，加速工況在內(nèi)界面波形與傳播規(guī)律均較為相近，均呈窄V 型傳播。在波峰尾端均存在漩渦，該漩渦尺度隨航速增加而增加，但加速時，內(nèi)界面波峰范圍內(nèi)幅值明顯強于勻速航行工況A2 和A3，且向后方衰減更慢。而減速時，工況B2 和B4 顯示艇體下方波谷幅值略小于勻速工況A2，后方波谷幅值僅在艇體近場略大，向后方衰減相比勻速、加速工況更快。其波峰幅值強于勻速16 kn 工況，波峰后方漩渦尺度也弱于勻速與加速工況。

圖7 內(nèi)界面興波Fig. 7 Wave-making of the internal interface

圖8 為t=60 s 時潛艇在不同勻速與非定常航行工況下自由面興波波形圖。潛艇開始運動時，潛艇周圍流體由于潛艇自身的體積排擠效應產(chǎn)生“伯努利水丘”，在艇艏與艇艉分別對應首波峰與尾波谷。隨著潛艇不斷運動，“伯努利水丘”不斷向后方及艇體兩側(cè)發(fā)展。低航速時（工況A1），艇體后方興波以橫波為主；隨著航速增加至16 kn，該工況下自由面橫波散波均較明顯；當潛艇航速增加至最高20 kn 時，自由面興波逐漸以散波為主。所有工況自由面興波均呈V 型，波谷值位置除工況B4 外，均位于潛艇指揮臺上方，勻速工況下興波波長、波高、傳播速度隨航速增加顯著增加，其興波散波包絡角隨航速增加緩慢增加。工況B1 顯示，當航速勻速增加4 kn，波谷值約為?0.392 m，較勻速16 kn 增加約2%，但波峰值降低約28%，該工況下艇體后方散波較16 kn 航速的情況尤為明顯。當航速由12 kn 加速至20 kn 時，自由面興波夾角小于工況A3，自由面興波以艇體近場“伯努利水丘”為主，艇體后方興波明顯較弱，其興波波高較工況A3 增加約2%，波峰波谷幅值均略有增加。當潛艇從16 kn 減速至12 kn時，自由面興波波高較工況A1 降低約21%；當從20 kn 減速至12 kn 時，其興波波高較工況A3 降低約55%。綜上，減速可在短時間內(nèi)將波高顯著降低。本文加速工況下，興波波長隨航速增加呈逐漸增加趨勢；減速工況下，隨航速減小而減小，波長及艇艏興波隨航速變化極為敏感。自由面興波夾角在加速時略微減小，但減速時卻明顯增加為寬V 型向艇體后方傳播，且減速工況下，興波由于前后期興波波速差距而疊加，顯得更加雜亂無章。

圖8 自由面興波Fig. 8 Free surface wave-making

為定量分析自由面興波特征，取y=0 處興波切面（圖9）。t=60 s 時，工況B1 潛艇近場自由面興波波形與航速16 kn 相似，最大波峰與航速16 kn一樣均位于第2 波峰處。波長較航速16 kn 時短，且后方興波波長不斷變化。潛艇最大波峰波谷所在位置相差較小，第3 波谷后較航速16 kn 波長減小較明顯，后方波形與航速12 kn 相似，且興波衰減同12 kn 航速相似；其余工況艇體近場自由面興波波高、波長與最終航速的勻速狀態(tài)較為接近，工況B1 和B2 顯示加速狀態(tài)下興波波長及波高隨潛艇航速的增加逐漸增加；工況B2 和B4 顯示減速工況興波波長逐漸減小，由于末端航速較低，初始高航速興波傳播速度快于末端低航速興波，從而產(chǎn)生興波疊加的現(xiàn)象，使興波波形無規(guī)律。綜上所述，定潛深狀態(tài)下進行潛艇變速機動航行時，艇體近場流體受到潛艇排開水的體積效應影響，而體積效應又受到潛艇航速、潛深影響，使得近場自由面興波波長、波高對航速變化較敏感，可從該角度判斷潛艇的運動狀態(tài)。

圖9 自由面興波截面Fig. 9 The cross-section of free surface wave-making

為體現(xiàn)自由面粗糙度，采用自由面三維輻聚輻散強度來表示。如圖10 所示，工況A1～A3，B1～B4 對應自由面輻聚輻散強度范圍分別為?0.045 4～0.032 s?1，?0.066～0.043 s?1，?0.128～0.083 s?1，?0.068～0.047 7 s?1， ?0.042 4～0.025 3 s?1，?0.126～0.088 s?1，?0.058 2～0.045 4 s?1。上述數(shù)據(jù)表明，無論勻速或加速、減速工況，其自由面輻聚流峰值均高于輻散流峰值，兩峰值差距隨航速增加而增加，且加速時差距明顯增加，減速時差距減小，其峰值均位于最大波谷兩側(cè)。工況B1 航速由12 kn加速至16 kn 時，其自由面粗糙度較勻速16 kn 工況略有增加，加速工況B3 的粗糙度也略高于勻速20 kn 的工況。相比減速工況B2，在第60 s 時其自由面粗糙度甚至小于勻速12 kn 的工況，且減速工況B4 由20 kn 高航速減速至12 kn 時，較勻速20 kn 工況，其自由面輻聚流強度降低約55%，輻散流強度降低約45%。綜上所述，在潛深不變的條件下，加速時會顯著提升潛艇被探測到的概率，減速可明顯降低潛艇被探測到的概率，潛艇在16 kn 航速以內(nèi)進行減速時，其隱身效果明顯較好。

圖10 自由面輻聚輻散強度Fig. 10 Free surface convergence-divergence intensity

為分析潛艇運動狀態(tài)改變對其尾流場的影響，在y=0 m 處取y方向為法向的縱截面，統(tǒng)一速度場區(qū)間為0～16 m/s，如圖11 所示。圖11（a）～圖11（c）顯示當潛艇航速逐漸增加時，潛艇指揮臺處對自由面流場速度擾動較大，且隨航速擾動逐漸增加，航速20 kn 時擾動尤為明顯，此時指揮臺對艇體后方數(shù)百米擾動均較大，衰減也較為緩慢。當工況B1 潛艇從12 kn 加速至16 kn 時，艇尾流場擾動甚至強于勻速16 kn 工況，且指揮臺對自由面速度場擾動也明顯強于工況A2。工況B2 從16 kn 減速至12 kn 時，指揮臺對潛艇近場自由面速度場擾動弱于勻速12 kn 工況，但遠場（艇身200 m 后）擾動略高，艇尾速度擾動略高于勻速12 kn，但明顯小于工況A2。工況B3 和B4的艇尾流場擾動距離相當，但工況B3 在近場擾動與工況A3 相當，末端由于初始低航速擾動明顯較小。而工況B4 與其相反，近場擾動略低于勻速12 kn 工況，遠場擾動明顯強于勻速12 kn。圖中明顯可見，相同加速減速區(qū)間內(nèi)，加減速狀態(tài)艇尾流場漩渦垂向尺度明顯高于加速狀態(tài)。

圖11 垂向速度場分布Fig. 11 Distribution of vertical velocity field

為研究潛艇運動狀態(tài)對內(nèi)波速度場的擾動，取內(nèi)界面速度場分布（圖12）。所有工況速度場擾動均在潛艇艏艉近場擾動極為劇烈。勻速狀態(tài)時，內(nèi)波水平方向傳播距離隨航速增加逐漸增加，航速為16 kn 時，艇體后方橫波逐漸明顯，且橫波隨航速的增加越發(fā)明顯。加速狀態(tài)時，潛艇艇身附近水平傳播的范圍略微增加，遠場傳播距離均較小。減速時，潛艇近場內(nèi)波水平方向的傳播范圍明顯高于其勻速最高航速工況，相比加速工況，減速時內(nèi)界面橫波呈前方弱、后方強的趨勢，該特征與加速時相反?？梢姺嵌ǔC動航行對近場內(nèi)波水平方向傳播較大，結(jié)合圖11 可見減速狀態(tài)尾流塌陷非常明顯。綜上所述，減速狀態(tài)雖然可顯著降低波高及海表粗糙度，但對于海洋內(nèi)部的擾動較為劇烈，且存在的范圍較廣。

圖12 內(nèi)界面速度場分布Fig. 12 Distribution of internal interface velocity field

3 結(jié) 論

本文通過CFD 方法，采用RANS 方程、SSTk?ω湍流模型驗證了數(shù)值模擬潛艇潛航狀態(tài)興波的準確性，在此基礎上通過歐拉多相流對強分層海洋模型進行兩層流分層，對潛艇在強分層海洋環(huán)境下進行勻速直航、變速機動的數(shù)值模擬。以潛艇非聲學探測與隱身為目的，研究3 種不同航速區(qū)間內(nèi)潛艇不同運動狀態(tài)下自由面、內(nèi)界面興波尾流場特征，并從遙感（輻聚輻散強度）角度對潛艇的隱身效果進行評估。相比常規(guī)速度場分析，本文從多個角度及層面對機動航行的潛艇展開較為全面的分析研究。具體結(jié)果如下：

1) 數(shù)值模擬的Rankine 卵形體潛行狀態(tài)興波結(jié)果與試驗相比誤差較小，本文各物理模型及網(wǎng)格劃分方法滿足潛艇潛行狀態(tài)自由面興波的精度要求，為求解不可壓縮密度分層流動提供了重要的方法。同時本文對實尺度潛艇數(shù)值計算網(wǎng)格無關(guān)性做出驗證，可為進一步數(shù)值模擬潛航體機動航行分層流動提供必要的數(shù)值計算驗證。

2) 強分層海洋環(huán)境中潛艇進行加速機動時，隨著航速的增加，首波峰、波谷隨之迅速增加，對應興波波長逐漸增加，但相對波高增加速度較慢，可見艇體的體積效應激發(fā)的內(nèi)波在自由面首波峰、波谷處響應較快，其體積效應受航速影響較大。無論加速還是減速工況，內(nèi)界面興波波峰幅值均有所降低，其內(nèi)界面興波波峰相比勻速狀態(tài)衰減較慢，減速時衰減較快。當加速度增加一倍時，輻散流強度增加約1 倍，與勻速20 kn 相比，航速從12 kn 加速至20 kn 高航速時更容易被發(fā)現(xiàn)，因此加速狀態(tài)可顯著增加潛艇被探測到的概率。加速狀態(tài)對整個尾流速度場擾動甚至強于其勻速最高速度航行工況。

3) 當潛艇以較低航速（16 kn）減速可顯著降低自由面興波幅值及其海表粗糙度，加速與減速時其波峰峰值相對潛艇位置不同，減速時位于散波波峰首部，加速時位于散波波峰尾部。以較高航速減速時，雖然相比最高航速（20 kn）時其興波波高、自由面粗糙度會顯著降低，但相比末端勻速12 kn 航速，其被探測的可能性仍較高。因此，當潛艇以較低航速（16 kn）減速可有效降低海表粗糙度。減速時對整個尾流速度場擾動有顯著降低。無論加速或減速工況，艇體近場內(nèi)波傳播距離均有所增加，減速狀態(tài)最明顯，且對應近場尾流塌陷最為劇烈。