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

劉 義, 鄒早建, 郭海鵬

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

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

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

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

1 數學模型

1.1 坐標系

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

圖1 坐標系Fig.1 Coordinate systems

1.2 流體流動控制方程

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

(1)

(2)

1.3 螺旋槳體積力模型

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

(3)

(4)

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

2 數值求解

2.1 研究對象

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

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

2.2 計算域及邊界條件

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

2.3 網格劃分

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

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

圖3 船體附近網格Fig.3 Grid structure around the ship

2.4 數值方法

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

3 數值結果分析

3.1 水動力結果

圖4 全附體KCS船模斜拖試驗水動力結果比較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 流場結果

進一步地，從流場細節進行比較分析.由于舵在螺旋槳后面，所以準確預報槳下游處流場是很重要的.圖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 象限.另一方面，船體的運動使得尾流場存在一對內旋反向對稱的舭渦， 導致橫向速度方向和槳 III、IV 象限的橫向速度方向相同，和I、II 象限相反，致使槳葉剖面在 III 和 IV 象限的攻角減小，相對應地在I、II 象限攻角增大，從而使得 III、IV 象限的軸向速度小于I、II 象限的軸向速度(圖5(b)).而對于BF法，反而是I、II 象限的軸向速度稍小于 III、IV 象限的軸向速度(圖5(a)).當β=6°時，尾流場與直航工況不同，在船體背風面，產生明顯的鉤狀尾流.在槳盤面內，軸向速度最大值并沒有因為漂角的影響改變，只是其分布由于入流速度分布的改變而稍有改變(圖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° 時船尾附近中縱剖面軸向速度和船、槳、舵表面的壓力系數Cp(Cp=p/(0.5ρU2))分布.從圖中可以看出，β=0° 時，在舵的前緣存在駐點，但是越往下游，由于舵的厚度減小，壓力逐漸減小，在舵后緣壓力又升高.在舵的左舷前緣上半部分可以看到高壓區，對應的右舷同樣位置是一個低壓區.在前緣的下半部分存在相反的壓力分布.β=6°時，由于左舷是迎風面，舵左舷的高壓區增大，低壓下降.相應地，舵右舷的高壓區減小，低壓區增大.不論是直航還是斜航，RP法的誘導速度稍微大于BF法的誘導速度，導致RP法得到的舵的低壓區增大，也就是說RP法得到的舵法向力偏大，從圖4(e)中也能看到.總體而言，兩種方法的壓力和軸向速度分布是類似的.

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

4 結論

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

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