劉 義, 鄒早建, 郭海鵬

(1. 中國(guó)船舶及海洋工程設(shè)計(jì)研究院， 上海 200011；2. 上海交通大學(xué) 船舶海洋與建筑工程學(xué)院， 上海 200240；3. 上海交通大學(xué) 海洋工程國(guó)家重點(diǎn)實(shí)驗(yàn)室， 上海 200240)

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

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

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

1 數(shù)學(xué)模型

1.1 坐標(biāo)系

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

圖1 坐標(biāo)系Fig.1 Coordinate systems

1.2 流體流動(dòng)控制方程

船舶周?chē)S黏性流場(chǎng)可由雷諾平均的連續(xù)性方程和動(dòng)量方程(RANS方程)來(lái)描述：

(1)

(2)

1.3 螺旋槳體積力模型

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

(3)

(4)

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

2 數(shù)值求解

2.1 研究對(duì)象

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

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

2.2 計(jì)算域及邊界條件

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

2.3 網(wǎng)格劃分

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

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

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

3.1 水動(dòng)力結(jié)果

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

β/(°)BF法線程求解器總計(jì)算時(shí)間/sRP法線程求解器總計(jì)算時(shí)間/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 流場(chǎng)結(jié)果

進(jìn)一步地，從流場(chǎng)細(xì)節(jié)進(jìn)行比較分析.由于舵在螺旋槳后面，所以準(zhǔn)確預(yù)報(bào)槳下游處流場(chǎng)是很重要的.圖5分別給出了基于BF法和RP法獲得的β=0°和6°時(shí)槳盤(pán)面后方(x/Lpp=0.488 3)流體的軸向速度u′(u′=u/U)分布；為討論方便，圖5(a)中定義了4個(gè)象限(I到 IV).當(dāng)β=0°時(shí)，對(duì)于RP法來(lái)說(shuō)，由于是右旋槳，槳葉從 III、IV 象限運(yùn)動(dòng)到I、II 象限.另一方面，船體的運(yùn)動(dòng)使得尾流場(chǎng)存在一對(duì)內(nèi)旋反向?qū)ΨQ的舭渦， 導(dǎo)致橫向速度方向和槳 III、IV 象限的橫向速度方向相同，和I、II 象限相反，致使槳葉剖面在 III 和 IV 象限的攻角減小，相對(duì)應(yīng)地在I、II 象限攻角增大，從而使得 III、IV 象限的軸向速度小于I、II 象限的軸向速度(圖5(b)).而對(duì)于BF法，反而是I、II 象限的軸向速度稍小于 III、IV 象限的軸向速度(圖5(a)).當(dāng)β=6°時(shí)，尾流場(chǎng)與直航工況不同，在船體背風(fēng)面，產(chǎn)生明顯的鉤狀尾流.在槳盤(pán)面內(nèi)，軸向速度最大值并沒(méi)有因?yàn)槠堑挠绊懜淖?，只是其分布由于入流速度分布的改變而稍有改?圖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í)船尾附近中縱剖面軸向速度和船、槳、舵表面的壓力系數(shù)Cp(Cp=p/(0.5ρU2))分布.從圖中可以看出，β=0° 時(shí)，在舵的前緣存在駐點(diǎn)，但是越往下游，由于舵的厚度減小，壓力逐漸減小，在舵后緣壓力又升高.在舵的左舷前緣上半部分可以看到高壓區(qū)，對(duì)應(yīng)的右舷同樣位置是一個(gè)低壓區(qū).在前緣的下半部分存在相反的壓力分布.β=6°時(shí)，由于左舷是迎風(fēng)面，舵左舷的高壓區(qū)增大，低壓下降.相應(yīng)地，舵右舷的高壓區(qū)減小，低壓區(qū)增大.不論是直航還是斜航，RP法的誘導(dǎo)速度稍微大于BF法的誘導(dǎo)速度，導(dǎo)致RP法得到的舵的低壓區(qū)增大，也就是說(shuō)RP法得到的舵法向力偏大，從圖4(e)中也能看到.總體而言，兩種方法的壓力和軸向速度分布是類(lèi)似的.

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

4 結(jié)論

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

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