深水密度層航行潛艇興波阻力的影響分析

劉雙， 何廣華, , ， 王威， 潘雁甲

(1.哈爾濱工業大學 機電工程學院，黑龍江 哈爾濱 150001； 2.山東船舶技術研究院，山東 威海 264209； 3.哈爾濱工業大學(威海) 船舶與海洋工程學院，山東 威海 264209)

密度分層廣泛存在于海洋環境中，其對潛艇航行性能會產生較大的影響[1]。密度分層流中航行的潛艇會同時在靜水面及分層界面興起波浪，興起波浪所需的能量以潛艇遭受興波阻力的形式展示出來。興波阻力為剩余阻力的主要成分，而快速性是潛艇設計、航行過程中著重需要考慮的性能。故研究潛艇在密度分層流這一特殊海況下航行時的受力特性具有十分重要的意義。

針對密度分層流中航行體的受力特性，各國學者進行了大量研究。Ekman[2]研究了由密度分層造成的“死水”現象，該現象在Grue[3]的研究中得到了進一步證實。此外，Hudimac等[4]建立了用于求解薄船表面波及內波興波阻力的公式。Motygin[5]基于勢流理論探討了二維物體分別在分流體中運動所產生的波浪阻力公式；Miloh等[6]根據格林函數法得到了半潛扁球體、船體所受的波阻；Esmaeilpour[7]采用CFD方法對密度分層流中艦船受力問題進行了研究，其他還有一些類似的研究[8-9]。此外，一些學者針對密度分層流中航行體運動及受力問題進行了試驗研究[10-12]。

綜上可知，對于密度分層流中潛艇航行性能研究大多關注其生成的內波模式，且對于分層流體存在引起潛艇受力的研究較少，尤其是位于深水層航行的模式(即在內波交界面以下航行)。由于海洋環境是多變、復雜的，對深水層航行潛艇的研究有著十分重要的意義，故本研究著重探討潛艇航行于深水層時的受力特性。

本研究所采用的數值模型基于RANS方程，采用Simple算法求解壓力-速度耦合項，采用Realizablek-ε湍流模型，用于密度分層流中求解潛艇水動力性能。

1 數值模型設置及網格劃分

本模型基于水動力分析軟件STAR-CCM+，求解雷諾時均的N-S方程，引入Realizablek-ε湍流模型對控制方程進行封閉。不可壓縮流體的控制方程及湍動能輸運方程可參考文獻[13]中的描述。

由于所求解問題及結構的對稱性，故采用半潛體以提高計算效率。計算域長為15L，寬為3.5L，其中L為潛艇長度。基于歐拉多相流模型，通過用戶自定義場函數的方式分別給定空氣、淡水、鹽水三相的密度(空氣密度ρair=1.184 15 kg/m3、淡水密度ρ1=997.561 kg/m3、鹽水密度ρ2=1 020 kg/m3)及初始體積分數。計算域示意圖見圖1，潛艇位于鹽水層中，圖中U為潛艇航速。d1指潛艇重心到內波交界面的垂向距離。h1為淡水層流體的厚度；ρ1、ρ2分別為淡水、鹽水層流體的密度。為保證有足夠的水深，本研究中鹽水層的厚度固定為5L。

圖1 分層示意Fig.1 Diagram of density stratified fluid

潛艇模型為SUBOFF標模，總長L=3 m，模型見圖2，具體尺寸見表1。

圖2 SUBOFF潛艇模型Fig.2 Submarine model of SUBOFF

表1 潛艇尺寸參數Table 1 Size parameters of SUBOFF

網格劃分采用切割體網格，可通過提供最大/最小網格單元尺寸來控制網格上、下限，能夠生成較高質量的網格，對交界面以及潛艇表面曲率變化處進行網格加密以準確捕捉流場特性，在潛艇前端進行網格定位。選取網格基礎尺寸為0.056 m，調節網格增長率、棱柱層厚度等參數控制網格質量，經過多次數值實驗，最終確定網格可保證對潛艇興波的較好捕捉，潛艇周圍計算域網格見圖3。此外，時間步長選取為0.01 s，可在保證計算收斂的基礎上盡量節約計算時間。

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

2 結果與分析

2.1 航速對潛艇受力的影響

采用A、B 2組工況來研究航速對潛艇受力的影響，計算域示意圖見圖1。由文獻[13]可知，潛艇越靠近交界面，其興波越劇烈，同時對潛艇受力的影響也越大。因本研究探討的是潛艇在深水層中航行時的受力特性，故取一偏小的淡水層厚度h1，以使本研究中的現象更加明顯。表2為工況A、B的模擬參數設置，潛艇均位于鹽水層中。

表2 模擬參數設置Table 2 Setting of the simulation

工況A中的潛艇距離2個交界面均較近，工況B中的潛艇與2個交界面的距離相比于工況A均較遠，兩者相比可觀察深水層中航行潛艇在靠近與遠離交界面時的阻力性能隨航速變化情況。本數值模型中通過式(1)計算總阻力，通過式(2)計算摩擦阻力，總阻力減去摩擦阻力即為剩余阻力，剩余阻力是由興波阻力以及粘壓阻力組成的，無法通過具體的公式將興波阻力分離出來。本研究中對于興波阻力的計算通過數值手段，即當每個算例計算完成后，將計算域拉大，使潛艇遠離兩交界面，保證此時潛艇在兩交界面上無興波，此時計算的總阻力即為除去興波阻力之后的結果，進而2次計算的總阻力相減即可得出興波阻力。

(1)

(2)

式中：Tf為作用于面f上的剪切應力；Pf為作用于面f上的壓力；af為面網格面積矢量。

2.1.1 航速對阻力的影響

圖4為A、B 2種工況下潛艇受力變化情況。由圖4(a)可知，2種工況下，隨著航速增加，潛艇的總阻力均呈現不斷上升的趨勢。整體來看，工況A中潛艇所受總阻力要大于工況B，在中速階段(0.4≤Fr≤0.8)，工況A的總阻力與工況B中總阻力相差較大，同時工況A的阻力增加更加平緩；在Fr<0.5與Fr>0.8時，兩者的阻力差別不大且增速基本一致。

由圖4(b)可見，隨著航速增加，2種工況中潛艇的摩擦阻力均不斷上升，且兩者摩擦阻力曲線基本重合，這是由于摩擦阻力主要由流體的密度及航速決定，與潛艇所處位置關系不大，也同時說明2種工況下的總阻力差別主要是由剩余阻力導致的，而剩余阻力中的主要部分為興波阻力。

由圖4(c)可知，當潛艇在鹽水層中航行時，2種工況下，興波阻力系數曲線均會呈現先上升后下降的走勢，曲線的峰值點均出現在Fr=0.6附近。在平面進行波理論中，興波阻力系數可近似表達為：

(3)

圖4 各工況下阻力隨航速的變化Fig.4 Resistances of submarine at different speeds under different conditions

式中：Cw為興波阻力系數；Rw為興波阻力；ρ為流體密度；λ為波長；mL為興波長度；Fr為長度弗勞德數(Fr=U/SQRT (gL))；g為重力加速度；S為濕表面面積；C、D為常數。

可見式(3)中存在余弦項cos(2πmL/λ)，其數值將根據潛艇艏、艉興波干擾的情況從-1～+1不斷變化，這導致興波阻力系數曲線將存在峰值點。結合圖4(c)說明當Fr=0.6附近時，潛艇艏、艉興波發生了不利干擾，而由前述研究得知，近水面航行潛艇的興波阻力系數曲線在Fr=0.5附近出現峰值，說明潛艇航行于深水層流體中會使興波阻力系數峰值點向高速推移。

此外，不同航速下工況A的興波阻力系數均大于工況B，這主要是由于工況A中潛艇距離交界面較近，會在交界面上激起較大的興波，進而產生了較大的興波阻力。且兩者之間的差別在0.4≤Fr≤0.8的中速階段較大，而低速及高速時差別不明顯，這與圖4(a)中總阻力的變化一致，進一步說明了：深水層航行潛艇總阻力的變化由興波阻力主導。

再回到圖4(a)中潛艇所受總阻力的變化趨勢，結合圖4(b)、(c)分析：當航速較低時，摩擦阻力與興波阻力均不斷增加，故總阻力不斷上升；當航速超過Fr=0.6附近后，興波阻力系數大幅下降，這也會導致總阻力的增速減緩；當航速繼續增加，興波阻力系數不再大幅下降，故總阻力再次加速上升。

為更加直觀地分析潛艇阻力的變化情況，展示工況A中Fr=0.3,0.6,0.9時潛艇在自由液面處誘發的波紋圖，圖中X、Y坐標代表計算域的尺寸，Z為波幅。由圖5可見，在航速Fr=0.6時，自由液面處的興波要比Fr=0.3,0.9時的波幅大、波紋更加明顯，這與圖4(c)中的興波阻力規律一致，也應證了興波阻力與興波波紋密切相關。隨著航速增加，在自由液面處的興波波長會不斷增加，內波交界面處的興波波形也會出現類似的變化趨勢。

2.1.2 航速對潛艇表面壓強分布的影響

分別選擇A、B 2組工況中Fr=0.3,0.5,0.6,0.9這4個航速工況來展示潛艇表面壓強分布情況。圖6中橫坐標X代表潛艇在航行方向上的位置，單位為m，縱坐標p為壓強，單位是kPa。其中艇艏a點橫坐標為-1.4，艇中b點橫坐標為0.15，艇艉c點橫坐標為1.25。

圖6 不同航速下各工況潛艇表面壓強分布Fig.6 Surface pressure distribution of submarine with different forward speeds under various conditions

由圖6可見，對于A、B 2種工況，潛艇表面壓強分布形狀基本相同。同一航速下，工況A中潛艇表面壓強值要整體低于工況B，這是由于工況A中潛艇所處位置要比工況B中要淺，其表面受到上方流體的靜壓能作用小。

對于同種工況，隨著航速的增加，潛艇艏、艉以及附體等潛艇表面曲率突變處的壓強顯著增大，而其余部位的壓強數值基本保持穩定。如工況A中，航速Fr=0.3時，艇艏a點處壓強p=8.533 4 kPa，當Fr=0.9時a點壓強p=19.249 0 kPa，兩者相差10.715 6 kPa；在艇艉c點處兩者分別為9.646 3 kPa、17.008 8 kPa，相差7.362 5 kPa；而在艇中b點處，兩速度下壓強分別為8.825 0 kPa、8.743 1 kPa，相差僅0.081 9 kPa。該原因是由于航速加大，在潛艇表面凸起引起了動壓的改變，進而使得潛艇表面部分位置處壓力變化。

2.2 位置對潛艇受力的影響

海況中由于密度分層流的存在，導致在淡水、鹽水之間產生了與靜水面對應的內波交界面，本節研究潛艇航行位置與內波交界面距離d1對其阻力特性的影響。采用工況C，使潛艇位于鹽水層中，模擬參數的設置見表3。

表3 模擬參數設置Table 3 Setting of the simulation

2.2.1 潛艇位置對阻力的影響

圖7為距離d1改變時潛艇受力的變化情況。由圖7(a)可見，隨著d1增加，潛艇逐漸遠離內波交界面，其總阻力不斷下降；由圖7(b)可見，由于流體密度及航速保持不變，故摩擦阻力幾乎沒有變化；由圖7(c)可見，興波阻力系數也同樣出現下降的趨勢，這說明當潛艇與內波交界面距離改變時，興波阻力主導了總阻力的變化。此外，隨著潛艇下潛，總阻力及興波阻力系數曲線下降趨勢逐漸變緩，可以預計，當潛艇下潛足夠深時，其興波阻力將接近于零，即在交界面處的興波劇烈程度會越來越弱，進而導致總阻力逐漸穩定于一個較小的值。

圖7 阻力隨d1改變的變化Fig.7 Resistances of submarine with different d3

圖8為潛艇在d1= 0.12L,0.15L,0.18L航行時自由液面處的波紋圖。可見隨著潛艇航行位置不斷下降，其在自由液面上的興波強度逐漸減弱，這解釋了圖7(a)中潛艇所受阻力的下降的情況。

圖8 工況C中不同航速下自由液面興波分布Fig.8 Wave distribution on free surface with different forward speeds under condition C

2.2.2 潛艇位置對潛艇表面壓強分布的影響

圖9為d1=0.12L,0.15L,0.18L時潛艇表面壓強分布。

由圖9可見，位于不同位置處航行的潛艇表面壓強分布形狀基本相同。潛艇上、下表面的壓差幾乎相同，在艇艏、艉及附體存在處壓強值較高。當d1增加時，潛艇表面壓強也逐漸增加，這是由于作用在潛艇單位表面上流體的靜壓能增大導致的。與圖6中航速增加帶來的潛艇表面壓強變化對比可以發現，航速改變會使潛艇表面壓強分布形狀發生變化，即艏、艉等一些拐點處壓強會隨航速增加逐漸變大，但主體部分壓強不變；而航行位置變化只會整體的改變潛艇表面壓強，但其分布形狀不變。通過具體數值來看，當d1=0.12L時，潛艇艏部d點壓強為p=10.887 5 kPa，當d1=0.18L時,壓強為12.709 8 kPa，兩者相差1.822 3 kPa；而對于艇中部e點，2個位置處壓強分別為8.684 5 kPa、10.528 8 kPa，相差1.844 3 kPa。可見，航行位置變化對潛艇表面壓強會帶來整體性的影響。

圖9 不同位置下潛艇表面壓強分布Fig.9 Surface pressure distribution of submarine with different positions

3 結論

1)潛艇航速對其受力特性的影響較大，隨著航速的增加，潛艇遭受阻力不斷增加，總阻力的變化主要來自于興波阻力的改變。

2)與前述研究中近水面航行潛艇對比，在深水層流體中航行的潛艇，其興波阻力系數曲線的峰值會由Fr= 0.5附近向Fr= 0.6附近推移。

3)當潛艇在深水層流體中不同位置處航行時，其受力特征有明顯不同。潛艇距內波交界面越近，其興波愈加劇烈，導致阻力增加。潛艇航行位置改變時，興波阻力主導了總阻力的變化。

4)航速會使潛艇表面壓強分布形狀產生變化，航速增加會使潛艇表面拐點處的壓強變大，但對主體部分壓強影響較小；而位置對潛艇表面壓強的影響是整體性的，由于靜壓能的變化，隨著潛艇遠離內波交界面，其表面壓強會整體增加，但分布形狀不變。

本研究給出了潛艇以不同工況在深水層航行時的受力特征，可為潛艇的航行策略選取提供參考，為密度分層流中潛艇受力研究及艇型設計提供一種可靠的分析手段。