宋磊，童駿，孔斌

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

0 引 言

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

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

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

1 研究對象

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

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

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

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

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

圖 2 計算模型Fig. 2Calculation model

2 控制方程

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

連續(xù)性方程

動量守恒方程

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

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

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

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

3 研究方法

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

圖 3 計算域示意圖Fig. 3Computational domain

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

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

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

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

4 阻力計算及試驗驗證

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

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

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

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

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

圖 7 船舶模型Fig. 7Ship model

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

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

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

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

5 阻力成分分析

通過上述分析可知，船舶在縱傾調(diào)節(jié)過程中阻力變化明顯，最大時與平浮相差5.449%。為了研究船舶阻力變化過程中各阻力成分的變化情況，采用同樣的網(wǎng)格和邊界條件選用無粘流模型計算船舶的興波阻力，計算所得總阻力減去興波阻力即為船舶的粘性阻力。將所得結(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

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

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

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

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

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

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

計算穩(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)化過程中阻力的變化情況，以及總阻力的變化量中不同阻力成分所占比例。通過計算阻力值與試驗值以及變化趨勢的對比，驗證網(wǎng)格和邊界條件的合理性，采用相同的網(wǎng)格計算興波阻力，總阻力中減去興波阻力余下部分為粘性阻力。通過比較得出如下結(jié)論：

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

2）通過不同阻力成分的比較可知，該散貨船在縱傾優(yōu)化過程中，總阻力的改變量中粘性阻力占很大比例，興波阻力變化量較小?？傻冕槍Ρ旧⒇洿裕v傾優(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°