基于CFD的斜航船舶浮態數值預報

邵江明，應業炬

(浙江海洋學院船舶與海洋工程學院，浙江舟山 316022)

基于CFD的斜航船舶浮態數值預報

邵江明，應業炬

(浙江海洋學院船舶與海洋工程學院，浙江舟山 316022)

目前涉及斜航船舶的數值模擬計算鮮有在模擬過程中考慮到船舶浮態的影響，然而浮態變化對水動力數值計算精度的影響是不容忽視的。基于此，文章以系列60船斜航運動（傅汝德數Fr=0.316，斜航角β=10°）為算例，采用重疊網格技術處理斜航船的大幅運動問題，通過耦合求解船舶六自由度運動方程和流體流動方程的方法，計算出了斜航船的升沉量、縱傾和橫傾角度。結果顯示，考慮浮態的數值預報結果與愛荷華大學的試驗數據相比誤差基本控制在10%以內，體現出該數值計算預報方法具有較強的工程實用性。

斜航；浮態計算；重疊網格；數值模擬

船舶斜航即指船舶的航行方向與前方來流形成一定夾角的航行狀態，水面船舶大多數情況下都處于斜航狀態。當船舶斜航時，受到不對稱流的影響，船體受到較直航情況更加復雜的水動力影響。比如，在流場方面：斜航船的尾流場中流體分離作用更加明顯，伴流分數變大[1]；浮態方面，斜航船除了會出現升沉、縱傾運動外往往還會出現較為明顯的橫傾。橫傾作用對高重心船（集裝箱船或者滾裝船）的影響是非常惡劣的，它可能導致貨物碰撞甚至傾覆。因此對斜航船浮態的研究顯得格外重要。

目前對船舶斜航的研究主要有約束模試驗[2-3]、基于勢流理論的理論計算[4-5]和和基于黏性流理論數值模擬[6-7]。約束模試驗方法具有一定的精度優勢，但一般費時費力?；趧萘骼碚摰挠嬎銓κ莛ば杂绊戄^大的運動很難模擬準確，比如橫傾運動。近些年來，隨著計算機科學技術的飛速發展，基于計算流體動力學（CFD）的數值模擬方法得到越來越廣泛的應用，模擬中引入湍流模型概念充分考慮到流體的黏性作用。目前，這種方法已成為與約束模試驗并駕齊驅、相輔相成的一種水動力預報方法?；仡櫼酝鶎Υ靶焙降臄抵的M，可以發現，在絕大部分的研究中很少在模擬中考慮到船舶浮態影響。然而，浮態對水動力計算精度的影響是不容忽視的。劉晗，馬寧等[8]在對比由Fluent軟件計算得到的斜航KVLCC2船模水動力的數值結果和由循環水槽平面運動機構試驗得到的試驗結果中發現，在斜航角為20度時，不考慮浮態情況下的水動力系數同試驗值誤差達到20%之多。鄒璐[9]在研究數值計算精度中提到，對比考慮斜航船舶升沉、縱傾后的水動力計算值和不考慮浮態改變的情況，誤差在20%-25%之間。

在本文的數值模擬中，采用Star-ccm軟件，通過耦合求解船舶的運動方程和流體運動方程，計算出了斜航船的浮態變化，其中包括升沉、縱傾以及橫傾。計算結果同愛荷華大學的實驗進行對比，吻合良好。

1 數學模型

1.1 流體運動方程

連續性方程：

動量方程：

對湍流的模擬采用雷諾平均法可得：

式中：ui，uj為速度分量；P為流體壓力；μ為動力粘性系數；fi為坐標變換引起的慣性力，其中是隨船坐標系下的線速度和角速度是隨船坐標系下的位置矢量(xs,ys,zs);Ui是斜航特征速度。自由液面采用流體體積函數VOF方法來追蹤自由液面的運動狀態：

式中：aq是某計算單元內第q相流體占的體積分數；ρq是第q相的密度。

文中選用SST k-ω湍流模型來封閉控制方程；微分方程的離散采用有限體積法；對流項采用二階迎風差分格式；湍動能及其耗散、體積分數都采用二階迎風格式；擴散項采用中心差分格式；壓力和速度耦合采用PISO方法。松弛因子除了壓力項都選用0.3，壓力項用0.1。時間積分采用二階歐拉隱士格式，時間步長取為0.002 s，內部迭代10次。

1.2 船舶運動方程

云南電力市場中，對具有年調節能力以上的機組的調節電量有一定的政策傾斜。具有年調節能力及以上水庫的水電廠(小灣、糯扎渡、龍江、馬鹿塘、普西橋、泗南江、小中甸)的調節電量(年設計電量的25%，枯水期每月分配年調節電量的10%，豐水期每月分配調節電量的6%)按照核定的上網電價。且接入500 kV電網的大型水電機組參與西電東送摘牌，價格在0.24～0.25元/kWh之間，高于市場電價。即云南已對省內有調節能力和大型的水電機組進行了一定的政策傾斜。

根據剛體六自由運動方程理論，推知船舶的六自由度運動包括重心的平動和繞重心的轉動。船舶的線加速度和角加速度可以分別通過動量定理和動量矩定理求得，如下：

2 算例描述

本文選取系列60船模作為算例。國際拖曳水池會議(ITTC)在1996年把系列60船列為CFD計算船舶阻力和推進的標準船型。它的主尺度如表1所示。采用雙重坐標系來描述斜航船舶的浮態，即慣性坐標系和隨船坐標系，如圖1所示。隨船坐標系的原點位于船模重心位置，x軸指向船首方向，y軸指向左舷，z軸指向垂直甲板方向向上。歐拉角φ、θ、ψ分別表示橫傾角、縱傾角和艏搖角，規定向右舷側轉向的橫傾角為正，抬艏的縱傾角為正，升沉量為正表示船舶上浮，斜航角β以偏向右舷側為正。

圖1 坐標系Fig.1 Coordinate system

表1 系列60實船、船模主要參數Tab.1 Parameters of series 60 ship and model

3 數值計算

3.1 計算方法

要計算出斜航船的浮態變化，這需要耦合求解船舶的六自由度運動方程和流體流動方程。大致步驟如下：

1）RANS求解器求解船舶周圍流場，通過積分船表面法向和切向應力得到作用在船上的力和力矩；

2）把第一步中求得力和力矩輸入到六自由度求解器中計算出船的加速度、速度和位移；

3）更新船舶浮態，再一次重復第一、二步，直至浮態不再發生變化。

3.2 計算域和邊界條件

計算域及其邊界條件設置如圖2所示。入流邊界位于距離船首1.5倍船長位置；左舷側的遠場邊界距離左舷側1.5倍船長；右舷側的出流邊界距離右舷側2.5倍船長；出流邊界距離船尾2.5倍船長，其上設置壓力出口邊界條件，且P=0。船表面設置為不可穿透無滑移邊界條件，為了避免甲板上浪，干舷被加高1.5倍吃水，上表面距離加高后的甲板1倍船長；文中模擬無限水域中的斜航情況，水底距離水面1.5倍船長。在入流邊界、遠場邊界、上表面和水底設置為速度進口邊界條件，且速度為u=-Ucosβ，v=-Usinβ，w=0，式中，u，v，w為速度分量，U為斜航特征速度。

圖2 計算域及其邊界條件Fig.2 The computational domain and boundary conditions

3.3 計算網格

重疊網格擁有網格邏輯關系簡單，流場計算精度高、效率高、壁面粘性模擬能力強等優點，更彌補了結構網格[10]對外形適應能力差的缺點，相比較于傳統的網格生成方法，重疊網格更易于生成高質量網格，對大幅度運動的斜航船浮態求解精度較高。本文中采用重疊網格來劃分流場，網格劃分如圖3所示，用重疊網格將流動區域分為兩個子區域，即背景網格區域和重疊網格區域，兩個子區域內網獨立生成，彼此存在重疊關系，在重疊區域內流場信息通過插值進行數據交換。重疊網格可以在背景網格中作大幅運動而不會出現網格變形和重生成，保證了斜航船大幅度運動時的網格質量。圖3（a）中貼近船舶周圍最內層顏色最深長方形部分為重疊網格區域，網格密度從里到外由高到低平滑過渡。在船舶斜航中，船兩側的波浪是不對稱的，文中在自由液面附近，采用兩邊不對稱的四面體加密區域（如圖3（a）所示）來捕捉船行波。此外，船附近生成邊界層網格以更加準確地求解流場。

圖3 網格劃分Fig.3 Grid generation

4 計算結果分析與討論

圖4，圖5分別展示了Fr=0.316,β=10°時升沉、縱傾、橫傾時歷曲線和總阻力系數（CT）、側向力系數(CS)、艏搖力矩(CM)系數時歷曲線，其中橫坐標均為時間（s），縱坐標為無量綱數值。為了同愛荷華大學的試驗數據進行直接對比，將總阻力、側向力和首搖力矩進行無量綱化，無量綱化方法如式（9）-（11）。

本文中定義一個比較誤差，具體為：E%D=100*(S-D)/D，式中，S為計算值，D為試驗值。從圖4看出除了橫傾角以外，其他考察的物理量都能在很短的時間內收斂到一個穩定值。橫傾角以一種衰減的形式逐步達到穩定。橫傾角、縱傾角，升沉量都為負，這表明斜航船處于向左舷側（迎風面）橫傾、埋首和下沉的狀態。橫傾角、縱傾角、升沉、總阻力系數、側向力系數和艏搖力矩系數同愛荷華大學試驗數據的比較結果如表2所示，誤差分別為-7.35、+6.8、-8.75、-4.93、-4.64和+4.81，浮態計算值基本控制在±10%以內，總體上可以滿足實際工程需要。

表2 浮態計算結果同實驗數據對比Tab.2 Comparisons between Computational result with experimental data

值得注意的是，側向力計算誤差為4.64%，之前的學者，比如CuraHochbaum[11]、Rodrigo Azcueta[12]在數值模擬斜航問題時，側向力的計算誤差都達到了25%量級，其中CuraHochbaum把側向力誤差偏大歸結于在他的計算中沒有考慮自由液面和沒有計及浮態的影響，本文的側向力計算精度有了很大的改善。

船模試驗中只有不考慮浮態的波高圖可供參考，本文中將不考慮浮態的計算波高、考慮浮態的計算波高分別同不考慮浮態的試驗波高進行對比，如圖6所示。不考慮浮態的計算波高（b）同試驗（a）相比很接近，出現在船首的最大波高同試驗值基本保持一致，出現在肩部的最小波高計算值比試驗值略大。在考慮浮態的后的波形變化，對比試驗波高圖，最大波高抬高了3%，而波谷降低了34%。最大波峰依然出現在船首迎風面，最大波谷出現在船的肩部。波形變化不大，這包括波包角、散波角和尾流場中的楔形角，只是波浪衰減的比較厲害，主要是由于網格離船越遠，密度越低造成的。

圖4 升沉、縱傾、橫傾時歷曲線圖Fig.4 Curves of sinkage、trim and heel

圖5 總阻力系數、側向力系數、艏搖力矩系數時歷曲線Fig.5 Coefficients curves of resistance、lateral force and yawing moment

圖6 計算波高同試驗波高對比Fig.6 Comparisons between Computationalwave height with experimental wave height

船體表面動壓分布情況如圖7所示。通過計算結果可知船首迎風面壓力大于背風面，而船尾兩側壓力變化不明顯，由于兩側壓力差產生了轉首力矩。迎風面底部壓力大于背風面壓力，這導致了船向迎風面橫傾。橫傾的同時也出現了首傾現象，這導致了船首底部壓力大于船尾壓力的情況。

圖7 船體表面動壓分布Fig.7 Dynamic pressure contours at hull surface

圖8展示了船舯橫剖面上Z向速度計算同試驗結果的對比圖。圖8(a)可以明顯地看出船有一個向左舷側橫傾的狀態，這造成了船舭部形成了明顯的渦。同試驗對比可以看出計算的負渦核心筒試驗比較接近，但是正的渦核心位置差別較大，試驗中的正渦核心遠離船體，出現流體分離現象，而計算中的正渦核心也有遠離船體的趨勢但還沒有遠離。這也可能是造成側向力以及升沉誤差較大的原因。此外，最大的速度偏低了11%，最小速度偏高了5%。總體而言，采用本文中的重疊網格技術可以較為準確地預報斜航船舶周圍流場。

圖8 計算船舯橫剖面速度云圖同試驗對比Fig.8 Amidships cross-sectional velocity contours compared with experiment data

5 結語

本文給出了一種基于CFD理論進行斜航船舶浮態數值預報的方法。采用重疊網格技術，通過耦合求解流體運動方程和船舶六自由度方程的方法，計算出了系列60船在傅汝德數Fr=0.316和斜航角β=10°情況下的斜航船浮態，包括升沉、縱傾和橫傾等。計算顯示采用重疊網格能很好的處理斜航船大幅運動問題，計算中網格不變形和重生成充分保障了網格質量，而這些都保證了計算精度。將計算結果同愛荷華大學的試驗數據進行對比，偏差基本控制在10%以內，吻合較好，體現出該方法具有較強的的工程實用性。文中展示了部分斜航船周圍的流場信息包括波高圖和舯橫剖面速度云圖，計算結果較為準確的預報了斜航船周圍流場變化，對理解斜航船浮態變化機理有一定的幫助。

[1]郭春雨,趙慶新,趙大剛.基于CFD仿真模擬的船模自航試驗數據處理[J].船海工程,2013,42(3):17-20.

[2]OKUNO T,TANAKA N,HASEGAWA Y.Flow field measurements around ship hull at incidence[J].J Kansai Soc Naval Arch, 1989,212:67-74.

[3]LONGO J,STERN F.Effects of drift angle on model ship flow[J].Experiments in Fluids,2002,32:558-569.

[4]INOUE S,HIRANO M,KIJIMA K.Hydrodynamic derivatives on ship manoeuvring[J].International Shipbuilding Progress,198, 28(321):112-125.

[5]YUKAWA K,KIJIMA K.An estimation of hydrodynamic forces acting on a ship hull in manoeuvring motion[J].Trans West-Japan Soc Nav Arch,1998,95:67-79.

[6]TOXOPEUS S,SIMONSEN C D,GUILMINEAU E,et al.Viscous-flow calculations for KVLCC2 in manoeuvring motion in deep and shallow water[C]//NATO AVT-189 Specialists Meeting on Assessment of Stability and Control Prediction Methods for NATO Air and Sea Vehicles.Portsdown West,UK,2011:151-169.

[7]WANG H M,ZOU Z J,TIAN X M.Computation of the viscous hydrodynamic forces on a KVLCC2 model moving obliquely in shallow water[J].Journal of Shanghai Jiaotong University,2009,14(2):241-244.

[8]劉 晗,馬 寧,鄧德衡,等.循環水槽在平面運動機構試驗中的應用及其數值驗證 [C]//2013年船舶水動力學學術會議, 2013:623-631.

[9]鄒 璐.淺水中低速斜航運動船舶水動力預報及誤差分析[C]//2013年船舶水動力學學術會議,2013:503-510.

[10]王留洋,陳曉亮.甌江口波流耦合作用下二維流場數值模擬研究[J].船海工程,2013,42(5):130-137.

[11]CURAHOCHBAUM A.Computation of the turbulent flow around a ship model in steady turn and in steady oblique motion[C] //Proc.22nd Symp on Naval Hydro,Washington,D.C.,USA,1998:198-213.

[12]AZCUETA R.Computation of turbulent free-surface flows around ships and floating bodies[J].Ship Technology Research, 2002,49(2):46-69.

CFD Theory Based Floating State Prediction of Oblique Ship Motion

SHAO Jiangming1，YING Yeju1，ZHOU Qi2,SUN Xiangjie3
（1.School of Naval Architecture&Ocean Engineering,Zhejiang Ocean University,Zhoushan,316022；2.Zhoushan Wonderful Marine Design CO.LTD.，Zhoushan 316021；3.School of Naval Architecture& Ocean Engineering，Harbin Engineering University，Harbin 150001，China）

The oblique ship motion calculation is rarelyconsidering the effects of ship floating state.However,the change of floating state can not be ignored due to it’s great effect on numerical calculation precision. In this paper,the oblique motion of Series 60 ship model with Fr=0.316,yaw angleβ=10°is taken as example. The running attitudes including the sinkage,trim and heel are obtained by coupling the 6-DOF equations of ship motion with the Reynolds Averaged Navier-Stokes(RANS)equations.The overset grid method is adopted to deal with the sinkage,trim and heel in large amplitude.Comparisons between the numerical results with IWOA (The University of Iowa)experimental data are made,and it is observed that the prediction error can be controlled bellow 10%,which successfully validates the applicability and efficiency of the proposed method in engineering applications.

Oblique Motion；Running Attitude；Six-DOF Equation；RANS

U661.3

A

1008－830X(2015)06－0559－06

2015-08-20

邵江明（1990-），男，碩士研究生，主要從事船舶水動力研究和船型優化設計.