方昭昭，趙丙乾，陳慶任

1 中國艦船研究設(shè)計中心，湖北武漢430064

2 中國船級社武漢規(guī)范研究所，湖北武漢430022

0 引 言

國際海事組織（IMO）于2011年7月通過了MARPOL 附則VI 有關(guān)船舶能效規(guī)則的修正案，確立了船舶能效指數(shù)（EEDI）及船舶效能管理計劃（SEEMP）等新要求，并將其分階段強制實施［1-2］。EEDI 等新要求的實施，必將促進海事界向節(jié)能減排及綠色環(huán)保方向發(fā)展［3-4］。雙艉船型作為一種優(yōu)秀的節(jié)能船型，在當(dāng)今世界范圍內(nèi)取得了良好的經(jīng)濟效益，特別是在我國內(nèi)河和沿海流域，以優(yōu)良的快速性和操縱性在內(nèi)河和江海直達運輸船型中獨占鰲頭，應(yīng)用頗為廣泛。

按照MARPOL 附則和船舶能效設(shè)計驗證指南對EEDI 檢驗的要求，在設(shè)計階段，船東或造船廠應(yīng)提供船舶滿載或試航狀態(tài)下的功率—航速估算曲線［5］。因此，要準(zhǔn)確計算船舶設(shè)計能效就必須準(zhǔn)確預(yù)報功率—航速曲線，而其中的關(guān)鍵就是阻力的準(zhǔn)確預(yù)報［6］。

近年來，隨著計算技術(shù)的飛速發(fā)展和計算數(shù)學(xué)理論的不斷完善，計算流體動力學(xué)（CFD）得到了蓬勃發(fā)展，成為了船舶水動力學(xué)性能分析的重要手段之一。基于CFD 的數(shù)值預(yù)報因具有費用低、無觸點流場測量、無尺度效應(yīng)、能消除物理模型試驗中由傳感器尺寸及模型變形等因素對流場的影響、可獲得較為詳細(xì)的流場信息等優(yōu)點而廣受關(guān)注，在船舶阻力性能預(yù)報方面的應(yīng)用也越來越廣泛［7-14］。由于雙艉船型艏艉線型復(fù)雜，船體表面曲率變化大，流場的數(shù)值模擬，尤其是艉部模型的生成及數(shù)值處理較困難，因此迄今為止，有關(guān)該船型基于CFD 的阻力預(yù)報研究較少［15］。

本文將基于計算流體動力學(xué)理論，提出一種雙艉船型的阻力數(shù)值預(yù)報方法，即采用SHIP?FLOW 軟件，基于勢流理論計算興波阻力，并基于粘性流理論計算粘性阻力。通過對某典型的散貨雙艉船與液貨雙艉船在不同航速下的阻力開展數(shù)值預(yù)報，并將數(shù)值結(jié)果與試驗數(shù)據(jù)［16］作比較，以驗證所提方法的正確性。

1 計算原理

1.1 計算區(qū)域劃分

運用SHIPFLOW 軟件進行阻力預(yù)報時，不同的阻力成分對應(yīng)的理論模型和計算原理各不相同。如圖1 所示，該軟件將計算流場劃分為了3 個區(qū)域：勢流區(qū)、邊界層區(qū)、粘性流區(qū)。

圖1 SHIPFLOW 中計算流場區(qū)域的劃分Fig.1 Divisions of the zones in SHIPFLOW

勢流區(qū)域（Zone I）：采用勢流面元法計算模塊計算興波阻力及波形。

邊界層區(qū)域（Zone II）：求解邊界層方程模塊，可計算船體的摩擦阻力。

尾部粘性流區(qū)域（Zone III）：求解RANS 方程的粘性流場計算模塊。該模塊由勢流區(qū)域提供流場邊界條件（如入流邊界），計算摩擦阻力和粘壓阻力等。

1.2 基于勢流方法的興波阻力計算

假定船舶以勻速V 在靜水中直線航行，采用大地固定坐標(biāo)系o-xyz。其中xoy 平面在靜水面上，x 軸正方向與船舶航向一致，y 軸正方向指向右舷，z 軸豎直向下，如圖2 所示。

圖2 坐標(biāo)系Fig.2 The coordinate system

假定流場中的流體為無粘、無旋、不可壓的理想流體，則存在一定常速度勢?，它滿足控制方程：

速度勢? 滿足自由面邊界條件、物面邊界條件以及遠(yuǎn)方輻射條件：

式中：Vs=(V，0，0)；g 為重力加速度；R 為流場中某點離擾動源的距離；n 為物面法向矢量。

上述式（1）～式（4）就構(gòu)成了流體速度勢的定解問題，可采用勢流Rankine 源法求解。即在船體表面和自由面布置源（匯），則流場中任意點P(x，y，z)的速度勢表達為

式中：P(x，y，z) 為場點；Q(x0，y0，z0) 和Q′(x0，-y0，z0)為源（匯）點；σ(Q)為源強分布；r和r′為場點與源（匯）點之間的距離。

將流體速度勢表達式（5）代入式（2）及式（4），求解源點源強分布σ(Q)。再由式（5）可得到流體域任意場點的流體速度勢?(P)，由式（2）可得到船體興波波形。

一旦求解出流體速度勢?，根據(jù)伯努利方程，便可得到流場中的壓力（單位：Pa）分布：

式中，ρ 為流體質(zhì)量密度。

將流體壓力沿船體濕表面進行積分，得到船體所受到的流體作用力與力矩：

式中：r 為原點至船體表面點P(x，y，z)的矢徑；n=(n1，n2，n3)，為點P(x，y，z)處的單位法向矢量。

船體興波阻力為

無因次的興波阻力系數(shù)

式中，S0為船體濕表面積。

1.3 基于RANS 方程的粘性阻力計算

當(dāng)船舶以航速V 在水中勻速航行時，船體受到的流體粘性阻力可通過模擬船體粘性疊模繞流場得到，控制方程為RANS 方程：

式中：Ui為流體質(zhì)點在i方向的速度分量；ν為運動粘性系數(shù)；fi為體積力在i 方向的分量；為雷諾應(yīng)力項。

上述控制方程并不封閉，需添加描述流體湍動特征的湍流模型方程。這里采用標(biāo)準(zhǔn)的k-ε模型湍流輸運方程，湍動能k 方程和湍動耗散率ε 方程分別為：

式中：Gk表示由平均速度梯度產(chǎn)生的湍流動能；Gb為由浮力產(chǎn)生的湍流動能；c1ε，c2ε和c3ε為經(jīng)驗常數(shù)；σk和σε分別為k 和ε 對應(yīng)的Prandtl 數(shù)；Sk和Sε為源項；μ 為流體的動力粘性系數(shù)；uτ為湍流粘性系數(shù)。

對式（11）～式（14）采用有限差分法進行數(shù)值求解，可計算得到船體粘性繞流場的壓力與速度分布。其中：船體表面的壓力沿船體濕表面的積分為船體受到的粘性流體總壓力，該力沿著與船舶航速方向相反的分量為粘壓阻力；船體表面的速度梯度乘以流體粘性系數(shù)為船體表面的流體切應(yīng)力，將其沿船體濕表面進行積分則得到船體受到的流體總摩擦力，該力沿著與船舶航速方向相反的分量即為船體摩擦阻力。

2 算例描述

本文以某散貨雙艉船與液貨雙艉船為例，分別對其船模阻力進行數(shù)值預(yù)報。液貨雙艉船的橫剖面型線圖如圖3 所示，從中可看出，這種船型艏艉部的曲面曲率變化非常復(fù)雜。

圖3 某液貨雙艉船橫剖面型線圖Fig.3 The plans of the liquid cargo ship with twin-skeg

散貨雙艉船船模的縮尺比為15.83，其主要參數(shù)如表1 所示，船模的計算航速如表2 所示。液貨雙艉船船模的縮尺比為16.0，其主要參數(shù)如表3所示，船模的計算航速如表4 所示。

表1 散貨雙艉船的主要參數(shù)Tab.1 Principal dimensions of the bulk cargo ship with twin-skeg

表2 散貨雙艉船船模的計算航速Tab.2 The speeds of the bulk cargo ship model with twin-skeg

表3 液貨雙艉船的主要參數(shù)Tab.3 Principal dimensions of the liquid cargo ship with twin-skeg

表4 液貨雙艉船船模的計算航速Tab.4 The speeds of the liquidcargo ship model with twin-skeg

3 數(shù)值計算

3.1 計算模型

SHIPFLOW 根據(jù)加密的分站型線建立船體表面模型。如通過Catia 和Rhino 等三維建模軟件建立船殼表面，然后導(dǎo)入SHIPFLOW 生成計算所需的型值文件。圖4 所示為散貨雙艉船的全船模型；圖5 所示為艉部局部模型。

圖4 散貨雙艉船數(shù)值模型Fig.4 The configuration of the bulk cargo ship model with twin-skeg

圖5 散貨雙艉船的艉部模型（一半）Fig.5 The half stern sketch of the bulk cargo ship with twin-skeg

3.2 網(wǎng)格劃分

SHIPFLOW 軟件中，是用勢流方法計算興波阻力，面元的分布影響結(jié)果精度；采用粘性流方法計算粘性阻力，體網(wǎng)格的劃分對結(jié)果精度有一定的影響。為了提高計算精度，結(jié)合船型自身特征及航行的特點，興波阻力計算的面網(wǎng)格劃分應(yīng)注意以下幾點：

1）在船體瞬時濕表面、瞬時自由面均需劃分面元網(wǎng)格，并分布速度源或匯。

2）在船體表面上，對于曲率變化較大的艏、艉區(qū)域，面元網(wǎng)格應(yīng)進行加密。

3）在自由面上，面元網(wǎng)格沿船長方向的分布應(yīng)盡量均勻，沿寬度方向船體附近的網(wǎng)格應(yīng)較密。

4）艉封板后的自由面需劃分面元網(wǎng)格。

5）不同航速下，面元網(wǎng)格的數(shù)量不同。一般較低航速（Fn＜0.2）的面元網(wǎng)格數(shù)量比較高航速（Fn＞0.2）的大。

對于粘性阻力計算的體網(wǎng)格劃分，應(yīng)注意以下幾點：

1）只在船體中、后部劃分體網(wǎng)格，網(wǎng)格類型為整體結(jié)構(gòu)化貼體網(wǎng)格。

2）在船體表面進行貼體網(wǎng)格的加密，第1 層網(wǎng)格厚度對應(yīng)的y+值約取為1～2。

3）粘壓阻力的計算精度與艉封板周圍的網(wǎng)格質(zhì)量關(guān)系較大，因此在生成真實、精確的艉封板形狀的基礎(chǔ)上，其周圍應(yīng)劃分加密的貼體結(jié)構(gòu)化網(wǎng)格。

興波阻力計算中面網(wǎng)格的劃分如圖6 所示，粘性流場計算中體網(wǎng)格的劃分如圖7 所示。

3.3 阻力的數(shù)值預(yù)報結(jié)果

3.3.1 散貨雙艉船的阻力計算結(jié)果

1）興波阻力計算結(jié)果。

表5 給出了散貨雙艉船船模在不同航速下興波阻力系數(shù)的計算結(jié)果，其中Cw為興波阻力系數(shù)。

圖6 興波阻力計算中面網(wǎng)格的劃分Fig.6 The mesh for calculating wave-making resistance in potential flow

圖7 粘性流場計算中尾部體網(wǎng)格的劃分Fig.7 The stern mesh for calculating viscosity resistance in viscous flow

表5 散貨雙艉船船模的興波阻力系數(shù)計算結(jié)果Tab.5 The coefficients of wave-making resistance for the bulk cargo ship model with twin-skeg

圖8 散貨雙艉船船模的興波阻力系數(shù)變化曲線Fig.8 Coefficients curve of the wave-making resistance for the bulk cargo ship model with twin-skeg

圖8 給出了興波阻力系數(shù)隨航速的變化曲線。從中可以看出，在該航速段內(nèi)，興波阻力系數(shù)隨著航速的增加而增大。在Fn=0.19 附近，興波阻力系數(shù)快速增加。圖9 給出了該船模在3 個不同航速下的興波波形等高線圖。從中可以看出，數(shù)值計算能較好地反映出船體興波的船艏波系和艉波系，其散波和橫波特征均能得到較好的描述。

圖9 不同航速下散貨雙艉船船模船體興波波形Fig.9 Wave contours for the bulk cargo ship model with twin-skeg advancing at different speeds

2）剩余阻力計算結(jié)果。

表6 給出了散貨雙艉船船模在不同航速下的剩余阻力計算結(jié)果，圖10 給出了剩余阻力系數(shù)計算結(jié)果與試驗數(shù)據(jù)隨航速的變化曲線。

表6 散貨雙艉船船模的剩余阻力系數(shù)計算結(jié)果Tab.6 The coefficients of residual resistance for the bulk cargo ship model with twin-skeg

圖10 散貨雙艉船船模的剩余阻力系數(shù)變化曲線Fig.10 Coefficients curves of residual resistance for the bulk cargo ship model with twin-skeg

3）摩擦阻力計算結(jié)果。

表7 給出了散貨雙艉船船模在不同航速下的摩擦阻力計算結(jié)果，并與1957 ITTC 公式估算結(jié)果進行了比較。從中可以看出，摩擦阻力系數(shù)的數(shù)值計算結(jié)果和1957 ITTC 公式的估算結(jié)果相比略小。其主要原因是：由于雙艉船型長寬比較小，船體表面縱向曲率大，艉部形狀復(fù)雜，一方面，船體表面邊界層流動情況與平板邊界層流動情況不同，在船體曲驟處，特別是較豐滿船的艉部易發(fā)生邊界層分離，產(chǎn)生漩渦，摩擦阻力下降；另一方面，雙艉船型艉部較肥大，去流段短，粘壓阻力較單艉船型大，剩余阻力系數(shù)和二因次換算試驗結(jié)果相比也大。

表7 散貨雙艉船船模的摩擦阻力系數(shù)計算結(jié)果Tab.7 The coefficients of frictional resistance for the bulk cargo ship model with twin-skeg

4）總阻力計算結(jié)果。

表8 給出了散貨雙艉船船模的總阻力計算結(jié)果，圖11 給出了總阻力系數(shù)變化曲線。這里的總阻力為興波阻力與粘性阻力之和，也即摩擦阻力與剩余阻力之和。由表8 與圖11 均可看出，總阻力的數(shù)值計算結(jié)果與試驗數(shù)據(jù)非常接近，總阻力系數(shù)曲線的趨勢也一致，最大誤差不超過2%，表明本文的數(shù)值計算方法具有較高的精度。

表8 散貨雙艉船船模的總阻力系數(shù)計算結(jié)果Tab.8 The coefficients of total resistance for the bulk cargo ship model with twin-skeg

圖11 散貨雙艉船船模的總阻力系數(shù)變化曲線Fig.11 Coefficients curves of total resistance for the bulk cargo ship model with twin-skeg

3.3.2 液貨雙艉船的阻力計算結(jié)果

1）興波阻力計算結(jié)果。

表9 給出了液貨雙艉船船模在不同航速下興波阻力系數(shù)的計算結(jié)果，圖12 給出了興波阻力系數(shù)隨航速的變化曲線。從中可以看出，在該航速段內(nèi)，興波阻力系數(shù)隨著航速的增加而增大，在Fn=0.184 附近，興波阻力系數(shù)快速增加。圖13 給出了該船模在3 個不同航速下的興波波形等高線圖。可以看出，數(shù)值計算能較好地反映出船體興波的船艏波系和艉波系，其散波和橫波特征均能較好地得到描述。

表9 液貨雙艉船船模的興波阻力系數(shù)計算結(jié)果Tab.9 The coefficients of wave-making resistance for the liquid cargo ship model with twin-skeg

圖12 液貨雙艉船船模的興波阻力系數(shù)變化曲線Fig.12 Coefficients curve of the wave-making resistance for the liquid cargo ship model with twin-skeg

圖13 不同航速下液貨雙艉船船模船體興波波形Fig.13 Wave contours for the liquid cargo ship model advancing at different speeds with twin-skeg

2）剩余阻力計算結(jié)果。

表10 給出了液貨雙艉船船模在不同航速下剩余阻力系數(shù)的計算結(jié)果，圖14 給出了剩余阻力系數(shù)計算結(jié)果與試驗數(shù)據(jù)隨航速的變化曲線。

3）摩擦阻力計算結(jié)果。

表11 給出了液貨船船模在不同航速下摩擦阻力系數(shù)的計算結(jié)果及其與1957 ITTC 公式估算結(jié)果的比較。

表10 液貨雙艉船船模的剩余阻力系數(shù)計算結(jié)果Tab.10 The coefficients of residual resistance for the liquid cargo ship model with twin-skeg

圖14 液貨雙艉船船模的剩余阻力系數(shù)變化曲線Fig.14 Coefficients curves of residual resistance for the liquid cargo ship model with twin-skeg

表11 液貨雙艉船船模的摩擦阻力系數(shù)計算結(jié)果Tab.11 The coefficients of frictional resistance for the liquid cargo ship model with twin-skeg

4）總阻力計算結(jié)果。

表12 給出了液貨雙艉船船模的總阻力系數(shù)計算結(jié)果，圖15 給出了總阻力系數(shù)的變化曲線。從表12 與圖15 均可看出，總阻力系數(shù)的數(shù)值計算結(jié)果與試驗數(shù)據(jù)非常接近，總阻力系數(shù)曲線的趨勢也一致，最大誤差為3.66%，表明文中數(shù)值計算方法具有較高的精度。

表12 液貨雙艉船船模的總阻力系數(shù)計算結(jié)果Tab.12 The coefficients of total resistance for the liquid cargo ship model with twin-skeg

圖15 液貨雙艉船船模的總阻力系數(shù)變化曲線Fig.15 Coefficients curves of total resistance for the liquid cargo ship model with twin-skeg

3.3.3 結(jié)果分析

分別比較兩種不同類型雙艉船船模的總阻力數(shù)值結(jié)果與試驗數(shù)據(jù)，可以看出：這兩種船型的總阻力計算結(jié)果與試驗數(shù)據(jù)吻合較好，各航速下的總阻力計算誤差均在4%以內(nèi)，且總阻力隨航速的變化趨勢與試驗結(jié)果也較一致。

就計算結(jié)果與船模阻力試驗二因次換算結(jié)果來看，兩種船型的摩擦阻力系數(shù)數(shù)值計算結(jié)果均比基于平板邊界層理論的ITTC 1957 估算公式的值略小，同時剩余阻力系數(shù)的數(shù)值計算結(jié)果比試驗的二因次換算的結(jié)果大。其主要原因是：由于雙艉船型長寬比較小，船體表面縱向曲率大，艉部形狀復(fù)雜，一方面，船體表面邊界層流動情況與平板邊界層流動情況不同，特別是雙艉型船艉部較豐滿，易發(fā)生邊界層分離；另一方面，雙艉型船艉部較肥大，去流段短，粘壓阻力較單艉船型大，因此計算的剩余阻力系數(shù)比二因次換算試驗的結(jié)果要大。

比較散貨船船模與液貨船船模的阻力數(shù)值計算結(jié)果，發(fā)現(xiàn)這兩種船型的主尺度幾乎相同，且形狀相似，但液貨船的方形系數(shù)較散貨船的大，因此液貨船的水線面系數(shù)較大。在相同傅汝德數(shù)情況下，液貨船的興波阻力比散貨船的大，例如，當(dāng)Fn=0.175 時，液貨船船模的興波阻力系數(shù)為1.286×10-3，粘性阻力系數(shù)為3.914×10-3；當(dāng)Fn=0.174時，散貨船船模的興波阻力系數(shù)為0.758×10-3，粘性阻力系數(shù)為3.822×10-3，粘性阻力系數(shù)非常相近。

4 結(jié) 語

本文給出了一種基于CFD 理論進行雙艉船型阻力數(shù)值預(yù)報的方法。采用非線性勢流方法計算船體興波阻力，基于粘性流方法計算船體粘性阻力。分別對典型散貨雙艉船與液貨雙艉船在不同航速下的阻力進行了數(shù)值預(yù)報，并將數(shù)值計算結(jié)果與試驗數(shù)據(jù)進行了比較。結(jié)果顯示，總阻力誤差均在4%以內(nèi)，且阻力的數(shù)值計算結(jié)果隨航速的變化趨勢與試驗結(jié)果吻合良好。研究表明，該方法計算效率較高、易于實現(xiàn)、經(jīng)濟性較好且預(yù)報精度能滿足工程需要，具有較強的工程實用性。

［1］國際海事組織. 防污公約2011 綜合文本［M］. 北京：人民交通出版社，2012.

［2］李路，芮曉松.論EEDI（能效設(shè)計指數(shù)）的強制實施的合理性［J］.中國造船，2011，52（增刊1）：33-37.LI Lu，RUI Xiaosong. Comments to the compulsory im?plementation of EEDI ［J］. Shipbuilding of China，2011，52（Supp1）：33-37.

［3］周偉新，李百齊，胡瓊，等. 中國船舶科學(xué)研究中心關(guān)于EEDI 的研究進展［J］.中國造船，2011，52（4）：13-22.ZHOU Weixin，LI Baiqi，HU Qiong，et al. EEDI—An important factor of green ship studied in China Ship Scientific Research Center［J］. Shipbuilding of China，2011，52（4）：13-22.

［4］劉飛，林焰，李納，等.我國船舶EEDI 分析研究［J］.中國造船，2012，53（4）：128-136.LIU Fei，LIN Yan，LI Na，et al. Research on EEDI analysis for the ships of China［J］. Shipbuilding of Chi?na，2012，53（4）：128-136.

［5］中國船級社. 綠色船舶規(guī)范［M］. 北京：人民交通出版社，2012.

［6］盛振邦，劉應(yīng)中.船舶原理（上冊）［M］.上海：上海交通大學(xué)出版社，2003.

［7］HUA Z L，XING L H，GU L. Application of modified quick scheme to depth-averaged k-epsilon turbulence model based on unstructured grids［J］. Chinese Journal of Hydrodynamics（Ser.B），2008，20（4）：514-523.

［8］KIM K J. Ship flow calculation and resistance minimi?zation［D］. Gothenburg，Sweden：Chalmers University of Technology，1989.

［9］吳曉蓮. 基于CFD 的船舶球艏/球艉低阻線型研究［D］.哈爾濱：哈爾濱工程大學(xué)，2008.

［10］陳偉，許輝，邱遼原，等. 基于SHIPFLOW 軟件的方尾艦船阻力快速預(yù)報［J］. 中國艦船研究，2012，7（4）：17-22.CHEN Wei，XU Hui，QIU Liaoyuan，et al. Fast re?sistance prediction for ship with transom stern based on SHIPFLOW［J］. Chinese Journal of Ship Re?search，2012，7（4）：17-22.

［11］徐力，陳作鋼.船體艏部水動力性能優(yōu)化［J］.中國艦船研究，2012，7（2）：37-64.XU Li，CHEN Zuogang. Hydrodynamic performance optimization of ship hull's forebody［J］. Chinese Jour?nal of Ship Research，2012，7（2）：37-64.

［12］倪崇本，朱仁傳，繆國平，等. 一種基于CFD 的船舶總阻力預(yù)報方法［J］. 水動力學(xué)研究與進展（A輯），2010，25（5）：579-586.NI Chongben，ZHU Renchuan，MIAO Guoping，et al. A method for ship resistance prediction based on CFD computation［J］. Chinese Journal of Hydrody?namics（Ser.A），2010，25（5）：579-586.

［13］NI S Y. Higher order panel method for potential flows with linear or non-linear free surface boundary condi?tions［D］. Gothenburg，Sweden：Chalmers Universi?ty of Technology，1987.

［14］方昭昭，朱仁傳，繆國平，等. 基于數(shù)值波浪水池的波浪中船舶水動力計算［J］. 水動力學(xué)研究與進展（A 輯），2012，27（5）：515-524.FANG Zhaozhao，ZHU Renchuan，MIAO Guoping，et al. Numerical calculation of hydrodynamic forces for a ship in regular waves based on numerical wave tank［J］. Chinese Journal of Hydrodynamics（Ser. A），2012，27（5）：515-524.

［15］張凈宙.雙尾（鰭）船型線特征及局部變形設(shè)計方法研究［D］.武漢：武漢理工大學(xué)，2007.

［16］張偉.內(nèi)河典型船型船機槳匹配優(yōu)化及推進效率研究報告［R］. 武漢：中國船級社武漢規(guī)范研究所，2012.