潛艇斜航工況下操舵水動力及繞流場數(shù)值研究

郭海鵬，祁江濤，袁文鑫，李廣年

(1.寧波大學 海運學院，浙江 寧波 315211；2.中國船舶科學研究中心，江蘇 無錫 214082)

0 引 言

操縱性是衡量潛艇綜合性能優(yōu)良的重要指標之一，良好的操縱性能是潛艇安全航行的重要保證[1-2]。潛艇操縱工況下繞流特性極為復雜，不但有主附體結合部和流水孔處各種形式渦的產生與演化，而且還涉及艇體、舵以及螺旋槳之間的復雜相互干擾，導致針對潛艇操縱性的研究具有很高的難度。早期研究主要依靠模型試驗方式進行，對人力物力要求高而且試驗設施復雜，難以在潛艇設計階段進行廣泛應用。

隨著計算機技術的發(fā)展，計算流體動力學（computational fluid dynamics, CFD）方法在水下航行體操縱性預報中發(fā)揮越來越重要的作用。張楠等[3-5]提出了適用于水下航行器流場預報的數(shù)值方法，并開展了帶自由液面、全附體模型的數(shù)值模擬。孫銘澤等[6]實現(xiàn)了全附體潛艇平面運動機構試驗數(shù)值模擬，并將計算得到的水動力導數(shù)與試驗值以及直線拖曳試驗或回轉試驗的數(shù)值計算結果進行了對比驗證。龐永杰[7]以橢球體為研究對象計算了其水動力導數(shù)，并通過與理論值進行對比驗證了所采用數(shù)值方法的有效性。焦玉超等[8]對不同X 形尾舵的潛艇繞流場進行了數(shù)值模擬，分析了X 舵角度對潛艇性能的影響。翟朔等[9]對某潛艇模型的尾部水平操縱面分別進行了共翼型設計和非共翼型設計，并計算了2 種操縱面產生的艇體水動力和尾流特征。

從目前已有研究可以看出，應用CFD 方法開展?jié)撏Р倏v性能預報及評估已經成為領域內的趨勢和研究熱點。本文基于CFD 平臺STAR-CCM+對潛艇斜航工況下的操舵水動力及繞流場進行數(shù)值研究，以國際上廣泛用于對比驗證研究的SUBOFF 潛艇模型為研究對象，分別應用雷諾平均（reynolds-averaged navier-stokes,RANS）與分離渦模擬（detached eddy simulation, DES）2 種方法對潛艇不同舵角下的斜拖試驗進行數(shù)值模擬，據(jù)此分析潛艇縱向力、橫向力、轉首力矩以及繞流場隨舵角的變化規(guī)律，揭示潛艇操縱工況下水動力特性與繞流場間的內在聯(lián)系，為水下航行體大舵角狀態(tài)操縱性能預報和評估提供技術支持。

1 數(shù)值計算方法

采用STAR-CCM+中的RANS 與DES 求解器開展相關數(shù)值計算工作，其中RANS 方法和DES 方法均采用剪切應力傳輸模型（shear stress transfer, SST）k-ω湍流模型進行構造。SSTk-ω兩方程湍流模型為：

式中：ρ為流體密度；uj為速度矢量；xj為位置矢量；k和 ω分別為湍流動能和湍流比耗散率；μ和 μt分別為層流和湍流黏性系數(shù)。源項P和函數(shù)F1的具體形式以及系數(shù) β*，γ，σk，σω和σω2可參考SSTk-ω模型的相關文獻[10]。

將RANS 湍流模型中的長度尺度使用DES 的長度尺度代替，就獲得了基于該湍流模型的混合RANS/LES方法即DES 方法。DES 的長度尺度由RANS 和LES 的長度尺度以如下方式混合得到，即

LES 長度尺度由網格間距得到，即

式中，Δ為局部網格間距，一般取網格單元與相鄰網格單元中心連線長度的最大值。

RANS 長度尺度為最近壁面距離，即

由DES 模型的構造可知，當RANS 長度尺度激活時，湍流模型為傳統(tǒng)的RANS 湍流模型，而當LES 長度尺度激活時，湍流模型轉換為經典的Smagorisnky 亞格子模型，從而實現(xiàn)RANS 和LES 方法的轉換。

2 數(shù)值計算設置

2.1 研究對象

以美國DARPA 潛艇模型SUBOFF-8 為研究對象，該潛艇由裸艇、指揮臺圍殼及4 個完全相同的尾翼（左右為升降舵，上下為方向舵）組成，主艇體總長Loa=4.356 m，垂線間長L=4.261 m，其中前體長1.016 m，平行中體長2.229 m，后體長1.111 m；艇體最大直徑D=0.508 m；指揮臺圍殼長0.368 m，高0.205 m，上部有2∶1 的橢圓形橫截面的頂蓋；4 個相同的尾翼截面為NACA0020，呈十字布置在尾部，頂端弦長為0.152 m，詳細幾何特征由文獻[11]提供，三維幾何形狀如圖1所示。

圖1 SUBOFF 幾何模型Fig.1 Geometric model of SUBOFF

2.2 計算工況

計算工況具體參數(shù)如表1 所示。考慮潛艇在漂角β=-4°、舵角0°～40°范圍內的水動力性能，其中δ在0°～15°范圍內數(shù)值計算結果可與Roddy[12]所做的公開實驗數(shù)據(jù)進行對比，從而確定數(shù)值方法的可靠性。在STAR-CCM+中，漂角的處理方法是將SUBOFF 以艇體質心所在的Z軸為旋轉軸，進行逆時針旋轉4°，其他保持不變。

表1 工況設定Tab.1 Setup of study condition

2.3 計算域設置

計算域總長為7 倍艇長，首部距離進口1.5 倍，尾部距離出口4.5 倍，左右及上下邊界距艇體均為2.0 倍艇長。對于邊界條件的設定，除計算域出流面采用壓力出口條件以外，其他邊界均采用速度入口條件，艇體所有表面采用無滑移壁面邊界條件。計算域尺度及邊界條件如圖2 所示。

圖2 計算域尺度及邊界條件Fig.2 Dimensions and boundary conditions of the computational domain

計算域采用STAR-CCM+中的切割體網格進行離散。切割體網格屬于非結構網格，適用于具有復雜外形模型的網格劃分。切割體網格在生成時可以通過控制體進行任意區(qū)域的網格加密，同時還可以對幾何模型中的控制線、控制面進行線加密及面加密。為了保證潛艇附近的網格質量，對其表面進行了加密；通過使用切割體網格單元，在艇體附近及其尾流進行網格密化，從而更好地捕捉流動的細節(jié)，獲得更準確的流場信息。同時，在近壁面生成棱柱層網格以便于邊界層流動的求解。通過調整棱柱層網格的總厚度、層數(shù)及增長比率來控制近壁面第一層網格高度，使近壁面Y+值在30～60 范圍內。網格數(shù)量在500 萬左右，網格劃分情況如圖3 所示。

圖3 計算域網格劃分情況Fig.3 Gird generation in the computational domain

3 結果分析與討論

3.1 水動力及力矩

圖4 給出了RANS 與DES 方法獲得的作用在潛艇上的縱向力X、橫向力Y及轉首力矩N。為了便于分析，均采用無因次形式表述：

圖4 數(shù)值結果與模型試驗數(shù)據(jù)對比Fig.4 Comparison between the numerical results and the experimental data

式中：ρ為流體密度；V為航行速度；L為潛艇總長。

相對X′而言，無論RANS 還是DES 對阻力的預報精度都在±3%以內；δ的變化會對阻力產生較大的影響，δ在0°～30°范圍內，隨著舵角的增大，艇體所受的阻力以拋物線的趨勢增大；隨著δ的進一步增大，RANS 計算結果在30°～40°范圍內基本呈線性增大，DES 計算結果先有略微減小然后再增大。相對Y′而言，RANS 計算結果與試驗值的吻合程度略好于DES 方法，RANS 計算的誤差均在10%以內，DES 計算除了δ=15°時誤差偏大，其余工況均在12% 以內；δ在0°～30°范圍內，DES 與RANS 計算得到的橫向力隨著舵角的增大而近似呈現(xiàn)線性增大；δ在30°～35°范圍內，RANS 及DES 計算的橫向力顯著下降，發(fā)生了明顯的舵失速現(xiàn)象；進一步增大舵角，艇體所受的橫向力有微幅的增加，這表明舵效并未隨舵角進一步增大而發(fā)生明顯改善。相對于N′而言，DES 計算在各工況下的計算精度均在±12% 以內，而RANS 方法在δ=10°和δ=15°計算誤差較大，均超過20%；N′與Y′的變化規(guī)律保持基本一致。通過對水動力分析可知，就對比現(xiàn)有的公開試驗數(shù)據(jù)而言，RANS 方法在X′與Y′相關的計算精度比DES 方法略好，但是DES 方法對N′的計算精度明顯優(yōu)于RANS 方法。

3.2 泄出渦系結構

圖5 給出了由RANS 及DES 方法獲得的不同舵角下艇尾區(qū)域的渦系結構，其中渦系結構采用Q準則等值面表示，并以無量綱縱向速度Ux/U0著色。可以看出，隨著方向舵舵角的變化，艇尾區(qū)域渦系發(fā)展呈現(xiàn)出顯著的差異。在小舵角（δ=10°）下，舵背風面主要以附著流為主，舵梢部有輕微渦泄出。隨著舵角增大（δ=20°），方向舵附近產生的渦系結構明顯增多。方向舵背風面靠近尾部區(qū)域開始產生泄出渦系。同時，由于偏轉方向舵的阻滯效應，在艇身與方向舵的交匯處形成了明顯的馬蹄渦結構。隨著舵角的進一步增大（δ=30°），舵背風面完全由泄出渦系覆蓋，而且馬蹄渦結構也更為顯著。當舵角達到δ=40°時，舵背風面的泄出渦完全脫離舵表面發(fā)展，在下游形成了復雜的渦系結構，與舵失速狀態(tài)相對應。通過對比可知，在小舵角工況下，RANS 方法與DES 方法得到的艇尾區(qū)域渦系結構差異不明顯，隨著舵角增大，二者的差異越發(fā)明顯。特別是發(fā)生大規(guī)模流動分離的工況下，DES方法捕捉到的渦系結構更加精細。

圖5 不同舵角工況下的Q=100 s-2 等值面圖Fig.5 Iso-surface of Q=100 s-2 under different rudder angles

3.3 流線分布

圖6 給出了DES 與RANS 方法獲得的舵葉中間位置處的舵附近流線分布。可知，在舵角δ=10°工況下，舵剖面附近流動主要以附著流為主，而且2 種方法得到的流場沒有明顯的差異。在舵角δ=20°工況下，舵剖面背流面發(fā)生了明顯的流動分離現(xiàn)象，其中RANS方法呈現(xiàn)出的流動分離更為明顯，導致該區(qū)域的速度明顯下降。在舵角δ=30°工況下，舵剖面的整個背流面均發(fā)生了流動分離，2 種數(shù)值方法獲得的流線分布差異也更加明顯。RANS 方法中背流面處的流動分離更加明顯，而DES 方法捕捉到的流動分離并未隨舵角明顯增強。另外，2 種方法捕捉得到流場均在舵下游呈現(xiàn)出明顯的低速區(qū)。在舵角δ=40°工況下，舵剖面背流面處的流動分離顯著增強，對附近流場產生的擾動也十分明顯，在下游形成了更為顯著的低速區(qū)，這一流場特征與舵失速工況下的阻力增加及升力下降相對應。綜上，DES 方法在流場捕捉方面具有明顯優(yōu)勢，能清晰地捕捉艇尾區(qū)域的復雜流場信息。

圖6 不同舵角工況下舵附近流線分布Fig.6 Streamlines around the rudder under different rudder angles

3.4 壓力分布

圖7 給出了RANS 方法與DES 方法獲得的不同舵角工況下的舵剖面壓力分布情況。壓力采用無因次系數(shù)Cp表示，Cp=p/(0.5ρV2)。橫坐標由無因次量x/C表示，坐標原點位于該舵前緣處，無因次量“1”表示舵尾緣。在舵角δ=10°工況下，舵剖面壓力呈現(xiàn)出明顯的吸力面和壓力面，正壓和負壓的峰值均位于前緣附近，而且2 種方法獲得的壓力分布較為一致。在舵角δ=20°工況下，DES 方法得到的吸力面負壓區(qū)較RANS方法更為顯著，但整體分布差異很小。在舵角δ=30°工況下，吸力面負壓區(qū)更為顯著，同時發(fā)生流動分離的區(qū)域壓力分布較為平緩。DES 方法得到的吸力面負壓區(qū)明顯強于RANS 方法，這也與該舵角工況下DES計算得到的舵橫向力高于RANS 計算結果相對應。在舵角δ=40°工況下，吸力面負壓區(qū)峰值明顯下降，而且發(fā)生流動分離的區(qū)域壓力分布呈平直狀，導致吸力面與壓力面之間的壓力差顯著下降，這與舵發(fā)生失速導致的舵效下降相對應。總體而言，2 種數(shù)值方法得到的壓力分布差異相對較小，與潛艇操舵水動力性能的一般規(guī)律較為吻合。

圖7 不同舵角工況下舵剖面壓力分布圖Fig.7 Pressure distribution on the rudder section under different rudder angles

4 結 語

本文針對潛艇斜航工況下的操舵水動力性能及繞流場進行數(shù)值模擬研究，主要得到以下結論：

1）潛艇斜航工況下的縱向力、橫向力及轉首力矩隨著舵角的變化呈現(xiàn)復雜的變化規(guī)律，在潛艇操舵水動力預報精度方面DES 方法整體上優(yōu)于RANS 方法；

2）DES 方法對于潛艇操舵工況下艇尾區(qū)域的繞流場捕捉更為精細，在潛艇操縱工況下的水動力學機理探討方面有明顯的優(yōu)勢，在潛艇復雜操縱運動分析和性能評估方面具有很好的工程應用潛力。