倪寶玉，胡 冰，王偉彬，狄少丞，薛彥卓

(1. 哈爾濱工程大學(xué)船舶工程學(xué)院，黑龍江 哈爾濱 150001；2. 上海交通大學(xué)海洋工程國家重點(diǎn)實(shí)驗(yàn)室，上海 200240；3. 中國遠(yuǎn)洋海運(yùn)集團(tuán)有限公司，上海 200127；4. 中遠(yuǎn)海運(yùn)特種運(yùn)輸股份有限公司，廣東 廣州 510623)

0 引 言

北極航道是聯(lián)系歐亞美三大洲最短航線，使用北極航道比繞行南部的蘇伊士或巴拿馬運(yùn)河節(jié)省至少40%航程，縮短傳統(tǒng)航道5千千米以上。極地蘊(yùn)藏豐富的自然資源，是未來重要的能源資源基地。“海洋強(qiáng)國”、“一帶一路”、“冰上絲綢之路”等提出，均表明北極航道的開發(fā)對(duì)我國具有重要的經(jīng)濟(jì)價(jià)值和戰(zhàn)略意義[1-2]。2013年，中國商船“永盛輪”首航北極，圓滿完成我國第1艘商船經(jīng)由北極東北航道到達(dá)歐洲的首創(chuàng)[3]，開啟了我國商船的極地航行之旅。

與敞水航行不同，極地航行的商船將不可避免地遭遇冰載荷。冰載荷與水載荷不同，具有幅值大、強(qiáng)間斷性和強(qiáng)非線性的特點(diǎn)，對(duì)極地航行商船的阻力計(jì)算、主機(jī)功率計(jì)算、操縱性計(jì)算以及結(jié)構(gòu)安全性等均提出巨大挑戰(zhàn)[4]。以主機(jī)功率的計(jì)算為例，以往的計(jì)算主要依賴船級(jí)社冰級(jí)規(guī)范的經(jīng)驗(yàn)公式[5-6]。經(jīng)驗(yàn)公式計(jì)算方法簡便，但研究表明，經(jīng)驗(yàn)公式計(jì)算得到的最小主機(jī)功率往往過于保守[7]。尤其對(duì)于常規(guī)商船配備中低速主機(jī)的船舶，此功率要求往往會(huì)大幅度高于日常使用所要求的經(jīng)濟(jì)值[7]。為此，正如芬蘭-瑞典冰級(jí)規(guī)范[6]中所指出的，對(duì)于某些船型，需要采用更加精確的計(jì)算方法或者基于模型試驗(yàn)來確定船舶主機(jī)功率。

在計(jì)算船舶主機(jī)功率之前，需要確定船舶遭遇的冰阻力。目前對(duì)于冰阻力的計(jì)算方法主要有經(jīng)驗(yàn)方法、試驗(yàn)方法和數(shù)值方法。經(jīng)驗(yàn)公式法[8]通過對(duì)實(shí)船測(cè)量、模型試驗(yàn)以及數(shù)值計(jì)算的結(jié)果進(jìn)行總結(jié)，歸納出經(jīng)驗(yàn)公式以預(yù)報(bào)船舶在冰區(qū)航行的冰載荷，但是經(jīng)驗(yàn)公式在預(yù)測(cè)冰載荷問題上經(jīng)常局限于某種特定的船型或者結(jié)構(gòu)形式。試驗(yàn)方法[9-10]是船舶冰阻力研究中比較直接有效的方法，但是試驗(yàn)需要比較苛刻的試驗(yàn)環(huán)境（如實(shí)測(cè)環(huán)境或者冰水池），同時(shí)需要比較昂貴的成本，使得試驗(yàn)方法具有很大的局限性。數(shù)值方法是通過數(shù)學(xué)建模，將海冰與船體的相互作用采用數(shù)值手段進(jìn)行求解并積分獲得冰阻力的方法。目前計(jì)算冰阻力主流的數(shù)值方法包含有限元法[11-12]（FEM）和離散元法[13-14]（DEM）。由于離散元法將海冰視為塊體或顆粒體的思想更接近海冰真實(shí)的物理力學(xué)特性，近年來離散元法在計(jì)算海冰與船舶結(jié)構(gòu)的碰撞和阻力方面的研究有了很大進(jìn)展[15-16]。

針對(duì)目前極地航行商船主機(jī)功率預(yù)報(bào)中缺乏有效計(jì)算方法與模型的問題，本文借鑒中國商船“永盛輪”首航北極的實(shí)際航行經(jīng)驗(yàn)，對(duì)航行東北航道至北歐地區(qū)的船舶進(jìn)行技術(shù)調(diào)研，基于離散元方法和開敞水域船舶的主機(jī)功率計(jì)算理論，建立冰級(jí)船主機(jī)功率的計(jì)算模型和預(yù)報(bào)方法。將預(yù)報(bào)方法應(yīng)用于中國遠(yuǎn)洋海運(yùn)集團(tuán)有限公司B1冰級(jí)的36 000 t冰級(jí)多用途船的主機(jī)功率預(yù)報(bào)中，計(jì)算結(jié)果與船級(jí)社規(guī)范方法進(jìn)行對(duì)比，總結(jié)相關(guān)規(guī)律，旨在為我國極地航行商船設(shè)計(jì)制造提供技術(shù)支撐。

1 冰阻力數(shù)值模型

1.1 航道內(nèi)冰況確定

本文采用離散元方法計(jì)算船舶遭遇的冰阻力。在計(jì)算冰阻力之前，需要給出船舶在極地航道航行時(shí)遭遇的海冰狀態(tài)，例如是平整冰、碎冰或冰脊甚至冰山等。根據(jù)“永盛輪”首航北極的實(shí)際航行經(jīng)驗(yàn)，我國極地航行商船多航行于破冰船開拓的或者其他商船往復(fù)航行的碎冰航道內(nèi)。碎冰航道內(nèi)碎冰的尺度較小，形狀更接近球形。在此基礎(chǔ)上，本文將建立基于球形顆粒的海冰離散元模型。

1.2 碎冰離散元模型

在采用離散單元模型計(jì)算碎冰間相互作用過程中，將碎冰看作為具有一定質(zhì)量和大小的顆粒單元，考慮單元間相對(duì)速度和相對(duì)位置引起的非線性粘彈性作用力，并采用Mohr-Coulomb摩擦定律確定單元間的剪切力，其接觸力模型[17]如圖1所示。圖1中MA和MB分別為圓球顆粒單元A，B的質(zhì)量，Kn和Ks分別為法向和切向剛度系數(shù)，Cn和Cs分別為法向和切向阻尼系數(shù)，μ為滑動(dòng)摩擦系數(shù)。

圖 1 顆粒單元間的接觸模型Fig. 1 Contact model between elements

在顆粒接觸的法線方向，顆粒單元間的法向力包括Hertz非線性彈性力和非線性粘滯力。彈性力模擬顆粒間相互接觸時(shí)的排斥力，而非線性粘滯力模擬顆粒接觸過程中因相對(duì)速度導(dǎo)致的能量耗散，可表述為：

式中：xn和x˙n分別為顆粒的法向重疊量和相應(yīng)速率；A為顆粒材料的粘滯性參數(shù)。

在顆粒接觸的切線方向，基于Mindlin理論和Mohr-Coulomb摩擦定律，并忽略切向粘滯力影響，則切向接觸力為：

式中,xs為切向變形。法向剛度Kn和切向剛度Ks分別為：

非線性離散元模擬中，計(jì)算時(shí)間步長一般通過由顆粒表面瑞雷波的傳播周期[18]來確定。首先定義臨界時(shí)間步長為[18]：

式中,為 顆粒最小半徑，ρ為材料密度。

在離散元計(jì)算時(shí)，實(shí)際時(shí)間步長要小于臨界時(shí)間步 長 Δtcrit， 可 寫 為 ：

式中,K為經(jīng)驗(yàn)系數(shù)，要求K≥1。通常當(dāng)顆粒材料的配位數(shù)大于4時(shí)，取K=5即 Δt=0.2Δtcrit；當(dāng)配位數(shù)小于4時(shí)，取K=2.5即 Δt=0.4Δtcrit。碎冰航道內(nèi)海冰顆粒為密集排列，其配位數(shù)要大于4，因此本文取K=5。

1.3 海冰與船體耦合作用模型

采用網(wǎng)格劃分前處理軟件將船體離散化為一系列的三角形單元，如圖2所示。建立三角形單元與球體海冰單元間的接觸模型，由此計(jì)算海冰與船體之間的相互作用。這里主要介紹海冰與船體間的作用力模型，具體的接觸判斷算法詳見文獻(xiàn)[19]。

圖 2 船體三角形離散單元模型Fig. 2 Discrete triangle elements on ship surfaces

通過接觸模型可以確定海冰球狀顆粒在船體三角形單元內(nèi)的嵌入量 ΔL，如下式：

式中，P代表顆粒點(diǎn)坐標(biāo)，點(diǎn)PC代 表P對(duì)應(yīng)的接觸點(diǎn)坐標(biāo)[19]，Rball為海冰球狀顆粒的半徑。顆粒與三角形單元接觸點(diǎn)之間的法向量可表示為：

相對(duì)位移Δx可表示為：

式中，vp為 顆粒速度，vw為三角形單元上接觸點(diǎn)的速度。

將相對(duì)位移沿接觸面分解為法向分量 Δxn和切向分量 Δxs， 分 別 為 ：

顆粒與邊界單元之間的相互作用力采用線性接觸模型計(jì)算，其中法向力可表示為：

切向力的計(jì)算采用增量形式，并服從Mohr-Coulomb定律：

式中：為顆粒與邊界單元之間的接觸剛度，為顆粒與邊界單元間的最大摩擦系數(shù)。求得作用在每個(gè)顆粒上的作用力后，即可更新每個(gè)顆粒的速度和位置坐標(biāo)。反之，將法向力與切向力沿著船體表面進(jìn)行積分，即可獲得船體遭受的冰阻力Ri，如下式：

式中，Si為船冰的接觸面積，x為船體運(yùn)動(dòng)方向，nx為接觸面積內(nèi)船體表面單位外法向沿著x方向的分量。

2 摩擦水阻力經(jīng)驗(yàn)公式

由于B1冰級(jí)極地加強(qiáng)型船舶需要具有航行在碎冰區(qū)的能力，所以船體將遭受冰阻力和水阻力兩大部分。首先忽略冰阻力和水阻力間的耦合作用，假設(shè)冰阻力和水阻力相互獨(dú)立且滿足線性疊加的關(guān)系。此外，考慮到海冰的存在將極大地消除自由液面的興波，并改變船舶尾流的變化等，故在計(jì)算水阻力時(shí)，僅考慮水與船體間的摩擦阻力，忽略比重較小的興波阻力和粘壓阻力的作用。

對(duì)于實(shí)船的摩擦阻力，有較多的計(jì)算方法，較為通用的是1957ITTC公式[20]。根據(jù)ITTC（1995）經(jīng)驗(yàn)公式摩擦水阻力近似為：

式中：為摩擦阻力系數(shù)，為粗糙度補(bǔ)貼系數(shù)；ρ為海水的密度；V為船舶航速；S為船舶的濕表面積，近似有：

式中：Re為雷諾數(shù)；；其中 ν為流體的運(yùn)動(dòng)粘性系數(shù)；LWL為船舶水線間長；B為船寬；T為吃水；CB為方形系數(shù)，有CB=M/(ρLWLBT) ，其中M為船舶排水量。

3 基于船體阻力-速度方法計(jì)算主機(jī)功率

通過計(jì)算船舶阻力與航速的曲線，獲得不同航速下的船體的有效功率為

式中：R為船舶總阻力，分別為冰阻力和摩擦水阻力，分別由式（16）和式（17）求得。考慮到有效功率與主機(jī)功率間的比值即推進(jìn)系數(shù)P.C，有

推進(jìn)系數(shù)P.C是多種效率相乘的綜合名稱，與軸系傳送效率、相對(duì)旋轉(zhuǎn)效率、螺旋槳的敞水效率和船身效率等多種因素有關(guān)，對(duì)于不同的螺旋槳和不同船體會(huì)有一定的變化。因?yàn)檫@里主要考慮冰區(qū)航行船舶的阻力問題，同時(shí)考慮到本船型螺旋槳以及槳機(jī)匹配情況，選取P.C為0.60進(jìn)行計(jì)算，暫不考慮P.C的變化。根據(jù)船體總阻力-速度變化曲線，即可獲得有效功率-速度變化曲線，再根據(jù)式（21）可獲得主機(jī)功率曲線，從而根據(jù)最大航速可獲得船舶所需的主機(jī)功率。

4 計(jì)算結(jié)果及分析

4.1 數(shù)值模擬結(jié)果

根據(jù)上述數(shù)值模型和計(jì)算方法，本小節(jié)根據(jù)中國遠(yuǎn)洋海運(yùn)集團(tuán)有限公司B1冰級(jí)36 000 t冰級(jí)多用途船為對(duì)象，計(jì)算船舶冰阻力以及對(duì)應(yīng)的主機(jī)功率等。在CUDA平臺(tái)上基于GPU并行算法構(gòu)建海冰及船-冰作用離散元模型。船體離散元模型如圖2所示，由2 008塊三角形邊界單元拼接而成。碎冰塊采用球體單元來模擬，主要計(jì)算參數(shù)如表1所示：

表 1 海冰離散元模型參數(shù)Tab. 1 Parameter of discrete elements of ice

用離散元方法模擬的船舶在碎冰航道內(nèi)以不同大小的恒定速度航行，船舶航行速率V分別設(shè)為1 m/s，3 m/s，5 m/s，7 m/s，冰厚H分別設(shè)定為0.8 m，0.5 m和0.3 m。這里選取以速率為5 m/s，冰厚0.8 m為例，船舶在碎冰航道內(nèi)航行過程如圖3所示。

由圖3中可知，在船舶的作用下，碎冰航道將在船舶過后形成具有自由水面的水域，隨著船舶繼續(xù)向前運(yùn)動(dòng)，碎冰在流體的作用下逐漸閉合，再次形成碎冰與水的混合流域?，F(xiàn)再次以冰厚0.8 m為例，獲得不同航速下的船舶冰阻力時(shí)程曲線，如圖4所示，其中橫直線表示航行過程中的平均冰阻力。

圖 3 離散元法模擬的船舶在碎冰航道中的航行過程（V=5 m/s，H=0.8 m）Fig. 3 Numerical simulation of ship moving in the brash ice channel based on DEM （V=5 m/s，H=0.8 m）

從圖4可見，隨著船舶的向前航行，碎冰將在船體上形成高頻、非線性冰載荷作用，這與船舶不同部位與碎冰的不斷接觸、分離以及碎冰沿著船體的滑移等現(xiàn)象有關(guān)。此外，隨著航速的增加，冰阻力均呈現(xiàn)非線性快速增加的趨勢(shì)。

根據(jù)式（17）可計(jì)算不同航速下的摩擦水阻力，將摩擦水阻力和計(jì)算得到的冰阻力疊加可獲得總阻力，采用式（20）可獲得有效功率。本文將不同航速、不同冰厚下阻力值與有效功率計(jì)算值列于表2中。

進(jìn)一步地，可獲得不同阻力值與航速變化以及有效功率與航速變化的曲線分別如圖5和圖6所示。

為了判斷有效功率隨航速變化的特性，這里采用三次多項(xiàng)式對(duì)圖6中有效功率與速度間的關(guān)系進(jìn)行擬合，如下：

借鑒中國商船“永盛輪”首航北極的實(shí)際航行經(jīng)驗(yàn)，對(duì)航行東北航道至北歐地區(qū)的船舶進(jìn)行技術(shù)調(diào)研可知，B1級(jí)冰區(qū)加強(qiáng)型船舶在最大0.8 m左右厚度的碎冰航道中航行的最大速度一般不超過8～10 kn。這里選取最厚冰厚0.8 m，最大航速10 kn計(jì)算，根據(jù)圖6和式（22）計(jì)算可知此時(shí)的有效功率PE約為3 689.73 kW，根據(jù)式（21）計(jì)算可知所需主機(jī)功率NE約為6 149.56 kW。

4.2 與規(guī)范法計(jì)算結(jié)果對(duì)比

圖 4 不同航速時(shí)船舶冰阻力（H=0.8 m）Fig. 4 Ice resistance of ship at different speeds（H=0.8m）

表 2 船舶在不同航速下的阻力及主機(jī)有效功率Tab. 2 Resistance and effective power of main engine of ship at different speeds

圖 5 船舶在不同航速下的摩擦水阻力Rw，冰阻力Ri和總阻力RFig. 5 Friction resistance of water Rw, ice resistance Ri and total resistance R of ship at different speeds

如前所述，冰級(jí)船的主機(jī)功率通常是采用規(guī)范法的經(jīng)驗(yàn)公式進(jìn)行計(jì)算。這里為了對(duì)比和校核，選用芬蘭-瑞典規(guī)范對(duì)上述船型的主機(jī)功率進(jìn)行計(jì)算。根據(jù)芬蘭-瑞典規(guī)范IA級(jí)（等同于B1級(jí)）冰級(jí)船的計(jì)算公式，主機(jī)輸出功率應(yīng)不小于下式所確定的值[6]：

圖 6 船舶在不同航速下的有效功率PE及其擬合曲線Fig. 6 Effective power PE of ship at different speeds and its fitting curves

式中,DP為推進(jìn)器的直徑;ke為推進(jìn)系數(shù)，對(duì)于單個(gè)定螺距螺旋槳選取為2.26；RCH為船舶航行在具有碎冰和凝結(jié)冰層航道中的阻力，由下式確定[6]：

將實(shí)船的參數(shù)代入到式（24）和式（23）中，可計(jì)算得到阻力 RCH約為753.82 kN，最小的主機(jī)輸出功率 NE為6785.87 kW。實(shí)際上，也可以采用CCS規(guī)范[5]計(jì)算上述主機(jī)輸出功率，計(jì)算結(jié)果為6 741.79 kW，與上述結(jié)果略有不同但相差不大。

通過將本文的數(shù)值計(jì)算結(jié)果（6 149.56 kW）與規(guī)范法的計(jì)算結(jié)果（6 785.87 kW）進(jìn)行對(duì)比，可見本文模型計(jì)算獲得的主機(jī)輸出功率比經(jīng)驗(yàn)公式小，說明經(jīng)驗(yàn)公式計(jì)算值相對(duì)保守，且由于經(jīng)驗(yàn)公式中沒有考慮航速的變化，本文模型考慮了航速的影響，從理論上更具合理性。對(duì)于目前這艘B1級(jí)冰區(qū)加強(qiáng)型船舶而言，本文數(shù)值結(jié)果與規(guī)范法相差不大，但是對(duì)于其他船型，本文的計(jì)算結(jié)果明顯小于經(jīng)驗(yàn)公式計(jì)算結(jié)果[7]。

5 結(jié) 語

極地強(qiáng)非線性冰載荷的存在對(duì)極地航行商船的阻力和主機(jī)功率計(jì)算均提出巨大挑戰(zhàn)。本文借鑒中國商船“永盛輪”首航北極的實(shí)際航行經(jīng)驗(yàn)，基于離散元的理論方法和數(shù)值模型，采用阻力-速度曲線法建立了冰區(qū)加強(qiáng)型商船主機(jī)功率計(jì)算方法，獲得主要結(jié)論如下：

1）通過規(guī)范法可以粗估冰級(jí)航行船舶在碎冰航道內(nèi)遭受的阻力以及所需的主機(jī)功率，規(guī)范法的公式簡單，計(jì)算較為簡便，但不同規(guī)范的計(jì)算略有不同，計(jì)算的阻力值以及主機(jī)功率與航速無關(guān)，計(jì)算結(jié)果相對(duì)保守。

2）基于離散元方法，建立了冰區(qū)加強(qiáng)型船舶在碎冰航道內(nèi)直航的運(yùn)動(dòng)模型，能夠預(yù)報(bào)不同冰情條件下的冰區(qū)加強(qiáng)型船舶的冰阻力，并與船舶摩擦水阻力疊加，獲得冰區(qū)加強(qiáng)型船舶的總阻力，從而可以獲得總阻力-航速曲線以及有效功率和主機(jī)功率-速度曲線。該方法考慮了航速和冰情對(duì)于主機(jī)功率的影響，理論上更具合理性。

3）根據(jù)“永盛輪”首航北極航道的冰情數(shù)據(jù)和航行經(jīng)驗(yàn)，將冰情與航速等數(shù)據(jù)輸入到本文數(shù)值模型中，可以應(yīng)用該數(shù)值模型獲得總阻力，從而獲得相應(yīng)的有效功率和主機(jī)功率，為冰區(qū)加強(qiáng)型船舶主機(jī)功率的選取提供創(chuàng)新性方法。計(jì)算結(jié)果表明，該方法計(jì)算獲得的主機(jī)功率較規(guī)范法小；

4）應(yīng)用此方法，對(duì)于已經(jīng)選定了主機(jī)的冰區(qū)船舶，從主機(jī)實(shí)際能夠持續(xù)輸出最大功率的能力角度，還可以給出不同冰情下船舶最大航速的建議值，可為我國商船極地水域航行提供技術(shù)支撐。