陳 科，賴明雁

（1.中遠(yuǎn)海運發(fā)展股份有限公司，上海 200135；2.上海船舶研究設(shè)計院，上海 201203）

0 引 言

隨著計算機性能的提升和數(shù)值方法的發(fā)展，計算流體力學(xué) (Computational Fluid Dynamics, CFD) 在船舶海洋工程領(lǐng)域的應(yīng)用越來越廣，目前對許多流動問題的預(yù)報精度已達(dá)到較高的水平。從Tokyo 2015 CFD Workshop的計算結(jié)果來看，多家機構(gòu)的阻力、自航預(yù)報精度已在5%的誤差范圍內(nèi)[1]，在某種程度上能滿足工程設(shè)計要求。同時，CFD在耐波性和操縱性數(shù)值模擬方面也顯示出較高的可靠性，并可為一些非線性流動現(xiàn)象的數(shù)值計算提供參考[2-6]。近幾年，隨著綠色船舶和節(jié)能減排概念的深入發(fā)展，船舶設(shè)計者越來越重視船舶推進效率的提升和新型節(jié)能裝置的開發(fā)，迫切需要得到能較為可靠地預(yù)報船/槳相互干擾問題的數(shù)值處理方法。CFD因具有較高的靈活性和巨大的潛力而成為人們研究的熱點。

然而，船/槳相互干擾是一種比較復(fù)雜的非線性流動現(xiàn)象。從船體周圍的流動特點來看，艉部附近流線曲率變化較快，常伴隨流動分離現(xiàn)象，導(dǎo)致槳盤面處伴流較為復(fù)雜，使螺旋槳的推進性能有別于敞水。同時，螺旋槳的旋轉(zhuǎn)會產(chǎn)生抽吸現(xiàn)象，影響船體附近的流動，造成阻力增額。船/槳干擾問題的復(fù)雜性給CFD計算帶來很大挑戰(zhàn)。為研究船/槳干擾，相關(guān)研究人員開發(fā)出不同的數(shù)值處理方法，這些方法各有優(yōu)缺點，具體區(qū)別體現(xiàn)在對螺旋槳的處理上，主要分為體積力法和船/槳整體建模法2類。體積力法不考慮螺旋槳的幾何形狀，而是利用虛擬圓盤代替螺旋槳，在圓盤所在區(qū)域的控制方程中施加體積力源項。體積力法建模簡單、計算效率高，因此得到廣泛應(yīng)用[7-8]。然而，體積力法的計算精度有限，不能很好地體現(xiàn)出螺旋槳幾何形狀改變對流場的影響，且對槳盤面所在區(qū)域的流場預(yù)報不準(zhǔn)。隨著研究的不斷深入，相關(guān)研究人員開始嘗試對船/槳整體建立模型來求解[9]。這類方法中常用的是Moving Reference Frame（MRF）法和滑移網(wǎng)格法。在利用整體建模法計算時，須建立一個包含螺旋槳的小域。在利用 MRF法計算時，螺旋槳所在的區(qū)域處在一個相對旋轉(zhuǎn)的坐標(biāo)系內(nèi)，但網(wǎng)格是不動的；在利用滑移網(wǎng)格法計算時，螺旋槳按照相應(yīng)的轉(zhuǎn)速旋轉(zhuǎn)，流場處于瞬態(tài)變化中。整體來說，MRF法和滑移網(wǎng)格法各有優(yōu)缺點，其中MRF法的計算效率較高且穩(wěn)定，而滑移網(wǎng)格法則能更真實地模擬船/槳之間的非定常流動現(xiàn)象。目前，已有一些學(xué)者對這2種方法進行探討，如ZHANG[10]和楊春蕾等[11]從船體阻力、螺旋槳受力及伴流等方面分析MRF法與滑移網(wǎng)格法的區(qū)別。這些研究豐富了對船/槳相互干擾現(xiàn)象的認(rèn)識，但對自航因子預(yù)報及螺旋槳梢渦等流動特征方面的分析還比較少。另外，船/槳相互干擾計算的可靠性還需更多的驗證和分析。

本文針對國際標(biāo)準(zhǔn)船模韓國船舶與海洋工程研究所的集裝箱船（KRISO Container Ship，KCS），采用MRF法和滑移網(wǎng)格法對船/槳相互干擾的流動問題進行計算，將計算結(jié)果與試驗結(jié)果和其他數(shù)值結(jié)果相比較，從伴流、船體表面壓力、螺旋槳梢渦、船體和螺旋槳受力及自航因子等方面分析MRF法和滑移網(wǎng)格法的預(yù)報精度及船/槳之間相互干擾的流動特點，為工程實際應(yīng)用提供參考。

1 控制方程及數(shù)值方法

數(shù)值模擬采用商業(yè)CFD軟件STAR-CCM+，計算中控制方程采用Reynolds平均Navier-Stokes方程（即RANS方程），湍流模型采用SST k-ω兩方程湍流模型，自由面基于VOF（Volume Of Fluid）方法求解。對流項離散采用二階精度的迎風(fēng)格式。

控制方程包括連續(xù)性方程和動量守恒方程，對于不可壓縮流體，有

式(1)和式(2)中：i和j分別為變量在坐標(biāo)系中第i方向及第j方向的分量；P為壓力；u為速度矢量；ρ為流體密度；f為質(zhì)量力。

SST k-ω模型是一種結(jié)合k-ε模型和原始k-ω模型的混合方法，湍動能k和湍動耗散率ω的求解方程為

VOF控制方程為

2 數(shù)值計算結(jié)果及分析

2.1 船型參數(shù)

計算模型為國際標(biāo)模KCS，該模型的原型是韓國KRISO（Koera Research Institute of Ships and Ocean Engineering）設(shè)計的集裝箱船，船體幾何形狀見圖1。KCS被國際學(xué)術(shù)界推薦用來進行CFD驗證和分析，有較為充分的試驗數(shù)據(jù)作為參考，包括靜水阻力、自航、耐波性及操縱性等。船體和槳的主要參數(shù)見表1。KCS實船長230m，設(shè)計航速為24kn，模型縮尺比為31.6。本文計算采用的航速為設(shè)計航速，對應(yīng)Fr=0.26。

圖1 KCS船體幾何形狀

表1 KCS船體和槳參數(shù)

2.2 邊界條件及網(wǎng)格

在 CFD計算中，計算域的選取不僅影響計算效率，還會影響計算精度。一般來說，計算域和邊界條件的指定依賴于計算問題本身。對于自航計算，由于船體左右流動不對稱，計算域?qū)挾缺瘸Ｒ?guī)阻力計算大。本文的計算域取為長方體，大小為：船前1倍船長，兩側(cè)和底部1.5倍船長。船前兩側(cè)和上下面均為來流邊界條件，取x軸負(fù)向為壓力出口邊界。計算網(wǎng)格為非結(jié)構(gòu)六面體網(wǎng)格，為更好地捕捉船體附近的流動特征，特別是艉部船/槳相互干擾區(qū)域的流動細(xì)節(jié)，對艉部和螺旋槳附近區(qū)域進行加密；為更好地模擬螺旋槳泄渦，對槳葉的導(dǎo)邊和隨邊均進行加密。計算域網(wǎng)格劃分見圖2，網(wǎng)格總數(shù)為271萬個，其中船體部分網(wǎng)格為220萬個，螺旋槳所在區(qū)域網(wǎng)格為51萬個。體網(wǎng)格類型主要以六面體為主，為更好地捕捉物面附近的流動，在船體和螺旋槳表面生成邊界層網(wǎng)格。

圖2 計算域網(wǎng)格劃分

2.3 計算結(jié)果及分析

為分析船/槳干擾問題及自航因子計算，首先對KCS在裸船條件下的靜水阻力進行計算，計算網(wǎng)格配置參數(shù)與“2.2”節(jié)中的敘述相似。表2給出阻力的計算結(jié)果和試驗結(jié)果，二者能較好地吻合。圖3為KCS靜水狀態(tài)航行時波形，對比了船體附近的自由面興波，可看到通過數(shù)值計算得到的自由面波形能與試驗結(jié)果較好地吻合，數(shù)值計算得到的波峰和波谷位置及波形都能與試驗結(jié)果相一致，表明本文采用的數(shù)值方法具有一定的可靠性。

表2 KCS靜水阻力計算結(jié)果

KCS自航試驗是在實船自航點下進行的，因此在船模尺度計算中需添加一項恒定的強制力，或稱摩擦力修正力（Skin Friction Correction，SFC，以Fc表示），用于修正尺度效應(yīng)引起的摩擦力誤差。Fc定義為

圖3 KCS靜水狀態(tài)航行時波形

式(6)中：分別為船模和實船的摩擦阻力系數(shù)；k為形狀因子，對于KCS而言，k+1=1.1；U0為航速；SW為濕表面積。

計算中螺旋槳的轉(zhuǎn)速與模型試驗一致，即n=9.50r/s。表3給出分別用MRF法和滑移網(wǎng)格法得到的推進及自航因子，并與試驗結(jié)果和其他數(shù)值結(jié)果相對比，其中CASTRO等[12]的計算采用IOWA-CFDSHIP求解器，并在計算中使用DES（Detached Eddy Simulation）模擬湍流流動，螺旋槳旋轉(zhuǎn)采用over-set網(wǎng)格方式來實現(xiàn)。需說明的是，本文的自航因子計算采用的是強制自航法，即：首先計算 2個轉(zhuǎn)速（n=9.50r/s和n=9.75r/s）下的自航，記錄相應(yīng)的強制力、螺旋槳推力系數(shù)KT和扭矩系數(shù)KQ；然后通過強制力插值得到自航時的螺旋槳轉(zhuǎn)速、推力和扭矩等信息，插值原則為保證強制力等于試驗的強制力（即30.25N）。自航因子采用等推力法計算得到，需說明的是，螺旋槳敞水?dāng)?shù)據(jù)采用試驗測量結(jié)果。

從轉(zhuǎn)速的計算結(jié)果來看，滑移網(wǎng)格法精度較高，MRF法給出的轉(zhuǎn)速偏高。這主要是因為利用MRF法計算得到的船體阻力較大，因此在進行強制力插值時得到的轉(zhuǎn)速較大。利用MRF法和滑移網(wǎng)格法得到的KT均與試驗結(jié)果相近，KQ相差略大，即通過數(shù)值模擬得到的螺旋槳推進效率較低，CASTRO等[12]的計算結(jié)果存在同樣的現(xiàn)象。這可能是因為艉流場流動比較復(fù)雜，存在邊界層的分離，數(shù)值模擬未能很好地捕捉這些流動特征。由表3中的自航因子可知，利用滑移網(wǎng)格法得到的自航因子與試驗結(jié)果能較好地吻合，相對而言，利用MRF法計算得到的推力減額偏大，推進效率偏低。對于伴流分?jǐn)?shù)和相對旋轉(zhuǎn)效率等因子，利用MRF法和滑移網(wǎng)格法得到的結(jié)果都能與試驗結(jié)果相一致。

表3 利用MRF法和滑移網(wǎng)格法得到的自航計算結(jié)果

圖4為靜水阻力計算和自航計算得到的船體表面動壓力分布對比。由圖4可知，在自航狀態(tài)下，受螺旋槳抽吸影響，船體表面的動壓力比靜水阻力低，導(dǎo)致自航時船體阻力較大，在功率換算中等價于推力減額。利用MRF法和滑移網(wǎng)格法得到的計算結(jié)果都能較好地反映該船/槳干擾的流動現(xiàn)象。

圖5為利用試驗、MRF法和滑移網(wǎng)格法得到的螺旋槳下游x/LPP=0.9941處的伴流對比。從整體上看，通過數(shù)值模擬得到的伴流場與通過試驗測量得到的結(jié)果相一致；但從細(xì)節(jié)上看，利用滑移網(wǎng)格法得到的結(jié)果與試驗結(jié)果比較吻合。這里需指出的是，在利用MRF法計算時螺旋槳不動，因此螺旋槳槳葉的相位會對結(jié)果有一定影響。圖6為艉部流場中Q等值面。對比采用MRF法和滑移網(wǎng)格法得到的計算結(jié)果可看到，由于MRF法僅在螺旋槳旋轉(zhuǎn)域內(nèi)部施加體積力，因此得到的螺旋槳梢渦在小范圍內(nèi)呈旋轉(zhuǎn)狀態(tài)，而在距離螺旋槳稍遠(yuǎn)處隨來流向下游水平泄出，在螺旋槳計算子域和外域交界處存在渦的非物理變形。相對而言，利用滑移網(wǎng)格法得到的渦結(jié)構(gòu)與物理直觀更為接近。

圖4 船體表面動壓分布

圖5 螺旋槳下游的軸向速度分布（x/LPP=0.9941)

圖6 艉部流場中Q等值面

3 結(jié) 語

本文基于MRF法和滑移網(wǎng)格法對國際標(biāo)模KCS的船/槳干擾問題進行數(shù)值研究。從計算結(jié)果來看，螺旋槳的存在會改變船體表面的壓力分布，造成船體阻力增大，即出現(xiàn)推力減額現(xiàn)象。利用MRF法和滑移網(wǎng)格法都能較好地反映該流動現(xiàn)象。從自航因子結(jié)果來看，滑移網(wǎng)格法給出的自航因子能與試驗結(jié)果和其他數(shù)值結(jié)果相吻合，而采用MRF法所得結(jié)果的誤差較大，這主要是因為利用MRF法得到的船體阻力偏大。從伴流結(jié)果來看，滑移網(wǎng)格法與試驗吻合較好。此外，在螺旋槳附近，利用MRF法和利用滑移網(wǎng)格法得到的槳葉梢渦形狀基本一致；在遠(yuǎn)離螺旋槳處，利用滑移網(wǎng)格法得到的計算結(jié)果更符合物理直觀。

【 參 考 文 獻 】

[1] LARSSON L, STERN F, VISONNEAU M, et al.Proceedings of CFD workshop Tokyo 2015[C].National Maritime Research Institute, 2015.

[2] CARRICA P M, FU H, STERN F.Computations of self-propulsion free to sink and trim and of motions in head waves of the KRISO Container Ship (KCS) model[J].Applied Ocean Research, 2011, 33: 309-320.

[3] SHEN Z R, WAN D C.RANS computations of added resistance and motions of a ship in head waves[J].International Journal of Offshore and Polar Engineering, 2013, 23 (4): 263-271.

[4] SHEN Z R, YE H X, WAN D C.URANS simulations of ship motion responses in long-crest irregular waves[J].Journal of Hydrodynamics, 2014, 26 (3): 436-446.

[5] WANG J H, ZHAO W W, WAN D C.Free maneuvering simulation of ONR Tumblehome using overset grid method in naoe-FOAM-SJTU solver[C]//31stSymposium on Naval Hydrodynamics, September 11-16, 2016, California.

[6] HU C H, SUEYOSHI M, KASHIWAGI M.Numerical simulation of strongly nonlinear wave-ship interaction by CIP/Cartesian grid method[J].International Journal of offshore and Polar Engineering, 2010: 81-87.

[7] 吳召華，陳作鋼，代燚.基于體積力法的船體自航性能數(shù)值預(yù)報[J].上海交通大學(xué)學(xué)報，2013, 47 (6): 943-949.

[8] SIMONSEN C D, STERN F.RANS maneuvering simulation of esso osaka with rudder and a body-force propeller[J].Journal of Ship Research, 2005, 49(2):98-120.

[9] CARRICA P M, SADAT-HOSSEINI H, STERN F.CFD analysis of broaching for a model surface combatant with explicit simulation of moving rudders and rotating propellers[J].Computers & Fluids, 2012, 53: 117-132.

[10] ZHANG Z R.Verification and validation for RANS simulation of KCS container ship without/with propeller[J].Journal of Hydrodynamics, 2010, 22 (5): 932-939.

[11] 楊春蕾，朱仁傳，繆國平，等.基于CFD方法的船/槳/舵干擾數(shù)值模擬[J].水動力學(xué)研究與進展，2011, 26 (6): 667-673.

[12] CASTRO A M, CARRICA P M, STERN F.Full scale self-propulsion computations using discretized propeller for the KRISO Container Ship KCS[J].Computers & Fluids, 2011, 51: 35-47.