柯梟冰，羅薇，趙小仨，徐海祥

1 武漢理工大學(xué)高性能船舶技術(shù)教育部重點實驗室，湖北武漢430063

2 武漢理工大學(xué)交通學(xué)院，湖北武漢430063

0 引 言

人類自進入工業(yè)時代以來，物質(zhì)財富極大地增長，但陸上資源和能源供應(yīng)卻日趨緊張，文明的可持續(xù)發(fā)展還有賴于陸地經(jīng)濟向海洋經(jīng)濟的轉(zhuǎn)型。而海洋環(huán)境復(fù)雜多變，海洋資源和能源的開發(fā)有其特殊的困難，需要先進的技術(shù)和設(shè)備作為支撐，其中動力定位技術(shù)就是海洋開發(fā)技術(shù)的一個典型代表。

供應(yīng)船作為動力定位技術(shù)的常用船型，其使用日益增加。在建立供應(yīng)船的運動數(shù)學(xué)模型時，需要確定模型中的水動力導(dǎo)數(shù)，本文將主要研究其中的位置水動力導(dǎo)數(shù)的求取方法。

在船舶設(shè)計階段，確定船舶水動力導(dǎo)數(shù)的方法一般有：經(jīng)驗公式法、約束模試驗法、理論和數(shù)值計算法，以及自航模試驗（實船試驗）加系統(tǒng)辨識法。其中經(jīng)驗公式法最方便、快捷，但使用范圍有一定的限制，對于合適的船型，所得到的結(jié)果能夠達到工程需要的精度。由于動力定位供應(yīng)船一般都裝有呆木，目前，還未見有關(guān)這一類船型的水動力導(dǎo)數(shù)經(jīng)驗公式方面的文獻。

鑒于此，本文將通過船型變換生成系列供應(yīng)船船型，采用CFD 方法計算系列供應(yīng)船型的位置水動力導(dǎo)數(shù)Yv和Nv，再采用多元回歸分析建立位置水動力導(dǎo)數(shù)關(guān)于船型參數(shù)的回歸公式，以為設(shè)計階段供應(yīng)船水動力導(dǎo)數(shù)的求解提供一定的參考。

1 數(shù)值計算方法

1.1 坐標(biāo)系

坐標(biāo)系如圖1 所示。其中，o 點位于船舯，縱軸x 的正方向指向船艏，橫軸y 的正方向指向左舷，按照右手法則，z軸正向垂直向上。

圖1 隨船坐標(biāo)系Fig.1 Coordinate system on board

1.2 控制方程

對三維粘性不可壓縮的N-S 方程和連續(xù)性方程進行時均化處理［1］，得到

式中：ui為平均速度分量，其中u1=u，u2=v，u3=w，i 和j 的取值范圍為（1，2，3）；p 為平均壓力；ρ 為流體密度；Fi為體積力；μ 為流體動力粘性系數(shù)；為雷諾應(yīng)力項。式（2）即為雷諾時均N-S 方程，簡稱RANS 方程。

1.3 湍流模型

2 KVLCC2 計算驗證

2.1 計算模型與網(wǎng)格劃分

為驗證本文計算方法的有效性，首先使用該方法計算ITTC 推薦的船型KVLCC2 在斜航狀況下所受的水動力［2-3］，并與INSEAN 水池斜航試驗結(jié)果［4］進行對比。KVLCC2 船型的主要幾何參數(shù)如表1 所示，其中船模縮尺比為45.714。

表1 KVLCC2 船模主要參數(shù)Tab.1 The main parameters of KVLCC2 model

計算域如圖2 所示。在上游1.5LPP處建立入流面；下游3LPP處建立出流面；左、右邊界距船模中縱剖面各為1.5LPP，水深吃水比h/T=8.3，傅氏數(shù)Fn=0.064。

圖2 計算域Fig.2 The computational domain

KVLCC2 船的艉部因有艉軸出口，因而在全船劃分結(jié)構(gòu)化網(wǎng)格難以實現(xiàn)?？稍诖禾巹澐址墙Y(jié)構(gòu)化的四面體網(wǎng)格，其余各處采用分塊結(jié)構(gòu)化網(wǎng)格，如圖3 所示。網(wǎng)格數(shù)共413.5 萬。在界層內(nèi)，使用Boundary Layer 命令劃分多層網(wǎng)格來捕捉船模邊界層內(nèi)的流動。

2.2 邊界條件及數(shù)值計算方法

所謂邊界條件，是指在計算域的邊界上，所求解的變量或者一階導(dǎo)數(shù)隨地點和時間變化的規(guī)律。本文邊界條件設(shè)置如下［1］：

1）入流面及左、右表面：VELOCITY_INLET；

圖3 艏、艉部網(wǎng)格Fig.3 Grids of bow and stern

2）流動出流表面：OUTFLOW；

3）船體表面：WALL；

4）計算域的上、下表面：SYMMETRY；

5）分塊劃分結(jié)構(gòu)化網(wǎng)格時，塊與塊之間重合而分屬不同的塊的面：INTERFACE。

湍流模型選取SST k-ω 模型，壓力速度耦合選取SIMPLEC 算法，壓力插值采用Standard 格式，動量、湍動能和湍流耗散率插值采用Second Order Upwind 格式。

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

保持x 方向速度不變，變化y 方向的速度來形成不同的漂角（v=U sin β，其中v 為橫向速度，U為縱向速度）。根據(jù)相對運動原理，船模位置保持不變，水流以一定的速度流向船模［5-6］?；贔LUENT 軟件平臺，計算不同漂角下船舶所受的橫向力Y 和轉(zhuǎn)艏力矩N。本文計算了0.268°，0.55°，1.05°，1.8°和3.8°這5 個不同漂角下各船型的橫向力Y 和轉(zhuǎn)艏力矩N，計算結(jié)果與INSEAN 水池約束模型試驗數(shù)據(jù)的對比如圖4 所示。

圖4 KVLCC2 船模CFD 結(jié)果與試驗數(shù)據(jù)對比Fig.4 The comparison of CFD data and experimental data of KVLCC2 model

由圖4 可以看到，橫向力Y 的計算值與試驗值吻合良好，誤差較小。在1.05°時，試驗值出現(xiàn)了微小的波動，而計算值則比較平穩(wěn)，這可能是因為數(shù)值模擬受到的干擾沒有試驗狀態(tài)下的復(fù)雜，比較符合實際情況。轉(zhuǎn)艏力矩N 在小角度（0.268°，0.55°，1.05°）時，其數(shù)值計算結(jié)果與試驗值吻合良好，而在稍大角度（1.8°，3.8°）時二者的間隔則較大，但總體趨勢還是一致的［7-8］。

3 系列船型計算方案

3.1 系列船型生成

變換的母型船為一艘75 m 的動力定位供應(yīng)船，船?？s尺比為20，主要幾何參數(shù)如表2 所示，供應(yīng)船側(cè)視圖如圖5 所示。

表2 供應(yīng)船的主要參數(shù)Tab.2 The main parameters of supply vessel

圖5 動力定位供應(yīng)船F(xiàn)ig.5 Dynamic positioning supply vessel

對50 條供應(yīng)船進行了統(tǒng)計，其主要參數(shù)取值范 圍 為：3.2＜L/B＜5，9.5＜L/d＜20，2.5＜B/d＜4.5，0.6＜Cb＜0.8 。在船型變換過程中，保持排水量和船長不變，呆木相對于船模中縱剖面的位置和面積也均保持不變，變換的模型利用Maxsurf 軟件生成。變換方案［9-10］如表3 所示。

表3 船型變換方案Tab.3 Hull form transformation scheme

在生成的25 個系列船型（除去母型船共24個）中，減去吃水過小和過大的船型，并將非常接近的方案取其中之一后，剩下的13 個船型即為需要進行計算的模型。

3.2 系列船型計算結(jié)果及無因次化

系列船型的斜航數(shù)值計算方法及步驟與驗證的船型KVLCC2 相同，僅在設(shè)置計算參數(shù)時考慮了各自船型的不同特點。需要說明的是，在劃分網(wǎng)格時，13 個船型均采用的是結(jié)構(gòu)化網(wǎng)格劃分方法，其分塊方式與網(wǎng)格節(jié)點分布一致，并保證y+值在一定的范圍內(nèi)以保證網(wǎng)格相似性。

下面，以方案A1 為例來說明求取水動力導(dǎo)數(shù)的步驟。

根據(jù)船舶水動力導(dǎo)數(shù)的物理意義，可以由表4所示的不同漂角β（橫向速度）下船模所受的橫向力和力矩來求取位置水動力導(dǎo)數(shù)Yv和Nv。以A1船為例，其在不同橫向速度下的橫向力和力矩曲線如圖6 所示。分別求取Y 關(guān)于v 的曲線和N 關(guān)于v 的曲線在零點的切線的斜率，該斜率值即為對應(yīng)的水動力導(dǎo)數(shù)。

表4 方案A1 的橫向力Y 與轉(zhuǎn)艏力矩N 計算結(jié)果Tab.4 Computational results of Y and N of hull form A1

圖6 A1 船的橫向力Y 和力矩NFig.6 Y and N of hull form A1

根據(jù)以上方法求取其它船型的位置水動力導(dǎo)數(shù)，并將結(jié)果進行無因次化，如表5 所示。其中，無因次化的規(guī)則為：

式中：U 為x 方向的速度；L 為船長。

4 位置水動力導(dǎo)數(shù)回歸公式

根據(jù)表5 中各船型的水動力導(dǎo)數(shù)無因次值，進行多元線性回歸，分析擬合位置水動力導(dǎo)數(shù)的回歸公式［11-12］。首先，進行相關(guān)性分析。直接選取B/d，L/B，L/d，Cb為自變量進行擬合，則各變量的相關(guān)性如表6 所示。

表5 系列船型的位置水動力導(dǎo)數(shù)無因次值Tab.5 Computational results of Yv′ and Nv′

表6 相關(guān)性分析（自變量為B/d，L/B，L/d，Cb）Tab.6 Correlation analysis（independent variable：L/d，L/B，B/d，Cb）

表7 相關(guān)性分析（自變量為d/L 和（d/L）Cb（B/d））Tab.7 Correlation analysis（independent variable：d/L，（d/L）Cb（B/d））

下面，將給出由全部13 個樣本點擬合的回歸公式：

為了檢驗回歸公式的擬合質(zhì)量，先取12 個樣本點進行公式擬合，然后再用余下的那個樣本點進行檢驗。余下的樣本點是參數(shù)較為居中的B3。表8 給出了B3 的和根 據(jù)12個樣本點回歸公式預(yù)測的值與CFD 計算的或值的比較。

表8 回歸公式質(zhì)量驗證Tab.8 Quality validation of regression formula

由表中的誤差項可以看出，最小二乘法擬合公式的預(yù)測值與CFD 計算值相比，誤差較小，可以認為將樣本點進行最小二乘法回歸是合適的。

5 結(jié) 語

本文對動力定位船舶位置水動力導(dǎo)數(shù)的CFD回歸方法進行了研究：

首先，基于FLUENT 軟件平臺計算了KVLCC2 船的斜航流場，得到了橫向力和轉(zhuǎn)艏力矩，經(jīng)與INSEAN 水池試驗數(shù)據(jù)進行對比，證明用CFD 方法計算船舶斜航流場是有效的。

然后，通過船型變換技術(shù)得到動力定位供應(yīng)船（裝有呆木）的不同船型，并采用與上述相同的數(shù)值方法計算其不同漂角下的斜航流場，通過最小二乘法擬合得到位置水動力導(dǎo)數(shù)Yv和Nv以及其無因次值。

最后，通過多元線性回歸，得到了動力定位供應(yīng)船位置水動力導(dǎo)數(shù)關(guān)于參數(shù)和的回歸公式。初步驗證回歸公式擬合質(zhì)量良好，證明對位置水動力導(dǎo)數(shù)的結(jié)果進行線性回歸是合理的。求取的回歸公式可為設(shè)計階段供應(yīng)船位置水動力導(dǎo)數(shù)的確定提供簡便的方法。

［1］王福軍.計算流體動力學(xué)分析：CFD 軟件原理與應(yīng)用［M］.北京：清華大學(xué)出版社，2004.

［2］田喜民，鄒早建，王化明. KVLCC2 船模斜航運動功粘性流場及水動力數(shù)值計算［J］.船舶力學(xué)，2010，14（8）：834-840.TIAN Ximin，ZOU Zaojian，WANG Huaming. Compu?tation of the viscous flow and hydrodynamic forces on a KVLCC2 model in oblique motion［J］. Journal of Ship Mechanics，2010，14（8）：834-840.

［3］楊勇.非定常操縱運動船體水動力數(shù)值計算［D］.上海：上海交通大學(xué)，2011.

［4］Workshop on Verification and Validation of Ship Ma?neuvering Simulation Methods［EB/OL］.（2008-04-16）.［2013-03-15］. http：//www.simman2008.dk/KV?LCC/KVLCC2/tanker2.html.

［5］李冬荔，楊亮，張洪雨，等. 基于CFD 方法的船舶操縱性能預(yù)報［J］. 武漢理工大學(xué)學(xué)報，2009，31（24）：120-123.LI Dongli，YANG Liang，ZHANG Hongyu，et al.Predic?tion of ship maneuverability based on CFD method［J］.Journal of Wuhan University of Technology，2009，31（24）：120-123.

［6］李冬荔，楊亮，聶武. 操縱運動船舶的水動力計算研究［J］.船舶工程，2009，31（2）：8-11，23.LI Dongli，YANG Liang，NIE Wu.Computational inves?tigation of hydrodynamic forces around maneuvering ship［J］.Ship Engineering，2009，31（2）：8-11，23.

［7］劉山. 基于CFD 技術(shù)數(shù)值模擬平面運動機構(gòu)試驗［D］.武漢：武漢理工大學(xué)，2012.

［8］楊波，萬林，王驍，等.純橫蕩和旋臂試驗的數(shù)值模擬［J］.艦船科學(xué)技術(shù)，2008，30（4）：138-141.YANG Bo，WAN Lin，WANG Xiao，et al. Numerical simulation of pure sway and rotating arm test［J］. Ship Science and Technology，2008，30（4）：138-141.

［9］王銀. 艦船方案設(shè)計中的耐波性預(yù)報模型研究［D］.哈爾濱：哈爾濱工程大學(xué)，2012.

［10］項久洋，毛筱菲. 基于主要尺度要素的船型變換［J］.中國艦船研究，2008，3（4）：15-18，25.XIANG Jiuyang，MAO Xiaofei. The hull variation based on the main hull parameters［J］. Chinese Jour?nal of Ship Research，2008，3（4）：15-18，25.

［11］袁益雷. 基于新細長船興波阻力理論的系列60 阻力圖譜之拓展研究［D］. 哈爾濱：哈爾濱工程大學(xué)，2008.

［12］楊鹽生.船舶阻力系數(shù)和推力系數(shù)計算的數(shù)據(jù)庫方法［J］.大連海事大學(xué)學(xué)報，1995，21（4）：14-17.YANG Yansheng. Calculating ship resistance and thrust coefficient based on database［J］. Journal of Dalian Maritime University，1995，21（4）：14-17.

［13］方崇智，蕭德云. 過程辨識［M］. 北京：清華大學(xué)出版社，1988.