基于CFD的潛艇模型自航仿真分析

楊 琴，王國棟，張志國，馮大奎，王先洲

華中科技大學船舶與海洋工程學院，湖北武漢430074

0 引 言

對潛艇周圍的流場信息特別是尾部流動特性進行研究一直是潛艇設計和流體力學領域的一項重要課題。當潛艇尾部螺旋槳盤面處的伴流場不均勻和不穩定時，會降低螺旋槳的推進效率。研究潛艇螺旋槳與尾部操縱面之間的互相影響，對提高潛艇的推進效率和機動性具有重要意義?？赏ㄟ^潛艇的自航試驗來獲得螺旋槳與船體之間的相互作用等諸因素，這是分析研究潛艇快速性與推進效率的重要手段。

數值水池的重要基礎是CFD 技術，即應用CFD 方法進行船舶流動數值模擬來實現船舶水動力學性能的數值計算與預報。數值水池船舶阻力與螺旋槳敞水性能計算發展較早。對全附體潛艇的數值模擬一般是采用RANS 方程結合k-ε 湍流模型求解［1-2］，或采用大渦模擬方法。目前，螺旋槳的敞水性能研究已發展得較為成熟，一般采用混合網格［3］來對大側斜螺旋槳進行數值計算，常用的湍流模型為k-ε，SST 和雷諾應力模型［4-5］，大渦湍流模型［6］也被運用到了大側斜螺旋槳的數值計算中，以上數值計算結果與實驗室值都已達到了較高的精度［7-8］。數值水池自航試驗研究需實現螺旋槳與船體的整體求解，Choi 等［9］分別采用體力槳和滑移網格方法對船—槳—舵之間的水動力相互影響進行了計算，并與模型試驗值進行了對比，獲得了比較準確的計算結果；劉祥珺等［10］通過求解RANS 方程，利用VOF 方法追蹤自由液面，在數值水池實現了船模的自航試驗研究。螺旋槳旋轉域與靜止域采用滑移網格技術。在既定航速下，船體阻力和螺旋槳推力為螺旋槳轉速的函數，通過變化轉速，可以得到自航點。采用數值方法得到的自航轉速與試驗值吻合良好，張楠等［11］采用數值模擬方法研究了潛艇近水面航行時的艇/槳干擾特性。流場采用RANS 方法結合RNG k-ε 湍流模型求解；自由液面捕捉采用VOF方法；螺旋槳采用滑移網格方法。

1 理論基礎

本文采用RNG k-ε 湍流模型求解Reynolds平均Navier-Stokes 方程，建立了相關的仿真數學模型，對帶七葉螺旋槳的SUBOFF 全附體模型的三維粘性流場進行了數值分析，并采用準靜態的方法對潛艇模型的阻力和螺旋槳推力進行了耦合，得到了指定速度所對應的螺旋槳推力和轉速。

1.1 控制方程

本文研究分析的“艇—槳”自航條件下螺旋槳的控制方程為不可壓縮牛頓流體流動的連續性方程和RANS 方程：

1.2 湍流模型

在渦粘模型中，根據Boussinesq 提出的渦粘假定，引入了湍動粘度，該假定建立了雷諾應力與平均速度梯度的關系，即

通常，渦粘模型有零方程模型、一方程模型和二方程模型。本文采用的RNG k-ε 湍流模型對應的方程組如下：

1.3 敞水特征曲線

螺旋槳模型單獨在均勻水流中試驗稱為敞水試驗，由螺旋槳敞水試驗得到的是螺旋槳的推力系數KT、扭矩系數KQ和敞水效率η0相對于進速系數J 的變化規律，即螺旋槳敞水特征曲線。與敞水試驗對應的船舶計算流體力學（CFD）計算工作稱之為螺旋槳敞水性能計算。采用CFD 方法計算螺旋槳的水動力性能，所采用的計算模型具有較高的可靠性，其結果也具有一定的準確性。

螺旋槳的敞水性能包括螺旋槳的推力系數KT、扭矩系數KQ和敞水效率η0，具體計算公式如下：

式中，

式中：J 為進速系數，是影響螺旋槳性能的重要參數，其相當于機翼理論中攻角的概念；Va代表螺旋槳的進速；T 為螺旋槳產生的軸向推力；Q 為螺旋槳運行中產生的扭矩；ρ 為工作介質（20 ℃水）的密度；D 為螺旋槳直徑；n 為螺旋槳轉速。

1.4 螺旋槳與艇的相互作用

伴流的大小通常用伴流速度u 與船速V 的比值ω 來表示。ω 稱為伴流分數，即

螺旋槳在船后工作時引起的船體附加阻力稱為阻力增額，螺旋槳發出的推力中只有一部分是用于克服阻力R 并推船前進的，稱這一部分力為有效推力T 。通常，將阻力增額稱為推力減額，并以ΔT 表示。推力減額ΔT 與推力T 的比值稱為推力減額分數：

2 計算對象

本文的計算對象為帶槳SUBOFF 全附體模型。該模型艇長4.356 1 m，螺旋槳采用改進的IN?SEANE 1619 型七葉螺旋槳。為配合全附體潛艇的尾端直徑，建立了直徑為302.6 mm 的七葉螺旋槳三維模型，其主要幾何參數如表1 所示。

表1 螺旋槳幾何參數Tab.1 Geometric parameters of propeller

本文采用SolidWorks 軟件對計算對象進行了三維建模，所建立的帶槳SUBOFF 全附體三維模型示意圖如圖1 所示。圖中，x 軸與螺旋槳旋轉軸一致，以指向船尾為正。

圖1 帶槳SUBOFF 三維模型Fig.1 3D model of SUBOFF with propeller

3 數值模擬方法

3.1 計算域與網格劃分

為了模擬帶槳艇周圍的流場，進行數值計算時要確定一個大小合適的流場控制體。本文將帶槳艇模放在了一個與螺旋槳同軸線的圓柱形流場區域中，進流面取為上游3 m，出流面取為下游5 m，徑向取為2.5 m。其中，螺旋槳被一圓柱小體包裹，此為旋轉區域，長0.18 m，直徑0.5 m，如圖2所示。

計算采用結構化網格與非結構化網格相結合的方式，在包裹螺旋槳的小體內采用非結構化網格，其余部分采用結構化網格。進行了合理的布置并選擇了適當的網格尺度，以期取得很好的計算精度和計算時間。圖3 給出了帶槳潛艇表面的網格分布，螺旋槳附近非結構化網格數均為50 萬左右，整個計算域的網格數約為250 萬。

圖2 數值模擬計算域Fig.2 Computation domain of numerical simulation

圖3 帶槳潛艇表面的網格劃分Fig.3 Surface mesh on the SUBOFF simulation model

3.2 邊界條件

本文對SUBOFF主艇體的數值模擬邊界如下：

入口條件（Velocity Inlet）：入口及計算域周向設置為速度入口，速度大小為Va= 8 kn，速度方向平行于槳軸，為x 軸正方向。

出口條件（Pressure Outlet）：認為流動在該處已充分發展，故邊界條件設置為壓力出口。

壁面條件（Wall）：螺旋槳表面和艇體外均設置為無滑移壁面條件。

Interface：由于使用滑移網格模擬螺旋槳的旋轉運動，故螺旋槳旋轉區域與外流場區域交界面均設置為interface。

坐標系：坐標原點在槳盤面圓心處，x 軸沿來流正向布置，坐標系符合右手定則。

3.3 自航點的確定

在已知不同航速下艇體阻力曲線和螺旋槳敞水性能特征的基礎上，對“艇—槳”整體進行數值模擬的自航試驗。自航點的尋找流程如圖4所示。

具體步驟如下：

1）由阻力曲線查得給定航速Va= 8 kn 下的艇體阻力Ds。

2）依據經驗給定初始推力減額系數t0=0.079，由全附體阻力Ds和t0求出螺旋槳初始推力T0，然后再由給定螺旋槳敞水性能曲線查得效率最高點對應的推力系數KT，最后，由T0和KT求得初始轉速n0。

圖4 確定自航點流程示意圖Fig.4 Flow chart of finding out self-propulsion point

3）給定Va和n0后，對“艇—槳”整體進行RANS 模擬。待迭代收斂后，讀取艇體阻力D 和螺旋槳推力T 值。當T ＞D 時，減小轉速n ；反之，則增加轉速n。改變轉速后，再次以初始計算流場為初值進行RANS 模擬，并重復該過程進行迭代計算，直至推力與阻力的平均值相等或達到允許的不平衡度為止。

4）找到自航點后，求取螺旋槳旋轉域進流面處的有效伴流系數ωe、螺旋槳推力減額系數t 和相對旋轉效率ηR。

4 數值計算結果與分析

4.1 艇帶槳水動力性能預報

對螺旋槳的定常性能模擬采用FLUENT 提供的滑移網格模型。采用RNG k-ε 湍流模型對SUBOFF 全附體模型進行阻力計算。在以上基礎上，對帶槳潛艇在來流速度Va= 8 kn 下進行數值模型的自航試驗，取螺旋槳轉速n = 13.2 為初始轉速，并預報其水動力性能。

本文以DTMB 4119 等螺旋槳為對象進行了研究。采用本文的數值模擬方法計算所得的結果與已知實驗值［12］的對比如表2 所示。SUBOFF 全附體模型阻力的數值模擬結果與試驗數據的比較如圖5 所示。以上誤差均控制在3%以內，充分證明了本文采用的數值模擬方法具有較高的準確性。

圖5 全附體Suboff實驗值與數值計算結果對比Fig.5 Comparison between experiment and numerical simulation of SUBOFF

在以上基礎上，對SUBOFF 全附體螺旋槳模型的計算網格，將其Y+值控制在了3.03～41.92 范圍內。圖6 給出了螺旋槳壓力面與吸力面上的壓力分布云圖。從圖中可以看出，7 個槳葉的壓力分布基本相同，且壓力面從隨邊至導邊，從葉根至葉梢，壓力是逐漸增加的，葉根處的壓力最小；在吸力面，壓力是中間低四周高，靠近葉梢的部分存在低壓區。對于一個葉切面而言，壓差最大點位于切面最大厚度處。由圖7 給出的帶槳潛艇表面壓力分布圖可看出，在艇前、圍殼前和機翼前均有較大的壓力梯度。由尾流處的軸向速度分布（圖8）可以看到，螺旋槳附近有較大的速度變化。而由于螺旋槳的旋轉運動，槳后尾流速度場呈繞軸線螺旋狀分布，其尾流流線如圖9 所示。

圖6 螺旋槳吸力面和壓力面壓力分布云圖Fig.6 Pressure distribution on the blade surface of suction side and pressure side

圖7 帶槳潛艇表面壓力分布圖Fig.7 Pressure distribution on the surface of SUBOFF

圖8 尾流處軸向速度分布圖Fig.8 Axial velocity distribution

圖9 帶槳艇流線圖Fig.9 Velocity streamlines

4.2 數值自航計算與結果分析

由給定的潛艇阻力曲線，得到對應的來流速度Va= 8 kn 下SUBOFF 全附體阻力Ds= 185 N。依據潛艇推力減額分數設計參考值（0.1～0.18）［13］，設定初始推力減額分數t = 0.14。在給定螺旋槳的敞水性能曲線上，查得最大效率點取KT值，得到初始轉速n = 0.221 7 r/min。根據自航點確定流程，得到帶槳的潛艇阻力與螺旋槳推力之間的關系如圖10 所示。

圖10 數值自航結果分析Fig.10 The analysis of numerical self-propulsion result

由圖10 可知，在推力與阻力相等時，即為潛艇自航點，螺旋槳轉速n=0.220 7 r/min(J=0.998)，潛艇阻力與螺旋槳的推力值均為213 N。求得螺旋槳旋轉域進流面處的有效伴流系數ωe= 0.17，螺旋槳推力減額系數t = 0.13，相對旋轉效率ηR=91%，根據數值模擬得到的上述結果，與常規設計手冊提供的參數選擇范圍進行比較，均在合理的取值范圍內。

5 結 語

本文采用數值模擬方法系統地研究了全附體潛艇+螺旋槳的水動力特性，清晰、形象地描述了帶槳潛艇表面的壓力分布情況，以及槳后尾流速度場的繞軸線螺旋狀分布規律。對螺旋槳在艇后非均勻進流條件下的推力和力矩的脈動特性進行了分析，得到了可靠的帶槳潛艇阻力和螺旋槳推力數據，潛艇自航點在螺旋槳轉速n = 0.220 7 r/min（J = 0.998）時，潛艇阻力與螺旋槳的推力值均為213 N。根據潛艇數值自航試驗結果，可以通過積分的方法準確計算伴流分數與推力減額分數，詳細分析螺旋槳葉片表面的壓力分布，梢渦特性和壓力脈動特性。該方法的開發可以顯著減少試驗研究的工作量，縮短研究周期，具有廣闊的應用前景。

［1］吳方良，吳曉光，馬運義，等.基于雷諾應力方程模型的全附體潛艇尾部伴流場三維粘性數值計算［J］.中國艦船研究，2007，2（6）：4-8.WU Fangliang，WU Xiaoguang，MA Yunyi，et al. Nu?merical calculation of 3D viscous wake field of subma?rine with full appendages based on Reynolds stress equation model［J］. Chinese Journal of Ship Research，2007，2（6）：4-8.

［2］吳方良，吳曉光，馬運義，等.潛艇實艇阻力預報方法研究［J］.中國艦船研究，2009，4（3）：28-32.WU Fangliang，WU Xiaoguang，MA Yunyi，et al. Meth?od of predicting the total submerged resistance of sub?marines［J］. Chinese Journal of Journal of Ship Re?search，2009，4（3）：28-32.

［3］RHEE S H，LEE C，LEE H B，et al. Rudder gap cavita?tion：Fundamental understanding and its suppression devices［J］. International Journal of Heat and Fluid Flow，2010，31（4）：640-650.

［4］楊瓊方，王永生，張志宏，等.大側斜螺旋槳敞水性能的RANSEs 模擬［J］. 湖南科技大學學報（自然科學版），2008，23（4）：66-70.YANG Qiongfang，WANG Yongsheng，ZHANG Zhi?hong，et al. RANSEs simulation of highly skewed pro?peller’s open water performances［J］. Journal of Hu?nan University of Science and Technology：Natural Sci?ence Edition，2008，23（4）：66-70.

［5］趙鵬飛. 船用螺旋槳敞水和空化性能CFD 預報［D］.大連：大連理工大學，2011.ZHAO Pengfei. CFD predictions of open water and cav?itation performances of marine propellers［D］. Dalian ：Dalian University of Technology，2011.

［6］FELICE F D，FELLI M，LIEFVENDAHL M，et al.Numerical and experimental analysis of the wake be?havior of a generic submarine propeller［C］//First Inter?national Symposium on Marine Propulsors（SMP'09）.Trondheim，Norway，2009.

［7］ANDERSEN P，KAPPEL J J，SPANGENBERG E.As?pects of Propeller developments for a submarine［C］//First International Symposium on Marine Propulsors（SMP'09）.Trondheim，Norway，2009.

［8］胡磊，譚廷壽.大側斜螺旋槳敞水性能分析［J］.船海工程，2010，39（1）：41-44.HU Lei，TAN Tingshou. Analysis of open water perfor?mance of the high-skew propeller［J］. Ship and Ocean Engineering，2010，39（1）：41-44.

［9］CHOI J E，MIN K S，KIM J H，et al. Resistance and propulsion characteristics of various commercial ships based on CFD results［J］. Ocean Engineering，2010，37（7）：549-566.

［10］劉祥珺，孫存樓. 數值水池船模自航試驗方法研究［J］.艦船科學技術，2011，33（2）：28-31.LIU Xiangjun，SUN Cunlou. Research on ship self-propulsion model test in numerical tank［J］. Ship Science and Technology，2011，33（2）：28-31.

［11］張楠，張勝利，沈泓萃，等.帶自由液面的艇/槳干擾特性數值模擬與驗證研究［J］. 水動力學研究與進展A 輯，2012，27（1）：94-99.ZHANG Nan，ZHANG Shengli，SHEN Hongcui，et al.Numerical simulation of hull/propeller interaction with free surface［J］.Chinese Journal of Hydrodynam?ics，2012，27（1）：94-99.

［12］盛振邦，劉應中. 船舶原理［M］. 上海：上海交通大學出版社，2004.

［13］BOSWELL R J. Design，cavitation performance，and open water performance of a series of research skewed propellers. NSRDC Report No. 3339 ［R］. 1971：4381-4383.