狄少丞，季順迎，薛彥卓

(1. 哈爾濱工程大學(xué) 船舶工程學(xué)院，黑龍江 哈爾濱 150001；2. 大連理工大學(xué) 工業(yè)裝備結(jié)構(gòu)分析國(guó)家重點(diǎn)實(shí)驗(yàn)室，遼寧 大連 116024)

船舶在平整冰區(qū)行進(jìn)過(guò)程的離散元分析

狄少丞1，季順迎2，薛彥卓1

(1. 哈爾濱工程大學(xué) 船舶工程學(xué)院，黑龍江 哈爾濱 150001；2. 大連理工大學(xué) 工業(yè)裝備結(jié)構(gòu)分析國(guó)家重點(diǎn)實(shí)驗(yàn)室，遼寧 大連 116024)

船舶在冰區(qū)航行中，冰荷載是船舶結(jié)構(gòu)設(shè)計(jì)和安全航行的重要影響因素。采用具有黏結(jié)破碎特征的離散單元法對(duì)船舶在平整冰區(qū)的航行過(guò)程進(jìn)行了數(shù)值分析。在船舶破冰過(guò)程的離散元模擬中，平整冰由球形顆粒黏結(jié)而成，同時(shí)考慮了海水浮力和拖曳力對(duì)海冰單元的作用，顆粒間的黏結(jié)作用可在海冰與船體的相互作用過(guò)程中發(fā)生破碎；船體由三角形單元構(gòu)造，分析了在海水浮力、推進(jìn)器推力及海冰共同作用下船體的六自由度運(yùn)動(dòng)特性。通過(guò)船體與海冰相互作用的離散元分析，確定了船舶在破冰過(guò)程中海冰作用在船體上的冰荷載，探討了船舶推進(jìn)器推力和海冰厚度對(duì)冰荷載的影響，并與Lindqvist船體冰阻力經(jīng)驗(yàn)公式進(jìn)行了對(duì)比。

平整冰；船體；離散元方法；冰荷載；破冰過(guò)程；海冰厚度

Abstract: During the navigation of ship in ice covered regions, the ice load on ship hull is a significant factor for the structure design and safe maneuvering of ships. In this study, the discrete element method (DEM) is applied to simulate the interaction between level ice and ship hull. The level ice is constructed with bonded spherical elements under the buoyancy and drag force of current. The bonded elements can be broken in the process of the interaction between level ice and ship hull. The ship hull is constructed with triangular elements and treated as a six-degree of freedom system. With the DEM simulations of the interaction between level ice and ship hull, the global ice resistance on ship hull in the process of ice breaking can be determined. The influences of propulsive force and ice thickness on ice load are also discussed based on the simulation results. The comparison of ice load between DEM simulation results and the empirical formula proposed by Lindqvist is presented. This study can be an aid to the structure design of ship and the navigation security in ice-covered fields.

Keywords: level ice; ship hull; discrete element method; ice load; ice breaking process; ice thickness

近年來(lái)，受全球氣候變暖的影響，北極海冰覆蓋面積持續(xù)減少，北極航道的夏季通航也成為可能[1-2]，北極東北航道和西北航道的開(kāi)通對(duì)我國(guó)海上運(yùn)輸有著重要的戰(zhàn)略意義[3-5]。為此，我國(guó)從2012年開(kāi)始在北冰洋的考察中針對(duì)北極航道開(kāi)通的可行性進(jìn)行了系統(tǒng)的調(diào)查[4]。同時(shí)，北極海域蘊(yùn)藏著豐富的油氣資源，其油氣儲(chǔ)量相當(dāng)于目前世界探明原油儲(chǔ)量的13%和天然氣儲(chǔ)量的30%[6]。因此，冰區(qū)航行的安全保障和抗冰船舶的設(shè)計(jì)是當(dāng)前亟待解決的工程問(wèn)題。目前我國(guó)在冰區(qū)船舶航行和船體冰力研究方面還處于相對(duì)薄弱階段。

船舶在冰區(qū)行進(jìn)過(guò)程中，受海冰類型、船型和航行方式等因素的影響，使得海冰與船體的作用過(guò)程極為復(fù)雜。與船體作用時(shí)海冰通常會(huì)經(jīng)歷斷裂、翻轉(zhuǎn)、滑動(dòng)清除三個(gè)階段[7-9]，同時(shí)會(huì)伴隨有海冰在船艏的堆積、重疊等過(guò)程[10-11]。此外，海冰對(duì)船體冰力的大小與船舶的錨泊、直行、轉(zhuǎn)向等航行方式密切相關(guān)[12-14]。在船體冰荷載研究中，通常采用現(xiàn)場(chǎng)冰荷載測(cè)量、室內(nèi)模型試驗(yàn)、理論分析與數(shù)值分析的研究方法來(lái)揭示海冰與船體的作用機(jī)理[15-17]。而數(shù)值分析方法從揭示海冰破壞機(jī)理到船體總冰力變化規(guī)律，是研究海冰與船體相互作用的重要途徑。

在船體冰荷載的數(shù)值分析中，人們發(fā)展了針對(duì)不同海冰形態(tài)特征和力學(xué)特性的數(shù)值模型[18-20]。由于海冰在與船體的相互作用過(guò)程中會(huì)呈現(xiàn)出由連續(xù)體向離散塊體轉(zhuǎn)變的破壞過(guò)程[21]，且在不同尺度下海冰又呈現(xiàn)出離散分布特性，而在碎冰區(qū)、重疊和堆積冰區(qū)離散特性更為明顯[22-24]，這使得離散元方法在確定船體冰荷載方面具有明顯的優(yōu)勢(shì)[25-27]。

為此，下文將針對(duì)海冰與船體相互作用的動(dòng)力特性，建立具有黏結(jié)—破碎特征的海冰離散元方法，對(duì)海冰與船體作用時(shí)的破壞模式、船體總冰力進(jìn)行數(shù)值模擬，進(jìn)而分析船舶在平整冰區(qū)行進(jìn)時(shí)的船舶航速及船舶破冰能力。

1 海冰與船體作用的離散元模型

1.1平整冰的離散元模型

圖1 顆粒間的平行黏結(jié)模型Fig. 1 Parallel bonding model between two spherical particles

1.2船體的離散元模型

在用離散單元法研究海冰與船體的相互作用過(guò)程中，先將船體離散成三角形單元的集合，離散后的各三角形單元通過(guò)頂點(diǎn)相連接，通過(guò)設(shè)定海冰顆粒與三角形單元之間的接觸準(zhǔn)則即可得到海冰與船體之間的作用力。

圖2 三角形單元Fig. 2 Triangular element

圖3 球心在三角形面內(nèi)的投影Fig. 3 Projection of center of sphere on triangle

1.2.1 接觸判斷

若要判定三角形單元f與某一球形顆粒p是否發(fā)生接觸，首先需判定球形顆粒是否與三角形單元所在的平面有重疊。定義球形顆粒p的球心P坐標(biāo)為Xp，顆粒半徑為Rball，球心P在三角形單元所在平面的投影點(diǎn)為Q，如圖3所示。

可以得到以下向量：AP、BP、CP、BA、AB、CA、AC、CB、BC、PQ，其中PQ=AP·n，如果PQ

如果球形顆粒與三角形單元所在的平面有重疊，需進(jìn)一步判斷球形顆粒是否與三角形單元本身接觸，并判斷它們的接觸類型，以正確地計(jì)算顆粒與三角形單元之間的接觸力。球形顆粒與三角形單元發(fā)生接觸的類型有三種：球-面接觸、球-邊接觸、球-頂點(diǎn)接觸，以下簡(jiǎn)稱為面接觸、邊接觸和角接觸。面接觸是指球形顆粒球心的投影點(diǎn)在三角形面內(nèi)；邊接觸是指球形顆粒與三角形單元的三條邊接觸；角接觸是指球形顆粒與三角形單元的三個(gè)頂點(diǎn)接觸，三種接觸類型如圖4所示。

圖4 球形顆粒與三角形單元的接觸類型Fig. 4 Contact types between spherical particle and triangular element

圖5 判定投影點(diǎn)是否在三角形單元內(nèi)部Fig. 5 Determination of the location of the projection point inside the triangle

采用重心法來(lái)確定球心的投影點(diǎn)是否在三角形單元內(nèi)部，如圖5所示，稱之為第二次判定。對(duì)于平面內(nèi)任意一點(diǎn)，都可以表示為：

如果投影點(diǎn)Q位于三角形ABC內(nèi)部，u和v必須同時(shí)滿足以下三個(gè)條件：

如果系數(shù)u和v滿足式(4)，則表明投影點(diǎn)Q在三角形單元內(nèi)部，進(jìn)而說(shuō)明球形顆粒與三角形單元發(fā)生面接觸。

如果投影點(diǎn)Q位于三角形邊內(nèi)部并且球心到投影點(diǎn)的距離小于顆粒半徑Rball時(shí)，表明二者發(fā)生邊接觸，如圖6所示，稱之為第三次判定。

以AB邊為例，點(diǎn)P在AB邊上的投影點(diǎn)為Q，則

因此，當(dāng)AP·BABP·BA≤0且QP≤Rball，顆粒與三角形單元發(fā)生邊接觸；當(dāng)AP·BABP·BA>0或QP>Rball，顆粒與三角形單元不發(fā)生邊接觸，每個(gè)三角形單元有三條邊，需做三次判定。

如果球形顆粒與三角形單元的接觸關(guān)系滿足第一次判定而不滿足第二、第三次判定，則表明顆?？赡芘c三角形單元的頂點(diǎn)發(fā)生接觸。如果顆粒球心P與頂點(diǎn)A、B、C之間的距離小于顆粒的半徑Rball，則發(fā)生角接觸，如圖7所示。

圖6 球形顆粒與三角形邊接觸Fig. 6 Contact between spherical particle and triangle edge

圖7 球形顆粒與三角形頂點(diǎn)接觸Fig. 7 Contact between spherical particle and triangle vertex

球形顆粒與三角形單元發(fā)生角接觸的判斷準(zhǔn)則為：AP≤Rball，與頂點(diǎn)A發(fā)生角接觸；BP≤Rball，與頂點(diǎn)B發(fā)生角接觸；CP≤Rball與頂點(diǎn)C發(fā)生角接觸。

1.2.2 接觸力計(jì)算

在確定球形顆粒與三角形單元的接觸類型后，即可計(jì)算出球體單元在三角形單元中的嵌入量ΔL：

式中：點(diǎn)PC代表球形顆粒與邊界單元的接觸點(diǎn)，分別為點(diǎn)Q(面接觸、棱邊接觸)；A、B、C(頂點(diǎn)接觸)。顆粒與三角形單元接觸點(diǎn)之間的法向量可表示為:

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

其中，vp為顆粒速度，vw為三角形單元與顆粒接觸點(diǎn)的速度。將相對(duì)位移沿接觸面分解為法向分量Δxn和切向分量Δxs：

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

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

1.2.3 船體浮力及拖曳力計(jì)算

在本文工作中，考慮海水浮力和海流拖曳力對(duì)船體運(yùn)動(dòng)的影響，浮力的計(jì)算采用阿基米德原理，計(jì)算公式如下

式中：ρw為海水密度，g為重力加速度，Vsub為船體浸入水中的體積。

在離散元模型中，船體由若干三角形單元拼接而成，在計(jì)算船體浸入水中的體積時(shí)采用如下方法：取水面上任一點(diǎn)Pw為頂點(diǎn)，以水線以下三角形為底面做四面體，所有四面體體積之和即為船體浸入水中的體積，四面體的體積(如圖8所示)可由下式計(jì)算：

圖8 四面體單元的體積Fig. 8 Volume of tetrahedron

式中：a=v2-v1、b=v3-v1、r=p-v1，并且四面體的體心為：

因此，浸入水中船體的體積可表示為：

浸入水中部分的體心可表示為：

三角單元離散化后的船體中會(huì)有部分露出水面的三角形板，在體積計(jì)算時(shí)要特殊考慮。對(duì)于圖9中所示的未完全浸入水中的三角形板，需要將露出水面的部分切除掉，只計(jì)算水下部分的體積。

圖9 切除四面體露出水面的部分Fig. 9 Clipping the tetrahedron along the water line

至此求得了船體浸入水中的體積Vsub，進(jìn)而由式(14)可求得船體在水中的浮力。船體在水中所受的拖曳作用由兩部分組成：船體發(fā)生平動(dòng)時(shí)的拖曳力Fd，以及旋轉(zhuǎn)時(shí)的拖曳力矩Md。二者可表示為：

式中：Cα、Cβ分別為對(duì)應(yīng)于拖曳力與拖曳力矩的拖曳系數(shù)，vw為水流速度，vs為船體平移的速度，ωs為船體轉(zhuǎn)動(dòng)的角速度，L為船體型長(zhǎng)。

對(duì)比分析2014年9月—2015年8月中國(guó)大陸IMERG和CGDPA日平均降水量,結(jié)果表明:兩者的空間分布比較一致,東部介于1～10 mm,西部介于0.1～5 mm,南部介于2～10 mm,北部介于0.1～5 mm,說(shuō)明IMERG日平均降水量能夠較好地反映中國(guó)大陸降水量西部少、東部多,南部多、北部少的特點(diǎn);IMERG日平均降水量比CGDPA的空間分布連續(xù)性更好,尤其是在西部、北部等地面雨量計(jì)觀測(cè)站點(diǎn)較少的地區(qū)(圖1)。西藏東南部、新疆天山,以及青海、廣西和廣東南部和海南等地區(qū),IMERG日平均降水量比CGDPA偏低。

2 船舶在平整冰區(qū)直航的離散元分析

2.1船體冰阻力的離散元模擬

以“雪龍”號(hào)極地考察船為研究對(duì)象，如圖10所示，建立船體的離散元模型。由于在航行過(guò)程中只有船體外殼與海冰發(fā)生接觸，因此只將船體外殼三角離散化作為“雪龍”船的離散元模型，建立的船體離散元模型如圖11所示，船體離散元模型由1 430塊三角形單元拼接而成。

圖10 “雪龍”極地考察船F(xiàn)ig. 10 “Snow Dragon” polar expedition vessel

圖11 船體的三角形單元離散化Fig. 11 Construction of ship hull with triangular elements

在“雪龍”船與平整冰相互作用的離散元模擬中，采用兩層規(guī)則排列的球形顆粒組合來(lái)構(gòu)造平整冰，平整冰排列方式如圖12(a)所示。海冰的水動(dòng)力特性只考慮了浮力和拖曳力，海水對(duì)海冰的浮力和拖曳力分別加到了每個(gè)顆粒上，如果顆粒全部浸入水中，浮力為整個(gè)顆粒受到的浮力；如果顆粒部分浸入水中，浮力為水面以下部分受到的浮力，拖曳力的計(jì)算與浮力類似。在模擬過(guò)程中平整冰保持固定不動(dòng)，船體在推進(jìn)器作用下向冰排移動(dòng)。船體在運(yùn)動(dòng)過(guò)程中會(huì)受到海水浮力和海流拖曳力、推進(jìn)器作用力及海冰的阻力，在這些力的共同作用下，船體會(huì)發(fā)生沿三個(gè)互相垂直方向上的平移和轉(zhuǎn)動(dòng)，規(guī)定船體行進(jìn)方向?yàn)閄方向，垂直于水面的方向?yàn)閆方向，X、Y、Z三個(gè)方向符合右手定則。模擬中采用的主要計(jì)算參數(shù)列于表1中。

表1 海冰與船體結(jié)構(gòu)作用離散元模擬的主要計(jì)算參數(shù)Tab. 1 Major computational parameters of DEM simulation of dynamic ice load on ship hull

離散元模擬的不同時(shí)刻破冰船在平整冰區(qū)的航行狀態(tài)如圖12所示，船舶在恒定推力作用下駛?cè)肫秸鶇^(qū)時(shí)，海冰在船體的沖擊下發(fā)生破碎，并在冰區(qū)內(nèi)形成一條狹長(zhǎng)的水道。船舶在駛?cè)牒婉偝霰鶇^(qū)過(guò)程中，海冰對(duì)船體在X、Y、Z三個(gè)方向上的冰阻力如圖13所示?？梢钥闯?，當(dāng)船體剛與平整冰接觸時(shí)，船體在X、Z方向上的冰阻力隨航行距離的增加而逐漸增大；隨著船體的駛?cè)耄枇χ饾u達(dá)到最大值并穩(wěn)定在一定范圍內(nèi)，呈現(xiàn)出很強(qiáng)的脈動(dòng)特性；隨著船體駛離冰區(qū)，冰阻力又開(kāi)始下降，最終降為0。

由于船體在Y方向上具有對(duì)稱性，使得在Y方向上冰阻力亦具有對(duì)稱特性。

圖12 離散元模擬的“雪龍”船在平整冰中的破冰行進(jìn)過(guò)程Fig. 12 The navigation of “Snow Dragon” ship in level ice-covered regions simulated with DEM method

將圖12中t= 300 s時(shí)刻的船舶破冰過(guò)程進(jìn)行放大顯示，放大后的破冰過(guò)程如圖14所示?？梢钥闯?，船舶在行進(jìn)過(guò)程中與船體接觸處的海冰發(fā)生擠壓破壞，而遠(yuǎn)離船體的海冰發(fā)生彎曲破壞，圖14(a)中的白色折線表示海冰發(fā)生彎曲破壞所處的位置。破碎后的海冰主要集中在船體兩側(cè)(圖14(b))，有部分海冰滑入船底(圖14(a))，船體兩側(cè)和滑入船底的海冰均會(huì)對(duì)船體產(chǎn)生摩擦作用。

在推進(jìn)器推力作用下船體從靜止開(kāi)始加速，船舶航速隨時(shí)間變化曲線如圖15所示。船舶從開(kāi)闊水域開(kāi)始加速運(yùn)動(dòng)，由于此時(shí)船體并未與海冰接觸，因此船速逐漸增大(0 ～ 22 s)；隨著船體與海冰發(fā)生接觸，船體冰阻力逐漸增大，此時(shí)船舶航速呈現(xiàn)先增大后減小的運(yùn)動(dòng)趨勢(shì)(22 ～ 300 s)；當(dāng)推進(jìn)器推力、海水阻力及冰阻力三者達(dá)到平衡時(shí)，船速開(kāi)始保持穩(wěn)定(300 ～ 500 s)，本文將此時(shí)間段內(nèi)的船體平均速度作為在給定推進(jìn)器推力和冰厚條件下船體在平整冰區(qū)航行的速率，對(duì)應(yīng)此工況下的船舶航速為1.03 m/s；最后船舶駛離冰區(qū)進(jìn)入開(kāi)闊水域，船速又進(jìn)一步增大(500 ～ 600 s)。從圖13(a)的冰阻力時(shí)程曲線可以看出，在船舶破冰過(guò)程中，冰阻力值在600 ～ 1500 kN范圍內(nèi)波動(dòng)，但在主要破冰時(shí)間段(50～500 s)內(nèi)，船舶冰阻力的平均值約為930 kN，與推力940 kN保持一致，可保證在后續(xù)時(shí)間段內(nèi)船舶破冰前行。

圖13 冰厚1.2 m，推力940 kN時(shí)船舶冰阻力時(shí)程曲線Fig. 13 Time history of ice resistance on ship hull with thickness of 1.2 m, propulsive force 940 kN

圖14 放大顯示的“雪龍”船破冰過(guò)程Fig. 14 The detail of ice breaking process simulated with DEM method

圖15 冰厚1.2 m，推力940 kN時(shí)船舶行進(jìn)速度時(shí)程曲線Fig. 15 Time history of navigational speed of ship with thickness of 1.2 m, propulsive force 940 kN

2.2船速、冰厚對(duì)冰阻力的影響

圖17 冰厚-航速-推進(jìn)力關(guān)系Fig. 17 Relationship between ice thickness, navigational speed and propulsive force

船舶在冰區(qū)航行時(shí)，航速和冰厚是影響冰阻力的重要因素。在同一推進(jìn)器推力作用下，在不同冰厚的平整冰區(qū)航行時(shí)船舶的航速也將不同。本節(jié)將分析在不同推力作用下船舶在不同冰厚的平整冰區(qū)航行時(shí)的航速，平整冰的厚度分別設(shè)為0.7、0.8、0.9、1.0、1.1和1.2 m，計(jì)算中顆粒之間的黏結(jié)強(qiáng)度設(shè)定為同一強(qiáng)度，海冰在厚度方向上的顆粒層數(shù)均設(shè)為兩層，采用相同的顆粒層數(shù)可保證海冰在破碎過(guò)程中具有相同的破碎特性[30]。推進(jìn)器推力的設(shè)定取決于平整冰的厚度，平整冰厚度較大時(shí)設(shè)定較大的推力值，厚度較小時(shí)設(shè)定較小的推力值。采用離散元法模擬的船舶在平整冰區(qū)航行時(shí)的冰厚-航速-推力關(guān)系如圖17所示，圖中橫坐標(biāo)Vship為船舶在不同厚度的平整冰中航行時(shí)最終達(dá)到的穩(wěn)定速度，縱坐標(biāo)Fp為對(duì)應(yīng)的推進(jìn)器推力。

從圖17中可以看出，“雪龍”船在較薄的平整冰區(qū)航行時(shí)只需較小的推進(jìn)力即可獲得較大的破冰速度；當(dāng)平整冰較厚時(shí)，只有增大推進(jìn)器的推力才能獲得更大的破冰速度。根據(jù)“雪龍”船的航行經(jīng)驗(yàn)，可知“雪龍”船可在1.2 m厚的冰層(含0.2 m厚的雪)中連續(xù)破冰前行，由圖17可知，當(dāng)推進(jìn)器推力為940 kN時(shí)，可獲得破冰速率為1.03 m/s(約2節(jié))；當(dāng)推進(jìn)器推力為1 100 kN時(shí)，可獲得破冰速率為2.44 m/s(約4.8節(jié))。

2.3與經(jīng)驗(yàn)冰阻力公式的比較

Lindqvist冰力公式[31]是目前使用較多且較為準(zhǔn)確的破冰船阻力計(jì)算經(jīng)驗(yàn)公式，公式中考慮了首柱處海冰擠壓破壞產(chǎn)生的阻力，涉及到的參數(shù)包括主尺度、船型、冰厚、摩擦以及海冰的彎曲強(qiáng)度。

破冰阻力Rc：

彎曲阻力Rb：

式中：hi為海冰厚度；E為楊氏模量；υ為泊松比；σf為海冰彎曲強(qiáng)度；B為船體型寬；ρw為海水密度。

浸沒(méi)阻力Rs：

式中：ρi為海冰密度；htot為冰與雪的厚度和；T為船體吃水深度；L為船體長(zhǎng)度。

總冰阻力Ri：

圖18 冰厚-航速-冰阻力關(guān)系Fig. 18 Relationship between ice thickness, navigational speed and ice resistance

離散元模擬的“雪龍”船在平整冰區(qū)航行時(shí)船體冰阻力情況如圖18所示，圖19為不同海冰厚度下船舶推進(jìn)力與破冰阻力之間的對(duì)應(yīng)關(guān)系，其中Vship為船體在不同厚度的平整冰中航行時(shí)的速度，縱坐標(biāo)Rice為船體受到的冰阻力。從圖中可以看出，船體冰阻力隨航速的增加呈增大趨勢(shì)，但并非線性增加。在本文船體冰阻力的離散元模型中只考慮了海冰與船體的浮力與拖曳力，未考慮二者的水動(dòng)力特性，這就導(dǎo)致了冰阻力隨航速的增大趨勢(shì)不明顯。因此本文在對(duì)比分析中不考慮船體航速的影響，而只將不同航速下的平均冰阻力與破冰阻力、彎曲阻力、浸沒(méi)阻力三者之和作比較，比較結(jié)果如圖20所示。

圖19 不同海冰厚度下船舶推進(jìn)力與冰阻力的關(guān)系Fig. 19 Relationship between propulsive force and ice resistance under different ice thicknesses

圖20 Lindqvist經(jīng)驗(yàn)計(jì)算冰阻力與DEM模擬冰阻力的對(duì)比Fig. 20 Comparison of the ice resistance computed with Lindqvist formula and simulated with DEM method

對(duì)比分析中用到的計(jì)算參數(shù)列于表2中，表中彎曲強(qiáng)度值通過(guò)海冰的懸臂梁數(shù)值模型試驗(yàn)得到[32]，離散元模型中顆粒間的黏結(jié)強(qiáng)度設(shè)為0.5 MPa(表1)，經(jīng)數(shù)值試驗(yàn)得到的宏觀海冰彎曲強(qiáng)度為0.93 MPa(表2)。

表2 Lindqvist冰力經(jīng)驗(yàn)公式中的主要計(jì)算參數(shù)Tab. 2 Major computational parameters used in Lindqvist empirical formula

從圖20可以看出，在不考慮船體水動(dòng)力特性的條件下，船體冰阻力隨所破海冰厚度的增加而增大；離散元計(jì)算的船體冰阻力與Lindqvist經(jīng)驗(yàn)冰力公式相比，二者較吻合，表明該離散元方法可合理模擬船舶在平整冰區(qū)的破冰過(guò)程，進(jìn)而確定船舶的冰阻力。

3 結(jié) 語(yǔ)

采用離散單元法對(duì)船舶在平整冰區(qū)的行進(jìn)過(guò)程進(jìn)行了數(shù)值模擬，分析了船體在不同推進(jìn)力、不同厚度的平整冰區(qū)內(nèi)行進(jìn)速度的變化規(guī)律，以及船體在行進(jìn)過(guò)程中的冰阻力。以上研究表明，上文建立的離散元法可合理計(jì)算船舶在冰區(qū)的行進(jìn)過(guò)程，可對(duì)海冰與船體作用過(guò)程中的破碎現(xiàn)象進(jìn)行模擬，進(jìn)而確定冰區(qū)船舶的破冰能力。在海冰破碎過(guò)程的離散元模擬中，模型細(xì)觀參數(shù)的選取對(duì)模擬結(jié)果有很大的影響，例如，單元尺寸、接觸模型、破壞準(zhǔn)則等因素都會(huì)對(duì)船體的冰荷載產(chǎn)生影響。因此在后續(xù)工作中需進(jìn)一步考慮這些因素的影響以完善海冰的離散元模型，在此基礎(chǔ)上進(jìn)一步考慮船體及碎冰的水動(dòng)力特性，以更精確地模擬船舶在冰區(qū)的破冰過(guò)程。以上工作有助于揭示船體冰阻力的變化規(guī)律，為后續(xù)研究不同航行狀態(tài)、不同冰況下船體冰阻力提供了有效的數(shù)值計(jì)算途徑。

[1] CHIRCOP A. Climate change and the prospects of increased navigation in the Canadian Arctic[J]. Journal of Maritime Affairs, 2007, 6(2): 193-205.

[2] HO J. The implications of Arctic sea ice decline on shipping[J]. Marine Policy, 2010, 34:713-715.

[3] 蘇潔, 徐棟, 趙進(jìn)平, 等. 北極加速變暖條件下西北航道的海冰分布變化特征[J]. 極地研究, 2010, 22(2):104-124. (SU Jie, XU Dong, ZHAO Jinping, et al. Features of northwest passage sea ice’s distribution and variation under arctic rapidly warming condition[J]. Chinese Journal of Polar Research, 2010, 22(2):104-124. (in Chinese))

[4] 王洛, 趙越, 劉建民, 等. 中國(guó)船舶首航東北航道及展望[J]. 極地研究, 2014, 26(2): 276-284. (WANG Luo, ZHAO Yue, LIU Jianmin, et al. China’s first trans-arctic voyage and related expectations[J]. Chinese Journal of Polar Research, 2014, 26(2): 276-284. (in Chinese))

[5] 何劍鋒, 吳榮榮, 張芳, 等. 北極航道相關(guān)海域科學(xué)考察研究進(jìn)展[J]. 極地研究, 2012, 24(2): 187-196. (HE Jianfeng, WU Rongrong, ZHANG Fang, et al. The progress of expeditions and research in the seas related to the arctic passages[J]. Chinese Journal of Polar Research, 2012, 24(2): 187-196. (in Chinese))

[6] CHIERMEIER Q. The great Arctic oil race begins[J]. Nature, 2012, 482:13-14.

[7] SU B, SKJETNE R, BERG T E. Numerical assessment of a double-acting offshore vessel's performance in level ice with experimental comparison[J]. Cold Regions Science and Technology, 2014, 106/107:96-109.

[8] KUJALA P, ARUGHADHOSS S. Statistical analysis of ice crushing pressures on a ship's hull during hull-ice interaction[J]. Cold Regions Science and Technology, 2012, 70: 1-11.

[9] SUYUTHI A, LEIRA B J, RISKA K. A generalized probabilistic model of ice load peaks on ship hulls in broken-ice fields[J]. Cold Regions Science and Technology, 2014, 97:7-20.

[10] HANSEN E H, LOSET S. Modelling floating offshore units moored in broken ice: comparing simulations with ice tank tests[J]. Cold Regions Science and Technology, 1999, 29: 107-119.

[11] HANSEN E H, L?SET S. Modelling floating offshore units moored in broken ice: model description[J]. Cold Regions Science and Technology, 1999, 29: 97-106.

[12] IBRAHIM R A, CHALHOUB N G, FALZARANO J. Interaction of ships and ocean structures with ice loads and stochastic ocean waves[J]. Applied Mechanics Review, 2007, 60(5): 246-290.

[13] KARULIN E B, KARULINA M M. Numerical and physical simulations of moored tanker behaviour[J]. Ship and Offshore Structures, 2011, 6(3): 179-184.

[14] ZHOU L, RISKA K, POLACH R V B U, et al. Experiments on level ice loading on an icebreaking tanker with different ice drift angles[J]. Cold Regions Science and Technology, 2013, 85(1): 79-93.

[15] GAGNON R, CUMMING D, RITCH R, et al. Overview accompaniment for papers on the barge bit impact trials [J]. Cold Regions Science and Technology, 2008, 52(1): 1-6.

[16] RALPH F, MCKENNA R, GAGNON R. Iceberg characterization for the barge bit impact study[J]. Cold Regions Science and Technology, 2008, 52(1): 7-28.

[17] KJERSTAD O K, METRIKIN I, L?SET S, et al. Experimental and phenomenological investigation of dynamic positioning in managed ice[J]. Cold Regions Science and Technology, 2015, 111: 67-79.

[18] GAGNON R E, WANG J. Numerical simulations of a tanker collision with a bergy bit incorporating hydrodynamics, a validated ice model and damage to the vessel[J]. Cold Regions Science and Technology, 2012, 81: 26-35.

[19] AKSNES V. A simplified interaction model for moored ships in level ice[J]. Cold Regions Science and Technology, 2010, 63: 29-39.

[20] SU B, RISKA K, MOAN T. Numerical simulation of local ice loads in uniform and randomly varying ice conditions[J]. Cold Regions Science and Technology, 2011, 65:145-159.

[21] 王鈺涵, 李輝, 任慧龍, 等. 連續(xù)破冰模式下破冰船的冰力研究[J]. 海洋工程, 2013, 31(4): 68-73. (WANG Yuhan, LI Hui, REN Huilong, et al. Study of ice force about icebreaker based on continuous breaking pattern[J]. The Ocean Engineering, 2013, 31(4): 68-73. (in Chinese))

[22] DEMPSEY J P. Research trends in ice mechanics[J]. International Journal of Solids and Structures, 2000, 37: 131-153.

[23] WEISS J, SCHULSON E, STERN H. Sea ice rheology from in-situ, satellite and laboratory observations: fracture and friction[J]. Earth and Planetary Science Letters, 2007, 255:1-8.

[24] HIBLER W D. Sea ice fracturing on the large scale[J]. Engineering Fracture Mechanics, 2001, 68: 2013-2043

[25] LAU M, LAWRENCE K P, ROTHENBURG L. Discrete element analysis of ice loads on ships and structures[J]. Ships and Offshore Structures, 2011, 6(3): 211-221.

[26] HOPKINS M A, SHEN H H. Simulation of pancake-ice dynamics in a wave field[J]. Annals of Glaciology, 2001, 33: 355-360.

[27] SUN S, SHEN H H. Simulation of pancake ice load on a circular cylinder in a wave and current field[J]. Cold Regions Science and Technology, 2012, 78: 31-39.

[28] POTYONDY D O, CUNDALL P A. A bonded-particle model for rock[J]. International Journal of Rock Mechanics & Mining Sciences, 2004, 41:1329-1364.

[29] JI S Y, DI S C, LONG X. DEM simulation of uniaxial compressive and flexural strength of sea ice: parametric study[J]. Journal of Engineering Mechanics, 2016: 10.1061/(ASCE)EM.1943-7889.0000996.

[30] FENG Y T, OWEN D R J. Discrete element modeling of large scale particle systems-I: exact scaling laws[J]. Computational Particle Mechanics, 2014, 1(2): 159-168.

[31] LINDQVIST G. A straightforward method for calculation of ice resistance of ships[C]//Proceedings of POAC, 1989:722-735.

[32] SCHWARZ J, FREDERKING R, GAVRILLO V, et al. Standardized testing methods for measuring mechanical properties of ice[J]. Cold Regions Science and Technology, 1981, 4(3):245-253.

Analysis of ship navigation in level ice-covered regions with discrete element method

DI Shaocheng1, JI Shunying2, XUE Yanzhuo1

(1. College of Shipbuilding Engineering, Harbin Engineering University, Harbin 150001, China; 2. State Key Laboratory of Structural Analysis for Industrial Equipment, Dalian University of Technology, Dalian 116024, China)

U674.21; U661.31

A

10.16483/j.issn.1005-9865.2017.03.08

1005-9865(2017)03-0059-11

2016-04-12

國(guó)家自然科學(xué)基金資助項(xiàng)目(51639004；51579054；41606213)

狄少丞(1986-)，男，山西五寨人，講師，從事工程海冰的數(shù)值模型研究。E-mail：dishaocheng@hrbeu.edu.cn

薛彥卓。E-mail：xueyanzhuo@hrbeu.edu.cn