附體對分層流中潛艇水動力特性的影響

劉雙， 何廣華,， 王威， 高云

(1.哈爾濱工業大學 機電工程學院, 黑龍江 哈爾濱 150001;2.哈爾濱工業大學(威海) 船舶與海洋工程學院, 山東 威海 264209)

0 引言

潛艇是各國發展軍事力量不可或缺的一環，研究潛艇的水動力特性對建立潛艇作戰系統至關重要。實際海水中存在的溫、鹽差會使水體出現分層現象，潛艇航行于多層流之間時會產生內波[1]，針對潛艇航行興起內波尾跡的研究具有重要意義。

各國學者采用理論及數值方法對航行體源致內波進行了研究[2-4]。在興波尾跡研究方面，Yeung等[5]通過求解Green函數，探討了運動潛體產生的表面波模式及內波模式對自由液面和界面處興波的影響；趙先奇等[6]對細長體在三層流體中運動生成內波的問題進行了模型實驗研究；Chang等[7]研究了雙層流中航行潛艇的尾跡特性；Chomaz等[8]研究了分層流體中球體的近場尾流情況。在分層流中物體受力研究方面，Motygin等[9]基于勢流理論探討了二維物體分別在上層、下層流體中運動所產生的波浪阻力公式；Grue[10]對“Fram”號極地船航行過程進行研究，得出的受力變化規律與Ekman[11]一致；勾瑩等[12]采用模型實驗方法研究了箱型結構在雙層流中拖動時的阻力特性，并與單層流中拖航阻力實驗結果進行了對比研究；Wei等[13]、Wang等[14]針對內孤立波與水下潛體相互作用問題進行了實驗研究。在潛艇附體研究方面，Bull[15]通過改變湍流模型、求解器及網格形式，系統地驗證了SUBOFF潛艇尾流場數值模擬方法；Bensow等[16]采用大渦模擬(LES)方法計算了SUBOFF潛艇的黏性流場；Posa等[17]研究了雷諾數對全附體SUBOFF潛艇尾跡特征的影響；Yao等[18]采用LES及邊界元方法對全附體潛艇周圍的流場及激流噪聲進行了研究；Seil等[19]研究了附體及其加裝方式對BB2型潛艇水動力載荷的影響。

綜上可知，尾跡研究常常僅關注界面波模式，對于航行體受力，尤其是航行體外形對其自身水動力特性的影響研究較少。針對密度分層流這一特殊海況，分析各個附體對潛艇水動力特征帶來的影響顯得十分必要。

本文采用黏性計算流體力學(CFD)理論，基于工程物理仿真軟件STAR-CCM+平臺建立了密度分層流數值仿真模型。首先進行收斂性驗證，并將數值結果與實驗結果進行對比，驗證本方法的準確性；再以SUBOFF潛艇為基礎研究其在裸艇、僅有駕駛艙附體、僅有尾翼附體及全附體工況下以不同航速在密度分層流中航行時的興波及阻力特征，分析各附體工況下航速對潛艇所受阻力的影響以及阻力系數峰值的出現與艇艏、艉興波干擾的關系；并結合壓力分布探討附體存在帶來的影響。本模型具有較好的計算精度，可為潛艇在密度分層流中航行時的水動力性能預報提供新的手段，同時對潛艇設計優化與規避策略提供參考。

1 數值計算模型

1.1 控制方程

本模型基于雷諾平均(RANS)方程，不可壓縮牛頓流體的控制方程為

(1)

(2)

(1)式為連續性方程，(2)式為動量方程，二者均為時均計算處理后的形式。

1.2 湍流模型

為使(2)式封閉，需要引入湍流模型。采用Rea-lizablek-ε湍流模型(k為湍動能，ε為耗散率)，該模型穩定性良好，壓力梯度求解精度高，工程應用比較廣泛，其湍動能及耗散率輸運方程為

(3)

(4)

式中：xt表示在x軸方向上對時間t求偏導；μ為分子擴散所造成的動力黏性；μt為湍流黏性系數；σk和σε分別為湍動能k和耗散率ε的湍流普朗特數；Gk為平均速度梯度產生的湍動能；Gb為由浮力產生的湍動能；YM表示湍流脈動膨脹對總耗散率的影響；C1=max[0.43,η/(η+5)]，系數η=Sk/ε，S為平均應變率；C2和C1ε為常數；C3ε為浮力對耗散率影響的函數，C3ε=tanh|up/uv|,up為平行于重力方向的速度分量，uv為垂直于重力方向的速度分量。參數取值為C2=1.9、C1ε=1.44、σk=1.0、σε=1.2，在文獻[20]中已證實是可靠的。

1.3 平面進行波理論

根據平面進行波理論分析興波阻力特性，對水下航行潛艇的快速性研究具有參考價值。寬度為b、長度為2λ的封閉區域中平面進行波波能Eb為

(5)

式中：g為重力加速度；H為波高。

平面進行波理論中興波阻力系數Cw可表示為

(6)

式中：Rw為興波阻力；v為航速；mL表示興波長度，m為系數，與船型和弗勞德數有關，L為艇長；A為濕表面面積；C、D為常數；Fr為弗勞德數。

由于cos(2πmL/λ)值是在1.0和-1.0之間變動，因此興波阻力系數Cw-Fr曲線存在波動。

1.4 流體體積方法

通過流體體積方法進行界面捕捉，引入變量相體積分數αq，定義每個控制單元中第q相的體積分數為

(7)

式中：Vq為第q相在控制單元中所占的體積；V為控制單元的體積。當αq=0時，表示在控制單元內不存在第q相流體；當αq=1時，表示在控制單元中充滿第q相流體；當0<αq<1時，說明在控制單元中存在第q相流體和其他相流體的界面。

2 數值模型與驗證

2.1 計算模型可靠性驗證

分層流中潛艇阻力及興波特性的CFD研究較少，研究基于STAR-CCM+軟件進行仿真，首先對數值模型進行驗證。根據文獻[12]中的實驗裝置參數進行建模如圖1所示。其中水深為0.6 m，淡水在上層，密度ρw=997.56 kg/m3；下層是密度ρs=1 024.00 kg/m3的鹽水；上層淡水厚hw為0.30 m，下層鹽水厚hs為0.30 m. 實驗模型為60 cm×45 cm×35 cm的方型浮箱，吃水深度定為0.2 m，速度分別為0.14 m/s、0.16 m/s、0.18 m/s、0.20 m/s、0.24 m/s.

圖1 實驗裝置示意圖[12]Fig.1 Schematic diagram of test device[12]

與相關實驗結果的對比見圖2，其中數值模擬結果為箱體在雙層流中拖航時所受的總阻力，用Rtot表示。由圖2可以看出，數值模擬結果與文獻[12]中的實驗數據基本吻合，說明本文所采用的CFD計算模型可靠。

圖2 本文數值模擬結果與實驗數據[12]對比Fig.2 Comparison of the simulated results in the present paper with the experimental data in Ref.[12]

2.2 數值模擬相關設置

采用對稱面條件來減少計算資源的消耗，提高計算效率，故模擬得到的是半艇阻力值，需乘以2來展示全艇阻力。為了防止數值模擬過程中回流的產生，適當延長計算域長度，取艇前3倍艇長，艇后12倍艇長。潛艇表面邊界條件類型為壁面，物面條件為無滑移?；跉W拉多相流模型，通過用戶自定義場函數(UDF)的方式分別給定空氣、淡水、鹽水3相的密度及初始體積分數，使流體域分層。入口邊界條件類型為速度入口，在此處指定流體速度；出口邊界條件類型為壓力出口，在此處通過UDF方式給定分層流體的壓力分布。

對各流體相分界面及潛艇附體處網格進行加密以準確捕捉潛艇在分層流中航行時的流場特性。此外，使用切割體網格并結合棱柱層網格生成器來處理邊界處的網格，其中棱柱層厚度為0.018 m，棱柱層數為6，增長率為1.2. 調節網格增長率等參數控制網格質量，如圖3所示，通過計算機(中央處理器i7-8700,主頻3.19 GHz)進行計算，一個算例耗時30 h左右。

圖3 潛艇周圍網格分布Fig.3 Grid distribution around submarine

研究采用的潛艇模型為標準SUBOFF潛艇模型，分別選擇4種附體工況(裸艇、僅有駕駛艙、僅有尾翼、全附體)，其中全附體SUBOFF模型見圖4，具體尺寸如表1所示，潛艇潛深為h，即潛艇重心到靜水面的垂向距離。

圖4 SUBOFF全附體潛艇模型Fig.4 Full-appendage SUBOFF model

表1 潛艇尺寸參數

計算域中流體密度沿垂向變化，分層情況如圖5所示。SUBOFF潛艇以不同速度航行于上層即淡水層，潛深h=0.12L,潛艇重心距淡水- 鹽水交界面的距離同樣為0.12L.

圖5 分層示意圖Fig.5 Schematic diagram of density stratified fluid

2.3 網格及時間步長收斂性驗證

采用加裝全附體的潛艇以航速Fr=0.5進行收斂性研究，由于網格與航速、波長、波陡等諸多因素有關，參考文獻[21]的建議，并經過一系列數值模擬最終找到適合本文的經驗方式。數值模型網格劃分的宗旨是：保證數值模型能夠在潛艇附體周圍、交界面等處有足夠的網格密度，使單位興波波長內大于80個網格、單位波高內大于15個網格，用以捕捉潛艇在各交界面處的興波特征。通過改變網格的基礎尺寸來控制網格數量，取時間步長為0.01 s，結果如表2所示。其中誤差是指采用某一網格量計算結果與采用上一個較少網格量所得結果之間的差別；表2中Rf代表摩擦阻力,在STAR-CCM+軟件中，摩擦阻力通過(8)式計算：

(8)

式中：Tn為面n上的應力張量；an為面網格面積矢量。

表2 網格收斂性分析

由表2可以看出，當網格基礎尺寸由網格3降低到網格4時，阻力計算結果幾乎沒有變化，證明可采用網格3方式進行高效的CFD研究。

為找到快速、穩定的時間步長方案，參考文獻[21]中的建議，時間步長Δt按照0.005L/U～0.01L/U，并考慮庫朗數的要求來給出時間步長的預估值，再通過大量的數值實驗最終確定時間步長。表3為采用網格3網格方案進行的時間步長收斂性研究，可見隨著時間步長降低，計算逐漸達到收斂。當時間步長達到0.01 s時，再降低時間步長到0.005 s對計算結果的影響不大，故本研究中均選擇時間步長為0.01 s.

表3 時間步長收斂性分析

3 數值結果與分析

3.1 不同附體類型對潛艇快速性的影響

3.1.1 阻力變化

圖6 阻力隨航速的變化Fig.6 Change of resistance with navigational speed

模擬中計算域如圖5所示，淡水層總厚度為0.24L，密度取為997.56 kg/m3，鹽水層厚為0.67L，密度為1 020.00 kg/m3. 根據文獻[20]，潛艇離交界面越近、興波情況越明顯，故使潛艇位于潛深h=0.12L處，同時潛艇中心距內波面距離同樣為0.12L，以便觀察潛艇附體改變時其在兩個界面處興波的變化。針對每一種附體工況均模擬了Fr為 0.3、0.4、0.5、0.6、0.7、0.8 6種速度，在研究附體影響同時分析潛艇水動力性能隨航速的變化。圖6為各附體工況下總阻力及摩擦阻力隨航速的變化。

由圖6可知，隨著航速增加，潛艇的總阻力和摩擦阻力不斷增加。低速時，附體對阻力的影響較??；高速時，附體對阻力有增大的作用；尾翼對摩擦阻力的影響要略大于駕駛艙，但二者對總阻力的影響幾乎一致。

3.1.2 阻力系數變化

密度分層流中各附體工況下潛艇興波阻力系數隨航速的變化趨勢如圖7所示。由圖7可以看出，曲線隨航速的增加出現先上升、后下降的趨勢，這是由于在興波阻力系數計算(6)式中存在一個余弦項，此項會根據潛艇艏、艉興波的相互干擾情況在-1～1之間變化。低速狀況下，僅加裝尾翼時興波阻力系數最低，加裝全附體工況時系數最高；高速時各附體對阻力系數的影響不大；整體來看，尾翼的存在會使興波阻力系數降低。

圖7 阻力系數隨航速的變化Fig.7 Change of resistance coefficient with navigational speed

各附體工況下：在Fr=0.5附近興波阻力系數均最高，說明此時艇艏興波與艇艉興波發生了不利干擾；Fr=0.5時僅加裝尾翼工況的阻力系數最低，而僅加裝駕駛艙的潛艇阻力系數最高。

可見，附體的不同會對阻力和阻力系數帶來一定影響，且這種影響與航速有關。

3.2 不同附體存在對潛艇興波特性的影響

潛艇在密度分層流中航行時，會在自由液面及內波面處產生興波，興波尾跡特征可以作為水下潛體探測的重要參考，興波會隨著潛艇航速的變化而改變。本節取Fr為0.3、0.5、0.7工況研究不同附體工況下潛艇以不同航速航行時在兩個分層界面上的興波及速度特性。

潛艇在裸艇與全附體工況下航行時在自由液面及內波面處的興波尾跡隨航速變化情況如表4和表5所示。由表4和表5可見：隨著航速的增加，裸艇與全附體自由液面均呈Kelvin波形，興波范圍越來越廣，波長逐漸增加；此外，相比低速(Fr=0.3)以及高速(Fr=0.7)情況，中速(Fr=0.5)附近的波峰峰值最大。內波波形與表面波存在較大區別，興波將不呈現Kelvin波形，而是波峰集中在整個興波區域的中間部分，并呈現峰谷交替情況，波形以狹長的V形向后方延伸，在整體興波的尾部出現環狀的興波區域。且自由液面處第1個興波必為波峰，稱為伯努利水丘，而內波面的第1個興波必為波谷。

表4 裸艇與全附體工況下自由液面興波云圖隨航速變化

表5 裸艇與全附體工況下內波面興波云圖隨航速變化

隨著航速的增加，內波面興波波長也會逐漸增加，并在Fr=0.5時波峰峰值較高。潛艇興波會消耗能量，波的能量與波長的平方呈正比，故隨著速度增加，潛艇所受阻力會出現上升的趨勢；在Fr=0.5時峰值較高，說明此時潛艇艏、艉興波發生了不利干擾，這與興波阻力系數曲線所展示的情況一致。由于淡水、鹽水之間的密度差遠小于淡水與空氣之間，且僅為其千分之一量級，故很小的擾動即會激發大振幅的內波，且其密度梯度小、恢復力弱，所以內波周期長、波長大，整體波形對速度變化的響應遲緩。

表6和表7分別為裸艇及全附體工況時自由液面和內波面上沿潛艇前進方向的水質點速度分布情況，與表4、表5一致，上側為裸艇工況、下側為全附體工況，潛艇航行方向為x軸負向。由表6和表7可見：裸艇及全附體自由液面與內波面速度云圖均呈現V形；自由液面處水質點速度要高于內波面，說明內波面上水質點運動對潛艇航速的改變不敏感。

表6 裸艇與全附體工況下自由液面水質點速度云圖隨航速變化

表7 裸艇與全附體工況下內波水質點速度云圖隨航速變化

將裸艇與全附體工況對比可知，潛艇航行于密度分層流中時，無論從興波尾跡、水質點速度的分布形狀還是從波形中峰、谷及波長的數值來看，附體存在帶來的影響均較小。經過分析研究，駕駛艙及尾翼附體單獨存在對潛艇興波尾跡及水質點速度分布同樣影響極小，故不在展示其興波及速度分布情況。

3.3 不同附體存在對潛艇上、下方波面抬升的影響

表8為潛艇在加裝不同附體工況下以航速Fr為0.3、0.5、0.7航行時，其上方自由液面及下方內波面在y=0 m剖面上的波切線。表8中兩條虛線之間為潛艇在x軸方向的位置，潛艇沿x軸負向航行。

由表8可以發現，隨著航速增加，興波波長不斷加大。相比內波面處，自由液面處的波形清晰，變化比較明顯，這與3.2節中闡述的內波面波形變化緩慢原因一致；自由液面上第1個波形為波峰，內波面上第1個波形為波谷；對比不同附體工況下的結果發現，低速時駕駛艙的存在會使艇艏上方波形產生波動(表8中方框)，高速時幾乎沒有影響；整體來看，附體改變對各個速度下交界面處興波的波面抬升影響不大。

表8 各附體工況不同航速下興波波面抬升

圖8 不同航速下各附體工況潛艇表面壓力分布Fig.8 Surface pressure distribution of submarine with different appendages at various speeds

3.4 不同附體存在對潛艇表面壓力分布的影響

在Fr分別為0.3、0.5、0.7時加裝不同附體工況下的潛艇表面壓力分布情況如圖8所示。

由圖8可見：潛艇表面壓力隨著速度的增加而逐漸加大，其中潛艇艏、艉壓力變化較大，而主體部分的壓力幾乎不變；當存在附體時，潛艇表面曲率突變處，如駕駛艙及尾翼附近(圖8中方框)的壓力會產生較大波動，且加裝全附體時的壓力要明顯高于其他附體工況。分析原因是由于加裝附體后潛艇表面出現拐點，改變了局部流體的速度，進而引起壓力變動。

下面將結合具體壓力數值進行定量分析，Fr=0.5時加裝全附體潛艇的表面壓力及上方自由液面興波如圖9所示，由3部分構成，最下方是潛艇位置，中間是潛艇表面壓力，最上方為潛艇航行時其中縱剖面與自由液面波切線上的波面抬升。

圖9 潛艇表面壓力分析圖Fig.9 Analysis diagram of submarine surface pressure

在水中航行的潛體，其壓力組成可分為靜壓、動壓兩部分，其中：靜壓是垂直于潛艇表面以上單位流體所具有的靜壓能，可以通過公式p=ρgh來計算；動壓是由局部速度變化所引起，可由總壓減去靜壓得到。表9為各個分析點處的壓力值。

表9中E、J為潛艇艏、艉兩點，其壓差為3 264.5 Pa，這是形成潛艇壓差阻力的主要原因，由圖8可知，此數值會隨著潛艇航速的增加而不斷增加。為便于對比，后續的壓力計算值將以此數值為基準，用百分比形式展示。F點靜壓與G點靜壓的差值為潛艇艏、艉壓差的40.5%，但總壓差值卻是潛艇艏、艉壓差的272.2%，說明動壓的變化占較大比例；再對比H、I兩點，靜壓差為16.5%，而總壓差為187.9%，H、I兩點在潛艇上處于同一水平位置，上方波面抬升帶來的靜壓變化很小，說明主要是由于動壓的變化導致了總壓的差別。分析以上現象的原因是由于F、I點分別處于駕駛艙及尾翼的拐點處，附體帶來的曲率突變導致此處局部流體速度改變，動壓急劇變化，進而影響了各點的壓力分布。

表9 各分析點壓力值

4 結論

本研究基于黏流理論建立了一種可用于求解密度分層流中水下潛體阻力及興波尾跡的CFD模型，可較好地模擬出潛艇在多個附體工況下以不同航速航行時在自由液面及內波面處的興波尾跡特性以及阻力特征。全面地分析了附體及航速的改變對潛艇水動力性能的影響，并結合壓力變化探討了阻力等參數改變的原因，得出結論如下：

1)潛艇的興波特性隨航速的改變發生較大變化。由于潛艇艏、艉興波產生了不利干擾，使得在Fr=0.5附近時潛艇的興波阻力系數最大。

2)在密度分層流中，各附體對潛艇水動力性能的影響與航速有關。阻力方面：高速時，附體對潛艇阻力有增大的作用，尾翼對摩擦阻力的影響要略大于駕駛艙；當航速較低時，僅加裝尾翼附體興波阻力系數最低。

3)興波方面：駕駛艙對于波面抬升的影響要高于尾翼。整體來看，附體對潛艇興波特性的影響不大。

4)附體的存在所帶來的潛艇表面曲率突變會產生較大的壓力變動，其中動壓的改變占主導作用。

5)本文可為密度分層流中潛艇水動力性能研究提供一種有效的數值模擬方法，并對潛艇附體的設計優化提供參考。