冰級船主機功率理論預報方法研究

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

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

0 引 言

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

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

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

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

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

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

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

1.2 碎冰離散元模型

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

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

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

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

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

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

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

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

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

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

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

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

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

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

式中，P代表顆粒點坐標，點PC代 表P對應的接觸點坐標[19]，Rball為海冰球狀顆粒的半徑。顆粒與三角形單元接觸點之間的法向量可表示為：

相對位移Δx可表示為：

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

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

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

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

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

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

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

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

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

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

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

3 基于船體阻力-速度方法計算主機功率

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

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

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

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

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

根據(jù)上述數(shù)值模型和計算方法，本小節(jié)根據(jù)中國遠洋海運集團有限公司B1冰級36 000 t冰級多用途船為對象，計算船舶冰阻力以及對應的主機功率等。在CUDA平臺上基于GPU并行算法構(gòu)建海冰及船-冰作用離散元模型。船體離散元模型如圖2所示，由2 008塊三角形邊界單元拼接而成。碎冰塊采用球體單元來模擬，主要計算參數(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ù)向前運動，碎冰在流體的作用下逐漸閉合，再次形成碎冰與水的混合流域。現(xiàn)再次以冰厚0.8 m為例，獲得不同航速下的船舶冰阻力時程曲線，如圖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)非線性快速增加的趨勢。

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

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

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

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

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

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

表 2 船舶在不同航速下的阻力及主機有效功率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

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

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

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

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

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

5 結(jié) 語

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

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

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

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

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