散貨船縱傾減阻及其成因分析

宋磊，童駿，孔斌

(1. 華中科技大學(xué)船舶與海洋工程學(xué)院，湖北 武漢 430074；2. 高新船舶與深海開發(fā)裝備協(xié)同創(chuàng)新中心，上海 200240；3. 船舶和海洋工程水動力湖北省重點(diǎn)實(shí)驗(yàn)室，湖北 武漢 430074；4. 武漢第二船舶設(shè)計(jì)研究所，湖北 武漢 430064)

0 引 言

全球氣候變暖最重要的一個因素就是人類活動造成的溫室氣體（GHG）排放，在所有的溫室氣體排放量中，航運(yùn)業(yè)占了相當(dāng)大的比例，其中CO2占2.9%，NOx占14%～15%，SOx占16%[1]。2018年4月國際海事組織（IMO）海上環(huán)境保護(hù)委員會第72次會議（MEPC 72）通過了關(guān)于減少船舶溫室氣體排放的“巴黎協(xié)定”，該協(xié)定要求到2030年，每個運(yùn)輸單位的二氧化碳排放量至少較2008年降低40%，并爭取在2050年實(shí)現(xiàn)70%的降低。對于船舶制造和營運(yùn)而言，如何降低溫室氣體排放已經(jīng)成為一個不得不考慮的問題。同時，隨著金融危機(jī)對航運(yùn)業(yè)的影響，運(yùn)價下跌、貨量減少，船東們?yōu)榱四軌颢@得利潤必須想方設(shè)法提高營運(yùn)效率，降低成本?？v傾調(diào)節(jié)技術(shù)在不改變船舶載重和航速的前提下，僅改變船舶航行過程中的縱傾角，進(jìn)而改變其航行阻力，能夠起到節(jié)能減排的作用。與其他節(jié)能減排方式相比，縱傾調(diào)節(jié)技術(shù)具有易于實(shí)現(xiàn)、成本低廉、效果明顯等優(yōu)點(diǎn)[2]。

隨著計(jì)算機(jī)技術(shù)的發(fā)展，計(jì)算流體力學(xué)（CFD）也得到了飛速的發(fā)展。相比模型試驗(yàn)，數(shù)值模擬具有費(fèi)用低、周期短、無尺度效應(yīng)等優(yōu)點(diǎn)，且無需考慮船模加工精度和船模變形的影響，并可獲得較為詳細(xì)的流場信息。徐杰等[3-4]采用直接求解RANS方程的方法對20艘散貨船的繞流場和阻力進(jìn)行了數(shù)值模擬，通過與試驗(yàn)結(jié)果的詳細(xì)對比，得出結(jié)論：只要選取合適的湍流模型，數(shù)值模擬結(jié)果能夠達(dá)到較高的精度,能滿足工程應(yīng)用。余建偉等[5-6]就漁政船在各航速下船體阻力進(jìn)行了理論計(jì)算和船模試驗(yàn)驗(yàn)證研究，對理論計(jì)算中的船體壓阻力、黏性阻力等數(shù)據(jù)分析，并提出分成分計(jì)算來簡化船舶阻力計(jì)算的思路。

縱傾優(yōu)化研究方面，涂海文等與中國船級社開發(fā)一套縱傾優(yōu)化軟件并用于指導(dǎo)實(shí)船壓載，通過對比得出4%～6%的節(jié)能效果；Salma Sherbaz等[7]通過對MOERI Container Ship進(jìn)行數(shù)值仿真，并與試驗(yàn)結(jié)果比較，分析了集裝箱船在縱傾調(diào)節(jié)過程中阻力各成分的變化趨勢。毛文雷等[8]應(yīng)用Rankine源自由面勢流理論計(jì)算不同縱傾下船舶興波阻力系數(shù)以及總阻力系數(shù)，應(yīng)用三維面元法計(jì)算規(guī)則波中的波浪增阻，并針對實(shí)海域進(jìn)行了増阻預(yù)報(bào)，得出一定的尾傾能減小興波阻力和總阻力的結(jié)論。本文運(yùn)用計(jì)算流體力學(xué)方法求解180000DWT散貨船在設(shè)計(jì)狀態(tài)不同縱傾角下阻力值和流場信息，通過模型試驗(yàn)驗(yàn)證阻力計(jì)算的準(zhǔn)確性，并比較同縱傾狀態(tài)下阻力值，得到船舶設(shè)計(jì)狀態(tài)的阻力最佳縱傾。進(jìn)一步計(jì)算出模型對應(yīng)狀態(tài)興波阻力，通過各個阻力成分流場細(xì)節(jié)的分析和比較，得到了船舶縱傾優(yōu)化過程中阻力變化的主導(dǎo)因素。

1 研究對象

文中計(jì)算模型選用180 000 t散貨船，該船具有肥大型球鼻首，實(shí)船設(shè)計(jì)吃水為16.5 m，設(shè)計(jì)航速v=15 kn，此狀態(tài)Froude數(shù)為Fn=0.1447。按照縮尺比λ=55進(jìn)行換算后各尺寸如表1所示。

表 1 船舶主要參數(shù)Tab. 1The main characteristics of ship

平浮時首尾對應(yīng)位置如圖1所示。在研究的過程中，保證排水量不變的前提下分別對船舶進(jìn)行不同的縱傾角調(diào)整，結(jié)合實(shí)船航行工況，分別為首傾0.804°，0.603°，0.402°，0.201°；尾傾 0.251°，0.503°，0.754°，1.005°。

圖 1 平浮時設(shè)計(jì)吃水位置Fig. 1The place of design draft in even keel condition

數(shù)值計(jì)算三維模型選用UG建模，如圖2所示。

圖 2 計(jì)算模型Fig. 2Calculation model

2 控制方程

采用RANS控制方程組，包括連續(xù)性方程和動量守恒方程[9]，分別如下：

連續(xù)性方程

動量守恒方程

式中：ui為流體時均速度分量；p為流體壓強(qiáng)；fi為流體體積力分量；ρ為流體密度；μ為流體的粘性系數(shù)；附加應(yīng)力記為稱為雷諾應(yīng)力。

針對肥大型球鼻艏船舶，選用SST κ-ω湍流模型封閉方程組數(shù)值預(yù)報(bào)都能達(dá)到比較高的精度[10-11]。其方程如下：

雷諾應(yīng)力的渦粘性模型

計(jì)算中兩相流模型選用VOF模型方程[12]

3 研究方法

考慮到船舶對稱性，建立半船模型。數(shù)值模擬計(jì)算域船長方向5倍船長，其中船前方取1倍船長，船后方約3倍船長；船寬方向取2倍船長，水深方向取1倍船長。如圖3所示。

圖 3 計(jì)算域示意圖Fig. 3Computational domain

計(jì)算域的離散通過ICEM網(wǎng)格劃分工具來完成，全局采用六面體結(jié)構(gòu)網(wǎng)格，網(wǎng)格結(jié)構(gòu)類型采用H-C型網(wǎng)格，即縱向?yàn)镠型網(wǎng)格，橫向?yàn)镃型網(wǎng)格。由于壁面附近的流動船體表面曲線變化較大，在靠近船體處網(wǎng)格適當(dāng)加密；自由面需要捕捉波形，靠近自由面附近網(wǎng)格較密，流場底部網(wǎng)格較稀。網(wǎng)格總數(shù)200萬左右，y+值在80左右。船體表面網(wǎng)格如圖4所示。

圖 4 船體表面網(wǎng)格Fig. 4The mesh of ship

計(jì)算過程中邊界條件設(shè)置如下：計(jì)算域前方為速度入口并設(shè)定來流速度，其中空氣域湍流強(qiáng)度為0.1%，湍動粘度比為1%；水域湍流強(qiáng)度為1%，湍動粘度比為1%。計(jì)算域后方設(shè)為壓力出口，計(jì)算區(qū)域的上方設(shè)為速度入口。底面與側(cè)面設(shè)定為壁面邊界條件，船體表面為壁面邊界條件，中間為對稱面邊界條件。

采用k-ω湍流模型封閉方程組，速度壓力的耦合方法為SIMPLE方法。壓力的插值方法采用體積力加權(quán)方式。動量項(xiàng)和體積分?jǐn)?shù)項(xiàng)均采用2階迎風(fēng)差分格式。自由表面的求解選用了歐拉隱式VOF法[12]。

4 阻力計(jì)算及試驗(yàn)驗(yàn)證

按照上述設(shè)置求解各縱傾角下船舶阻力，計(jì)算100 s左右趨于穩(wěn)定，穩(wěn)定后殘差收斂曲線和阻力收斂曲線如圖5和圖6所示。

模型試驗(yàn)選擇華中科技大學(xué)船模試驗(yàn)水池完成，該水池長175 m，寬6 m，水深4 m，模型如圖7所示。

保證船舶排水量相等前提下對船舶進(jìn)行不同的角度調(diào)整，分別為首傾 0.804°，0.603°，0.402°，0.201°；尾傾 0.251°，0.503°，0.754°，1.005°，為方便表述，將首傾角度定義為負(fù)，尾傾角度定義為正。將不同傾角下模型試驗(yàn)所得數(shù)據(jù)與計(jì)算數(shù)據(jù)對比如表2及圖8所示。

圖 5 殘差收斂曲線Fig. 5Residual convergence curve

圖 6 阻力收斂曲線Fig. 6Resistance convergence curve

圖 7 船舶模型Fig. 7Ship model

將計(jì)算值與試驗(yàn)值比較，可得計(jì)算值與實(shí)驗(yàn)值吻合度較高，最大誤差為3.034%，并且阻力隨傾角變化趨勢與試驗(yàn)一致。通過總阻力的比較，證明選用上述網(wǎng)格以及邊界條件的合理性。

表 2 計(jì)算數(shù)據(jù)與模型試驗(yàn)數(shù)據(jù)對比Tab. 2Comparison of calculations and model test data

圖 8 計(jì)算數(shù)據(jù)和模型試驗(yàn)數(shù)據(jù)Fig. 8Data of calculations and model test

通過對CFD仿真結(jié)果分析可知，船舶在尾傾1.005°時阻力值最小。船舶從平浮向首傾調(diào)整過程中，阻力隨著首傾角的增大不斷增大，在首傾角為-0.804°時達(dá)到最大，相對平浮增大3.353%；在平浮向尾傾調(diào)整過程中，阻力隨著縱傾角增大而減小，在尾傾角為1.005°時效果最好，相對平浮減小5.449%。

5 阻力成分分析

通過上述分析可知，船舶在縱傾調(diào)節(jié)過程中阻力變化明顯，最大時與平浮相差5.449%。為了研究船舶阻力變化過程中各阻力成分的變化情況，采用同樣的網(wǎng)格和邊界條件選用無粘流模型計(jì)算船舶的興波阻力，計(jì)算所得總阻力減去興波阻力即為船舶的粘性阻力。將所得結(jié)果如表3及圖9所示。為了得到船舶縱傾優(yōu)化傾角變化過程中阻力變化的主導(dǎo)因素，將結(jié)果表述為不同傾角下總阻力變化量中分阻力變化量所占的比例：

表 3 不同傾角阻力成分比較Tab. 3Comparison of resistance component on different trim condition

圖 9 不同傾角阻力成分比較Fig. 9Comparison of resistance component on different trim condition

上述表述中分阻力指興波阻力或粘性阻力，計(jì)算所得結(jié)果如表3所示。各分阻力變化量所占比例如圖10所示。

圖 10 各阻力成分變化量比Fig. 10Change of different resistancecomponent on trim

分析數(shù)據(jù)可以看出，在模型總阻力不斷變?。ㄓ墒變A向尾傾變化）過程中，粘性阻力變化占比越來越大，可以說明本散貨船縱傾造成阻力變化主要是由于粘性阻力發(fā)生較大變化引起。本船吃水已超過球鼻艏以及尾封板底部，在縱傾過程中濕面積變化不明顯，因此可以推測粘性阻力變化可能是由于縱傾造成模型水下形狀發(fā)生變化，引起粘壓阻力變化造成。

從Fluent中導(dǎo)出不同傾角計(jì)算穩(wěn)定后波高沿船長方向分布數(shù)值，如圖11和圖12所示。圖中船首為橫坐標(biāo)軸原點(diǎn)處，船尾沿著橫坐標(biāo)增大方向。比較計(jì)算所得興波阻力值與波高圖，可得其大小趨勢吻合。

圖 11 不同首傾波高變化Fig. 11Wave profile of trim by bow

圖 12 不同尾傾波高變化Fig. 12Wave profile of trim by stern

計(jì)算穩(wěn)定后各狀態(tài)興波云圖如圖13～圖21所示。

6 結(jié) 語

圖 13 θ=-0.804°興波云圖Fig. 13Wave making contour of θ=-0.804°

圖 14 θ=-0.603°興波云圖Fig. 14Wave making contour of θ=-0.603°

圖 15 θ=-0.402°興波云圖Fig. 15Wave making contour of θ=-0.402°

圖 16 θ=-0.201°興波云圖Fig. 16Wave making contour of θ=-0.201°

圖 17 θ=0°興波云圖Fig. 17Wave making contour of θ=0°

本文主要研究散貨船縱傾優(yōu)化過程中阻力的變化情況，以及總阻力的變化量中不同阻力成分所占比例。通過計(jì)算阻力值與試驗(yàn)值以及變化趨勢的對比，驗(yàn)證網(wǎng)格和邊界條件的合理性，采用相同的網(wǎng)格計(jì)算興波阻力，總阻力中減去興波阻力余下部分為粘性阻力。通過比較得出如下結(jié)論：

1）比較不同縱傾狀態(tài)總阻力可知，對于本散貨船設(shè)計(jì)吃水而言，模型尾傾狀態(tài)有明顯減阻效果，且在一定范圍內(nèi)隨著尾傾角增大減阻效果加強(qiáng)。

2）通過不同阻力成分的比較可知，該散貨船在縱傾優(yōu)化過程中，總阻力的改變量中粘性阻力占很大比例，興波阻力變化量較小?？傻冕槍Ρ旧⒇洿?，縱傾優(yōu)化過程中粘性阻力的改變?yōu)榭傋枇Ω淖兊闹鲗?dǎo)因素。

圖 18 θ=0.251°興波云圖Fig. 18Wave making contour of θ=0.251°

圖 19 θ=0.503°興波云圖Fig. 19Wave making contour of θ=0.503°

圖 20 θ=0.503°興波云圖Fig. 20Wave making contour of θ=0.503°

圖 21 θ=1.005°興波云圖Fig. 21Wave making contour of θ=1.005°