邵江明，應(yīng)業(yè)炬

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

基于CFD的斜航船舶浮態(tài)數(shù)值預(yù)報(bào)

邵江明，應(yīng)業(yè)炬

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

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

斜航；浮態(tài)計(jì)算；重疊網(wǎng)格；數(shù)值模擬

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

目前對(duì)船舶斜航的研究主要有約束模試驗(yàn)[2-3]、基于勢(shì)流理論的理論計(jì)算[4-5]和和基于黏性流理論數(shù)值模擬[6-7]。約束模試驗(yàn)方法具有一定的精度優(yōu)勢(shì)，但一般費(fèi)時(shí)費(fèi)力?；趧?shì)流理論的計(jì)算對(duì)受黏性影響較大的運(yùn)動(dòng)很難模擬準(zhǔn)確，比如橫傾運(yùn)動(dòng)。近些年來，隨著計(jì)算機(jī)科學(xué)技術(shù)的飛速發(fā)展，基于計(jì)算流體動(dòng)力學(xué)（CFD）的數(shù)值模擬方法得到越來越廣泛的應(yīng)用，模擬中引入湍流模型概念充分考慮到流體的黏性作用。目前，這種方法已成為與約束模試驗(yàn)并駕齊驅(qū)、相輔相成的一種水動(dòng)力預(yù)報(bào)方法?；仡櫼酝鶎?duì)船舶斜航的數(shù)值模擬，可以發(fā)現(xiàn)，在絕大部分的研究中很少在模擬中考慮到船舶浮態(tài)影響。然而，浮態(tài)對(duì)水動(dòng)力計(jì)算精度的影響是不容忽視的。劉晗，馬寧等[8]在對(duì)比由Fluent軟件計(jì)算得到的斜航KVLCC2船模水動(dòng)力的數(shù)值結(jié)果和由循環(huán)水槽平面運(yùn)動(dòng)機(jī)構(gòu)試驗(yàn)得到的試驗(yàn)結(jié)果中發(fā)現(xiàn)，在斜航角為20度時(shí)，不考慮浮態(tài)情況下的水動(dòng)力系數(shù)同試驗(yàn)值誤差達(dá)到20%之多。鄒璐[9]在研究數(shù)值計(jì)算精度中提到，對(duì)比考慮斜航船舶升沉、縱傾后的水動(dòng)力計(jì)算值和不考慮浮態(tài)改變的情況，誤差在20%-25%之間。

在本文的數(shù)值模擬中，采用Star-ccm軟件，通過耦合求解船舶的運(yùn)動(dòng)方程和流體運(yùn)動(dòng)方程，計(jì)算出了斜航船的浮態(tài)變化，其中包括升沉、縱傾以及橫傾。計(jì)算結(jié)果同愛荷華大學(xué)的實(shí)驗(yàn)進(jìn)行對(duì)比，吻合良好。

1 數(shù)學(xué)模型

1.1 流體運(yùn)動(dòng)方程

連續(xù)性方程：

動(dòng)量方程：

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

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

式中：aq是某計(jì)算單元內(nèi)第q相流體占的體積分?jǐn)?shù)；ρq是第q相的密度。

文中選用SST k-ω湍流模型來封閉控制方程；微分方程的離散采用有限體積法；對(duì)流項(xiàng)采用二階迎風(fēng)差分格式；湍動(dòng)能及其耗散、體積分?jǐn)?shù)都采用二階迎風(fēng)格式；擴(kuò)散項(xiàng)采用中心差分格式；壓力和速度耦合采用PISO方法。松弛因子除了壓力項(xiàng)都選用0.3，壓力項(xiàng)用0.1。時(shí)間積分采用二階歐拉隱士格式，時(shí)間步長(zhǎng)取為0.002 s，內(nèi)部迭代10次。

1.2 船舶運(yùn)動(dòng)方程

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

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

2 算例描述

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

圖1 坐標(biāo)系Fig.1 Coordinate system

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

3 數(shù)值計(jì)算

3.1 計(jì)算方法

要計(jì)算出斜航船的浮態(tài)變化，這需要耦合求解船舶的六自由度運(yùn)動(dòng)方程和流體流動(dòng)方程。大致步驟如下：

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

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

3）更新船舶浮態(tài)，再一次重復(fù)第一、二步，直至浮態(tài)不再發(fā)生變化。

3.2 計(jì)算域和邊界條件

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

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

3.3 計(jì)算網(wǎng)格

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

圖3 網(wǎng)格劃分Fig.3 Grid generation

4 計(jì)算結(jié)果分析與討論

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

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

表2 浮態(tài)計(jì)算結(jié)果同實(shí)驗(yàn)數(shù)據(jù)對(duì)比Tab.2 Comparisons between Computational result with experimental data

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

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

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

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

圖6 計(jì)算波高同試驗(yàn)波高對(duì)比Fig.6 Comparisons between Computationalwave height with experimental wave height

船體表面動(dòng)壓分布情況如圖7所示。通過計(jì)算結(jié)果可知船首迎風(fēng)面壓力大于背風(fēng)面，而船尾兩側(cè)壓力變化不明顯，由于兩側(cè)壓力差產(chǎn)生了轉(zhuǎn)首力矩。迎風(fēng)面底部壓力大于背風(fēng)面壓力，這導(dǎo)致了船向迎風(fēng)面橫傾。橫傾的同時(shí)也出現(xiàn)了首傾現(xiàn)象，這導(dǎo)致了船首底部壓力大于船尾壓力的情況。

圖7 船體表面動(dòng)壓分布Fig.7 Dynamic pressure contours at hull surface

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

圖8 計(jì)算船舯橫剖面速度云圖同試驗(yàn)對(duì)比Fig.8 Amidships cross-sectional velocity contours compared with experiment data

5 結(jié)語(yǔ)

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

[1]郭春雨,趙慶新,趙大剛.基于CFD仿真模擬的船模自航試驗(yàn)數(shù)據(jù)處理[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]劉 晗,馬 寧,鄧德衡,等.循環(huán)水槽在平面運(yùn)動(dòng)機(jī)構(gòu)試驗(yàn)中的應(yīng)用及其數(shù)值驗(yàn)證 [C]//2013年船舶水動(dòng)力學(xué)學(xué)術(shù)會(huì)議, 2013:623-631.

[9]鄒 璐.淺水中低速斜航運(yùn)動(dòng)船舶水動(dòng)力預(yù)報(bào)及誤差分析[C]//2013年船舶水動(dòng)力學(xué)學(xué)術(shù)會(huì)議,2013:503-510.

[10]王留洋,陳曉亮.甌江口波流耦合作用下二維流場(chǎng)數(shù)值模擬研究[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-），男，碩士研究生，主要從事船舶水動(dòng)力研究和船型優(yōu)化設(shè)計(jì).