基于CFD方法求取供應船位置水動力導數回歸公式

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

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

2 武漢理工大學交通學院，湖北武漢430063

0 引 言

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

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

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

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

1 數值計算方法

1.1 坐標系

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

圖1 隨船坐標系Fig.1 Coordinate system on board

1.2 控制方程

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

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

1.3 湍流模型

2 KVLCC2 計算驗證

2.1 計算模型與網格劃分

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

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

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

圖2 計算域Fig.2 The computational domain

KVLCC2 船的艉部因有艉軸出口，因而在全船劃分結構化網格難以實現。可在船艉處劃分非結構化的四面體網格，其余各處采用分塊結構化網格，如圖3 所示。網格數共413.5 萬。在界層內，使用Boundary Layer 命令劃分多層網格來捕捉船模邊界層內的流動。

2.2 邊界條件及數值計算方法

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

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

圖3 艏、艉部網格Fig.3 Grids of bow and stern

2）流動出流表面：OUTFLOW；

3）船體表面：WALL；

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

5）分塊劃分結構化網格時，塊與塊之間重合而分屬不同的塊的面：INTERFACE。

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

2.3 計算結果及分析

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

圖4 KVLCC2 船模CFD 結果與試驗數據對比Fig.4 The comparison of CFD data and experimental data of KVLCC2 model

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

3 系列船型計算方案

3.1 系列船型生成

變換的母型船為一艘75 m 的動力定位供應船，船模縮尺比為20，主要幾何參數如表2 所示，供應船側視圖如圖5 所示。

表2 供應船的主要參數Tab.2 The main parameters of supply vessel

圖5 動力定位供應船Fig.5 Dynamic positioning supply vessel

對50 條供應船進行了統計，其主要參數取值范 圍 為：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 系列船型計算結果及無因次化

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

下面，以方案A1 為例來說明求取水動力導數的步驟。

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

表4 方案A1 的橫向力Y 與轉艏力矩N 計算結果Tab.4 Computational results of Y and N of hull form A1

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

根據以上方法求取其它船型的位置水動力導數，并將結果進行無因次化，如表5 所示。其中，無因次化的規則為：

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

4 位置水動力導數回歸公式

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

表5 系列船型的位置水動力導數無因次值Tab.5 Computational results of Yv′ and Nv′

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

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

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

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

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

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

5 結 語

本文對動力定位船舶位置水動力導數的CFD回歸方法進行了研究：

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

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

最后，通過多元線性回歸，得到了動力定位供應船位置水動力導數關于參數和的回歸公式。初步驗證回歸公式擬合質量良好，證明對位置水動力導數的結果進行線性回歸是合理的。求取的回歸公式可為設計階段供應船位置水動力導數的確定提供簡便的方法。

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

［2］田喜民，鄒早建，王化明. KVLCC2 船模斜航運動功粘性流場及水動力數值計算［J］.船舶力學，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］楊勇.非定常操縱運動船體水動力數值計算［D］.上海：上海交通大學，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 方法的船舶操縱性能預報［J］. 武漢理工大學學報，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 技術數值模擬平面運動機構試驗［D］.武漢：武漢理工大學，2012.

［8］楊波，萬林，王驍，等.純橫蕩和旋臂試驗的數值模擬［J］.艦船科學技術，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］王銀. 艦船方案設計中的耐波性預報模型研究［D］.哈爾濱：哈爾濱工程大學，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］. 哈爾濱：哈爾濱工程大學，2008.

［12］楊鹽生.船舶阻力系數和推力系數計算的數據庫方法［J］.大連海事大學學報，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］. 北京：清華大學出版社，1988.