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

楊 琴，王國(guó)棟，張志國(guó)，馮大奎，王先洲

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

0 引 言

對(duì)潛艇周圍的流場(chǎng)信息特別是尾部流動(dòng)特性進(jìn)行研究一直是潛艇設(shè)計(jì)和流體力學(xué)領(lǐng)域的一項(xiàng)重要課題。當(dāng)潛艇尾部螺旋槳盤面處的伴流場(chǎng)不均勻和不穩(wěn)定時(shí)，會(huì)降低螺旋槳的推進(jìn)效率。研究潛艇螺旋槳與尾部操縱面之間的互相影響，對(duì)提高潛艇的推進(jìn)效率和機(jī)動(dòng)性具有重要意義。可通過潛艇的自航試驗(yàn)來獲得螺旋槳與船體之間的相互作用等諸因素，這是分析研究潛艇快速性與推進(jìn)效率的重要手段。

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

1 理論基礎(chǔ)

本文采用RNG k-ε 湍流模型求解Reynolds平均Navier-Stokes 方程，建立了相關(guān)的仿真數(shù)學(xué)模型，對(duì)帶七葉螺旋槳的SUBOFF 全附體模型的三維粘性流場(chǎng)進(jìn)行了數(shù)值分析，并采用準(zhǔn)靜態(tài)的方法對(duì)潛艇模型的阻力和螺旋槳推力進(jìn)行了耦合，得到了指定速度所對(duì)應(yīng)的螺旋槳推力和轉(zhuǎn)速。

1.1 控制方程

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

1.2 湍流模型

在渦粘模型中，根據(jù)Boussinesq 提出的渦粘假定，引入了湍動(dòng)粘度，該假定建立了雷諾應(yīng)力與平均速度梯度的關(guān)系，即

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

1.3 敞水特征曲線

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

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

式中，

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

1.4 螺旋槳與艇的相互作用

伴流的大小通常用伴流速度u 與船速V 的比值ω 來表示。ω 稱為伴流分?jǐn)?shù)，即

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

2 計(jì)算對(duì)象

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

表1 螺旋槳幾何參數(shù)Tab.1 Geometric parameters of propeller

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

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

3 數(shù)值模擬方法

3.1 計(jì)算域與網(wǎng)格劃分

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

計(jì)算采用結(jié)構(gòu)化網(wǎng)格與非結(jié)構(gòu)化網(wǎng)格相結(jié)合的方式，在包裹螺旋槳的小體內(nèi)采用非結(jié)構(gòu)化網(wǎng)格，其余部分采用結(jié)構(gòu)化網(wǎng)格。進(jìn)行了合理的布置并選擇了適當(dāng)?shù)木W(wǎng)格尺度，以期取得很好的計(jì)算精度和計(jì)算時(shí)間。圖3 給出了帶槳潛艇表面的網(wǎng)格分布，螺旋槳附近非結(jié)構(gòu)化網(wǎng)格數(shù)均為50 萬左右，整個(gè)計(jì)算域的網(wǎng)格數(shù)約為250 萬。

圖2 數(shù)值模擬計(jì)算域Fig.2 Computation domain of numerical simulation

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

3.2 邊界條件

本文對(duì)SUBOFF主艇體的數(shù)值模擬邊界如下：

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

出口條件（Pressure Outlet）：認(rèn)為流動(dòng)在該處已充分發(fā)展，故邊界條件設(shè)置為壓力出口。

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

Interface：由于使用滑移網(wǎng)格模擬螺旋槳的旋轉(zhuǎn)運(yùn)動(dòng)，故螺旋槳旋轉(zhuǎn)區(qū)域與外流場(chǎng)區(qū)域交界面均設(shè)置為interface。

坐標(biāo)系：坐標(biāo)原點(diǎn)在槳盤面圓心處，x 軸沿來流正向布置，坐標(biāo)系符合右手定則。

3.3 自航點(diǎn)的確定

在已知不同航速下艇體阻力曲線和螺旋槳敞水性能特征的基礎(chǔ)上，對(duì)“艇—槳”整體進(jìn)行數(shù)值模擬的自航試驗(yàn)。自航點(diǎn)的尋找流程如圖4所示。

具體步驟如下：

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

2）依據(jù)經(jīng)驗(yàn)給定初始推力減額系數(shù)t0=0.079，由全附體阻力Ds和t0求出螺旋槳初始推力T0，然后再由給定螺旋槳敞水性能曲線查得效率最高點(diǎn)對(duì)應(yīng)的推力系數(shù)KT，最后，由T0和KT求得初始轉(zhuǎn)速n0。

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

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

4）找到自航點(diǎn)后，求取螺旋槳旋轉(zhuǎn)域進(jìn)流面處的有效伴流系數(shù)ωe、螺旋槳推力減額系數(shù)t 和相對(duì)旋轉(zhuǎn)效率ηR。

4 數(shù)值計(jì)算結(jié)果與分析

4.1 艇帶槳水動(dòng)力性能預(yù)報(bào)

對(duì)螺旋槳的定常性能模擬采用FLUENT 提供的滑移網(wǎng)格模型。采用RNG k-ε 湍流模型對(duì)SUBOFF 全附體模型進(jìn)行阻力計(jì)算。在以上基礎(chǔ)上，對(duì)帶槳潛艇在來流速度Va= 8 kn 下進(jìn)行數(shù)值模型的自航試驗(yàn)，取螺旋槳轉(zhuǎn)速n = 13.2 為初始轉(zhuǎn)速，并預(yù)報(bào)其水動(dòng)力性能。

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

圖5 全附體Suboff實(shí)驗(yàn)值與數(shù)值計(jì)算結(jié)果對(duì)比Fig.5 Comparison between experiment and numerical simulation of SUBOFF

在以上基礎(chǔ)上，對(duì)SUBOFF 全附體螺旋槳模型的計(jì)算網(wǎng)格，將其Y+值控制在了3.03～41.92 范圍內(nèi)。圖6 給出了螺旋槳壓力面與吸力面上的壓力分布云圖。從圖中可以看出，7 個(gè)槳葉的壓力分布基本相同，且壓力面從隨邊至導(dǎo)邊，從葉根至葉梢，壓力是逐漸增加的，葉根處的壓力最?。辉谖γ?，壓力是中間低四周高，靠近葉梢的部分存在低壓區(qū)。對(duì)于一個(gè)葉切面而言，壓差最大點(diǎn)位于切面最大厚度處。由圖7 給出的帶槳潛艇表面壓力分布圖可看出，在艇前、圍殼前和機(jī)翼前均有較大的壓力梯度。由尾流處的軸向速度分布（圖8）可以看到，螺旋槳附近有較大的速度變化。而由于螺旋槳的旋轉(zhuǎn)運(yùn)動(dòng)，槳后尾流速度場(chǎng)呈繞軸線螺旋狀分布，其尾流流線如圖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 數(shù)值自航計(jì)算與結(jié)果分析

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

圖10 數(shù)值自航結(jié)果分析Fig.10 The analysis of numerical self-propulsion result

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

5 結(jié) 語

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

［1］吳方良，吳曉光，馬運(yùn)義，等.基于雷諾應(yīng)力方程模型的全附體潛艇尾部伴流場(chǎng)三維粘性數(shù)值計(jì)算［J］.中國(guó)艦船研究，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］吳方良，吳曉光，馬運(yùn)義，等.潛艇實(shí)艇阻力預(yù)報(bào)方法研究［J］.中國(guó)艦船研究，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］楊瓊方，王永生，張志宏，等.大側(cè)斜螺旋槳敞水性能的RANSEs 模擬［J］. 湖南科技大學(xué)學(xué)報(bào)（自然科學(xué)版），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 預(yù)報(bào)［D］.大連：大連理工大學(xué)，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，F(xiàn)ELLI 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］胡磊，譚廷壽.大側(cè)斜螺旋槳敞水性能分析［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］劉祥珺，孫存樓. 數(shù)值水池船模自航試驗(yàn)方法研究［J］.艦船科學(xué)技術(shù)，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］張楠，張勝利，沈泓萃，等.帶自由液面的艇/槳干擾特性數(shù)值模擬與驗(yàn)證研究［J］. 水動(dòng)力學(xué)研究與進(jìn)展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］盛振邦，劉應(yīng)中. 船舶原理［M］. 上海：上海交通大學(xué)出版社，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.