基于兩種螺旋槳建模方法的全附體船模斜拖試驗數(shù)值模擬

劉 義, 鄒早建, 郭海鵬

(1. 中國船舶及海洋工程設(shè)計研究院， 上海 200011；2. 上海交通大學 船舶海洋與建筑工程學院， 上海 200240；3. 上海交通大學 海洋工程國家重點實驗室， 上海 200240)

船舶操縱性是與船舶航行安全性密切相關(guān)的非常重要的船舶航行性能，準確預報船舶操縱性對保證船舶航行安全具有重要意義.在船舶設(shè)計階段，準確預報船舶操縱性是至關(guān)重要的.十多年來，計算流體動力學(CFD)方法在船舶水動力學學科得到日益廣泛的開發(fā)和應(yīng)用，其發(fā)展異常迅猛.在船舶操縱性方面，CFD方法也已得到成功的應(yīng)用，包括數(shù)值模擬約束模試驗以獲取操縱運動數(shù)學模型中的水動力，以及對操縱運動進行直接數(shù)值模擬以實現(xiàn)對操縱性的直接預報.這方面的研究已成為了國際上的熱點和前沿課題[1].

數(shù)學模型加計算機模擬的方法是一種在船舶設(shè)計階段預報操縱性的最常用的方法，而基于CFD方法對約束模試驗進行數(shù)值模擬是現(xiàn)階段獲取船舶操縱運動數(shù)學模型中的水動力導數(shù)的一種有效方法.Yao等[2]基于OpenFOAM開發(fā)RANS求解器，對大型油輪KVLCC2船模的斜拖試驗、旋臂試驗和PMM(平面運動機構(gòu))純橫蕩試驗進行了數(shù)值模擬.王建華等[3]基于非定常RANS方程和重疊網(wǎng)格技術(shù)，對軍船DTMB5415裸船體模型的PMM純首搖試驗進行了數(shù)值模擬.馮松波等[4]應(yīng)用CFD軟件Fluent對KVLCC2船-舵系統(tǒng)模型的斜拖試驗進行了數(shù)值模擬.當采用CFD方法數(shù)值模擬全附體約束模試驗時，重點和難點之一在于螺旋槳的建模.目前應(yīng)用CFD方法對真實幾何形狀的螺旋槳(實槳)和船體相互干擾進行全附體船舶操縱運動數(shù)值模擬已經(jīng)開展了一些研究.Badoe等[5]采用實槳建模方法，分析了斜流下的船-槳-舵干擾問題；類似地，Wang等[6]數(shù)值模擬船舶的斜拖試驗，分析了螺旋槳在斜流中的負荷.除實槳建模外，另一種方法是以力場模擬代替真實螺旋槳，即體積力法.國內(nèi)外已有眾多學者應(yīng)用體積力法研究過船-槳-舵干擾和操縱性等問題[7-9].但在約束模試驗數(shù)值模擬中，基于實槳建模和基于體積力建模之間的差別卻少有涉及.

本文應(yīng)用STAR CCM+軟件，采用雷諾平均Navier-Stokes(RANS)方程求解方法并選取Realizablek-ε湍流模型封閉黏性流體流動控制方程，對全附體的集裝箱船KCS船模斜拖試驗的黏性流場進行了數(shù)值模擬，其中對螺旋槳建模分別采用了體積力模型和實體螺旋槳模型.通過將數(shù)值計算結(jié)果與基準試驗數(shù)據(jù)進行對比，驗證了數(shù)值方法的正確性；同時，對比分析了兩種螺旋槳建模方法在計算效率和計算精度方面的差異.

1 數(shù)學模型

1.1 坐標系

如圖1所示，建立兩個右手坐標系：空間固定坐標系O0-x0y0z0和隨船運動坐標系O-xyz.空間固定坐標系原點位于靜水面內(nèi)船舯運動初始位置.O0-x0y0平面與靜水面重合，O0z0軸垂直向下.隨船運動坐標系原點選擇在船舯處，隨船體一起運動，Ox軸從船尾指向船首，Oy軸指向右舷.作用在船體上的水動力在x軸和y軸的分量分別為X和Y；繞船體z軸的轉(zhuǎn)首力矩為N.船體的平動速度U在隨船運動坐標系中的縱向和橫向分量為u和v，其中：u=Ucosβ,v=-Usinβ，U=|U|，β是漂角.舵的法向力記為FN.

圖1 坐標系Fig.1 Coordinate systems

1.2 流體流動控制方程

船舶周圍三維黏性流場可由雷諾平均的連續(xù)性方程和動量方程(RANS方程)來描述：

(1)

(2)

1.3 螺旋槳體積力模型

本文采用Hough等提出的體積力模型[11]，將體積力以源項的形式加入控制方程，即式(2)中.該方法通過對螺旋槳所在區(qū)域施加推力和轉(zhuǎn)矩代替槳葉表面載荷，計算公式如下：

(3)

(4)

式中：fbx和fbθ分別為軸向力和切向力在螺旋槳計算域內(nèi)的分布；rc為螺旋槳區(qū)域內(nèi)任意一點至螺旋槳軸線的距離；RH為槳轂半徑；RP為螺旋槳半徑；TP為螺旋槳推力；QP為螺旋槳轉(zhuǎn)矩；ΔP是螺旋槳盤面厚度，一般取值為槳直徑的2%.在實際計算中，施加體積力的過程為將推力和轉(zhuǎn)矩施加給螺旋槳作用區(qū)域，進而對控制方程進行求解計算.

2 數(shù)值求解

2.1 研究對象

選取KCS船模為研究對象，該船型是國際專題研討會SIMMAN2014[12]所采用的標準船型之一，船后舵和槳的幾何數(shù)據(jù)也由SIMMAN2014提供.數(shù)值研究采用的縮尺比和主尺度如表1所示.選取日本海上技術(shù)安全研究所(NMRI)提供的模型試驗數(shù)據(jù)(縮尺比1∶75.50)[13]進行比較，以驗證本文的數(shù)值方法.

表1 KCS船主尺度Tab.1 Main particulars of KCS

2.2 計算域及邊界條件

計算域和邊界條件設(shè)置如圖2所示，計算域總長度為5Lpp，寬度為2.5Lpp，水深為2.5Lpp，自由面以上高度為1.2Lpp.邊界條件：船體和舵表面設(shè)置為固壁邊界條件；除了出口設(shè)置為壓力出口外，其他計算域外面的5個面均設(shè)置為速度入口.在兩個側(cè)面、入口和出口處，設(shè)置人工數(shù)值阻尼來降低數(shù)值振蕩.

2.3 網(wǎng)格劃分

本文應(yīng)用CFD軟件STAR CCM+[14]，選取切邊網(wǎng)格生成器(Trimmed Mesher)，通過切割幾何表面，自動生成以六面體為主的網(wǎng)格.通過棱柱層生成器(Prism Layer Mesher)在壁面邊界上產(chǎn)生正交棱柱網(wǎng)格，從而精確求解流場.采用網(wǎng)格質(zhì)量修正數(shù)值模型(Cell Quality Remediation)來提高求解精度.采用混合壁面函數(shù)來處理壁面邊界條件.設(shè)置計算域網(wǎng)格初始目標尺寸(Target Size)為Lpp/30，在自由面、船首、船尾、舵以及槳周圍進行網(wǎng)格加密.棱柱層網(wǎng)格有6層，總棱柱層厚度為Lpp/150，拉伸比為1.5，無量綱化的壁面到第一層網(wǎng)格的距離y+≈25.由于采用混合壁面函數(shù)，生成的壁面邊界層厚是可接受的.在應(yīng)用體積力法(記為BF法)進行螺旋槳建模時，對激勵盤位置處的網(wǎng)格進行適當加密.在采用實槳法(記為RP法)進行螺旋槳建模時，對實槳進行網(wǎng)格離散，實現(xiàn)螺旋槳旋轉(zhuǎn)采用的方法是在槳附近產(chǎn)生一個圓柱形的運動域，定義這個運動域的旋轉(zhuǎn)速度，在其表面設(shè)置滑移壁面邊界條件，運動域與外部區(qū)域之間使用交界面邊界條件，使得螺旋槳域和整體域耦合，這種方法稱作滑移網(wǎng)格法.圖3給出了船體附近的網(wǎng)格結(jié)構(gòu).基于BF法和RP法得到整個計算域的網(wǎng)格總數(shù)分別為1.5×106和1.8×106.

圖2 計算域和邊界條件Fig.2 Computational domain and the boundary conditions

圖3 船體附近網(wǎng)格Fig.3 Grid structure around the ship

2.4 數(shù)值方法

本文應(yīng)用STAR CCM+進行流場數(shù)值模擬，基于有限體積法離散方程，其中瞬變項采用一階全隱式方案進行離散，對流項采用二階迎風格式，耗散項離散后存在一個二次梯度(交叉擴散)項，采用二階的高斯-最小二乘混合法求解梯度.采用分離流模型結(jié)合SIMPLE算法和Rhie-Chow修正方法分別求解壓力和速度項.通過VOF(Volume of Fluid)方法構(gòu)造自由面.船體的運動由動態(tài)流固耦合運動模塊(DFBI)計算獲得.在基于RP法的非定常計算過程中，為保證求解精度，時間步長小于1/(50nP)，其中nP是螺旋槳轉(zhuǎn)速，而BF法的螺旋槳用體積力代替，可以適當放大時間步長.所有計算都基于兩臺戴爾高性能計算服務(wù)器(36核/72線程，2.32 GHz).

3 數(shù)值結(jié)果分析

3.1 水動力結(jié)果

圖4 全附體KCS船模斜拖試驗水動力結(jié)果比較Fig.4 Comparison of hydrodynamic forces in oblique-towing tests for the fully appended KCS model

β/(°)BF法線程求解器總計算時間/sRP法線程求解器總計算時間/s065.06×104301.00×1052108.67×104142.49×105488.69×104201.48×105688.68×104162.28×105861.84×105201.57×1051261.84×105301.85×1051661.84×105301.76×10520367.68×104501.07×105平均10.751.18×10526.251.69×105

圖5 螺旋槳尾流軸向速度比較 (x/Lpp=0.488 3)Fig.5 Comparison of axial velocity contours in propeller slipstream (x/Lpp=0.488 3)

3.2 流場結(jié)果

進一步地，從流場細節(jié)進行比較分析.由于舵在螺旋槳后面，所以準確預報槳下游處流場是很重要的.圖5分別給出了基于BF法和RP法獲得的β=0°和6°時槳盤面后方(x/Lpp=0.488 3)流體的軸向速度u′(u′=u/U)分布；為討論方便，圖5(a)中定義了4個象限(I到 IV).當β=0°時，對于RP法來說，由于是右旋槳，槳葉從 III、IV 象限運動到I、II 象限.另一方面，船體的運動使得尾流場存在一對內(nèi)旋反向?qū)ΨQ的舭渦， 導致橫向速度方向和槳 III、IV 象限的橫向速度方向相同，和I、II 象限相反，致使槳葉剖面在 III 和 IV 象限的攻角減小，相對應(yīng)地在I、II 象限攻角增大，從而使得 III、IV 象限的軸向速度小于I、II 象限的軸向速度(圖5(b)).而對于BF法，反而是I、II 象限的軸向速度稍小于 III、IV 象限的軸向速度(圖5(a)).當β=6°時，尾流場與直航工況不同，在船體背風面，產(chǎn)生明顯的鉤狀尾流.在槳盤面內(nèi)，軸向速度最大值并沒有因為漂角的影響改變，只是其分布由于入流速度分布的改變而稍有改變(圖5(c)、(d)).

圖6 船尾表面壓力和中縱剖面軸向速度比較Fig.6 Comparison of pressure on the stern region and axial velocity contours on the central longitudinal section

圖6給出了β=0° 和6° 時船尾附近中縱剖面軸向速度和船、槳、舵表面的壓力系數(shù)Cp(Cp=p/(0.5ρU2))分布.從圖中可以看出，β=0° 時，在舵的前緣存在駐點，但是越往下游，由于舵的厚度減小，壓力逐漸減小，在舵后緣壓力又升高.在舵的左舷前緣上半部分可以看到高壓區(qū)，對應(yīng)的右舷同樣位置是一個低壓區(qū).在前緣的下半部分存在相反的壓力分布.β=6°時，由于左舷是迎風面，舵左舷的高壓區(qū)增大，低壓下降.相應(yīng)地，舵右舷的高壓區(qū)減小，低壓區(qū)增大.不論是直航還是斜航，RP法的誘導速度稍微大于BF法的誘導速度，導致RP法得到的舵的低壓區(qū)增大，也就是說RP法得到的舵法向力偏大，從圖4(e)中也能看到.總體而言，兩種方法的壓力和軸向速度分布是類似的.

雖然試驗方法具備較高的可靠性，但由于試驗費用高、周期長而受到限制，目前船-槳-舵系統(tǒng)斜拖試驗中的尾流軸向速度等流場細節(jié)是很難獲得的.由于缺乏流場細節(jié)測量數(shù)據(jù)，所以本文無法進行流場信息的驗證.

4 結(jié)論

本文分別采用體積力法和基于實槳的滑移網(wǎng)格法對螺旋槳進行建模，數(shù)值模擬了全附體KCS船模斜拖試驗的黏性流場.通過和國際發(fā)布的基準試驗數(shù)據(jù)對比，驗證了本文提出的數(shù)值方法的可靠性.此外，對比基于體積力法和滑移網(wǎng)格法的計算結(jié)果發(fā)現(xiàn)，體積力法雖然無法捕捉由于槳葉旋轉(zhuǎn)導致的流場細節(jié)變化，但對船-槳-舵整體的水動力和力矩的預報均比較準確.而采用基于實槳的滑移網(wǎng)格法，需要的計算時間要比前者長得多，因而計算成本很大.故對于船舶操縱運動水動力計算及操縱性預報問題，如果不需要詳細的螺旋槳流場信息，可采用體積力法進行螺旋槳建模.

數(shù)值模擬KCS船模斜拖試驗的計算結(jié)果表明，槳推力、船體的縱向力和橫向力、轉(zhuǎn)首力矩以及舵的法向力對應(yīng)的平均力同試驗結(jié)果的誤差是可接受的；但對于應(yīng)用體積力法進行螺旋槳建模得到的槳推力和船體縱向力的數(shù)值計算精度以及相應(yīng)的數(shù)值不確定度需要進一步開展研究.