任奕舟， 鄒早建,2

(1. 上海交通大學(xué) 船舶海洋與建筑工程學(xué)院,上?！?00240； 2.上海交通大學(xué) 海洋工程國家重點實驗室,上?！?00240)

?

破冰船在冰層中連續(xù)破冰過程的數(shù)值模擬

任奕舟1， 鄒早建1,2

(1. 上海交通大學(xué) 船舶海洋與建筑工程學(xué)院,上海200240； 2.上海交通大學(xué) 海洋工程國家重點實驗室,上海200240)

介紹了一種用于模擬破冰船在冰層中連續(xù)破冰的冰材料數(shù)值模型。通過與冰錐受壓實驗數(shù)據(jù)進行對比，對該數(shù)值模型進行了驗證。將該模型應(yīng)用于破冰船在無限冰區(qū)與冰層碰撞的數(shù)值模擬，對破冰船的破冰阻力進行了計算，并將不同船首、不同冰層厚度下計算所得的破冰阻力與經(jīng)驗公式計算結(jié)果進行了對比。

冰材料模型；破冰阻力；冰層；數(shù)值模擬

近十幾年來，隨著全球氣候變暖，北極冰層逐年融化，各國對連接歐洲、東亞和北美的北極航線的探索不斷深入，同時對北極地區(qū)油氣和礦產(chǎn)資源的關(guān)注不斷升溫。極地冰區(qū)的通航和科學(xué)考察、資源勘探需依靠破冰船開辟航道和破冰航行。因此，準(zhǔn)確模擬破冰船的破冰過程并預(yù)報該過程中破冰船所受的阻力，具有重要的理論和現(xiàn)實意義。

在國外，早期關(guān)于破冰阻力的研究主要采用經(jīng)驗公式方法。LINDQVIST[1]總結(jié)實船測量、模型實驗的結(jié)果，歸納出直接計算船舶在冰區(qū)航行阻力的經(jīng)驗公式。近些年來，大量學(xué)者使用數(shù)值方法研究船冰相互作用問題。SU等[2]在經(jīng)驗數(shù)據(jù)的基礎(chǔ)上，建立了研究船冰接觸時船體所受冰載荷的數(shù)值模型，對全局和局部冰載荷進行了系統(tǒng)分析；LIU等[3]采用非線性有限元軟件LS-DYNA分析了船舶和冰山碰撞的內(nèi)部機理，用單元的侵蝕模擬裂縫的增長；ZHOU等[4]建立數(shù)值模型，模擬了破冰船破冰過程中的動態(tài)冰載荷，并與試驗結(jié)果進行了對比；KIM等[5]進行了冰錐受壓實驗，并使用LS-DYNA進行數(shù)值仿真，驗證了一組冰材料模型。TAN等[6]使用離散元數(shù)值方法，模擬船在冰層中連續(xù)破冰的受力過程。

在國內(nèi)，數(shù)值模擬船舶破冰方面的研究方興未艾。何菲菲等[7]應(yīng)用商業(yè)軟件MSC.DYTRAN研究了破冰船的破冰載荷和破冰能力；張健等[8]應(yīng)用非線性有限元軟件研究了船-冰碰撞下球鼻艏結(jié)構(gòu)和船舶肩部的結(jié)構(gòu)動態(tài)響應(yīng)。李紫麟等[9]采用離散元模型對碎冰區(qū)浮冰與船舶的相互作用進行了數(shù)值研究。王健偉等[10]采用非線性有限元法建立散貨船與冰層的三維有限元模型，對散貨船與冰層的碰撞進行了數(shù)值模擬。上述方法對船舶破冰過程的數(shù)值模擬進行了一定嘗試，但尚缺乏對具體冰材料模型的驗證。

本文旨在使用非線性有限元軟件LS-DYNA建立冰材料數(shù)值模型，并結(jié)合實測數(shù)據(jù)進行驗證，進而采用該模型對破冰船在冰區(qū)航行時的連續(xù)破冰過程進行了數(shù)值模擬，對破冰阻力進行了預(yù)報。

1　冰錐受壓實驗的數(shù)值模擬

參考KIM的研究，冰材料選取LS-DYNA軟件中的可破碎泡沫模型(MAT_63: Crushable Foam Model)，使用見圖1的階梯狀正應(yīng)力-體積應(yīng)變曲線，并采用“截斷應(yīng)力”失效標(biāo)準(zhǔn)。材料屬性和失效標(biāo)準(zhǔn)參數(shù)(見表1)。冰錐的直徑為10 cm，錐角為120，有限元模型使用六面體單元，最大單元長度為5 mm。

圖1　體積應(yīng)變-正應(yīng)力曲線圖Fig.1　Volumetric strain-stress relation

表1　冰材料屬性和失效標(biāo)準(zhǔn)

為了驗證上述冰材料模型的準(zhǔn)確性，通過數(shù)值模擬鋼板擠壓圓錐狀冰體實驗，將計算所得的相互作用力與實驗結(jié)果進行比較。冰錐固定在實驗臺上，鋼板位于其上方并以恒定的速度向下運動，與冰錐發(fā)生碰撞。圖2給出了數(shù)值模擬該實驗的有限元模型。冰錐數(shù)值模型和實驗中所用樣品的形狀和尺寸保持一致，約束冰錐底部平面內(nèi)所有節(jié)點的運動。在實驗中，鋼板與冰體接觸的表面無明顯變形，因此，數(shù)值模擬中將鋼板設(shè)置為剛體，其材料參數(shù)見表2。

圖2　冰錐和鋼板的有限元模型Fig.2　FE model of ice cone and steel plate

表2　鋼材料模型參數(shù)

為使驗證所得結(jié)果適用性更全面，須考察冰體材料模型在不同應(yīng)變率下的適用情況。因此，分別令鋼板以1 mm/s和100 mm/s的速度向冰錐運動，模擬不同速度下碰撞的情境，并使用帶有侵蝕算法的點-面接觸模型(CONACT_ERODING_NODES_TO_SURFACE)模擬鋼板與冰層的接觸，使達到失效標(biāo)準(zhǔn)的冰體單元從模型中刪除。

圖3給出了數(shù)值模擬所得的不同速度下接觸力-位移曲線。從圖3可知，數(shù)值模擬所得的接觸力在整體走勢上與實驗值相符，二者誤差較小，不論是低速運動還是高速運動的計算值，都能夠反映接觸力上升的趨勢，同時也伴隨著一定程度的波動，這是由于鋼板在擠壓冰體時會產(chǎn)生冰體單元的破碎，導(dǎo)致瞬時接觸力突然下降。由此可以確定，上述驗證的冰材料與實際實驗的冰體力學(xué)性能相符，可進一步推廣。

另一方面，通過對比不同速度下的冰荷載過程可見，鋼板速度為100 mm/s時，冰錐與鋼板的接觸力比1 mm/s時更大。產(chǎn)生這一差異的原因是，鋼板高速運動時，對冰錐產(chǎn)生更強的沖擊力，二者碰撞更為劇烈，因此產(chǎn)生的接觸力更大。

圖3　接觸力-位移曲線Fig.3　Contact force-displacement curve

2　破冰船破冰的數(shù)值模型

基于上述冰錐受壓實驗與數(shù)值模擬，所選取的冰材料模型得到了驗證。將該模型進一步應(yīng)用于破冰船在冰層中連續(xù)破冰過程的數(shù)值模擬，并通過與經(jīng)驗公式計算值對比，進行進一步驗證。

2.1船與冰層數(shù)值模型

選擇某破冰船[11](A船)為計算對象，其主要船型參數(shù)見表3，橫剖面型線圖見圖9。該船船首為傳統(tǒng)的楔形船首，首部橫剖線呈V型。

破冰船首部在破冰過程中變形很小，可以忽略，因此在數(shù)值模擬中視船體外板為剛體。船體使用厚度為35 mm的殼單元進行模擬，材料屬性采用表2中的各項參數(shù)。對于船體的運動狀態(tài)，不考慮破冰過程對船體運動的影響，設(shè)定船體以恒定速度向前運動，約束除X方向(定速運動方向)外的其他5個自由度的運動。

表3　 A船主要參數(shù)

在船舶破冰過程中，當(dāng)冰層應(yīng)力達到一定值時，冰體單元將發(fā)生破壞。因此，不同于船體使用殼單元進行模擬，考慮到冰層的厚度以及對單元失效、侵蝕過程的仿真，采用體單元模擬冰層。冰層材料采用表1中各項參數(shù)。

由于研究的是船舶在無限冰區(qū)中航行的破冰過程。不同于在水面可自由運動的浮冰與流冰，無限大的冰層可視為位置固定。為模擬無限大冰層條件，取冰層沿船長方向為60 m，垂直于船長方向為30 m(船體縱中剖面左右各15 m)，與船舶接觸碰撞的邊界為自由端，遠場邊界為剛性固定。

2.2船與冰層接觸模型

在船舶破冰的過程中，冰層表面達到失效標(biāo)準(zhǔn)的單元須要不斷被侵蝕，隨之產(chǎn)生新的接觸面，繼續(xù)與船舶接觸。為保證在數(shù)值模擬船舶破冰過程中，當(dāng)單元的最大主應(yīng)力超過設(shè)定值后，單元失效，失效后的單元從模型中刪除，研究中使用帶有侵蝕算法的點-面接觸模型(CONTACT_ERODING_NODES_TO_SURFACE)模擬破冰船與冰層的接觸。在模擬兩種力學(xué)屬性差異較大的材料時，使用這一接觸模型可獲得較為精確可靠的模擬結(jié)果。冰層水線面與船體滿載吃水水線面的高度保持一致(見圖4)。

圖4　船舶破冰的有限元模型Fig.4　FE model of ship breaking level ice

3　船舶破冰過程數(shù)值模擬

3.1收斂性分析

為了獲得較可靠的數(shù)值模擬結(jié)果，選擇合適的冰體單元尺寸是至關(guān)重要的。為此首先進行收斂性分析，以為合理選擇冰體單元尺寸提供依據(jù)。

設(shè)置船速為3 m/s，冰層厚度為0.4 m，碰撞前船體與冰層之間的距離取為0.1 m；計算時間設(shè)為8.0 s。

破冰船與冰層碰撞過程中，二者之間的接觸力始終處于周期性的振蕩狀態(tài)中。因此，將一定周期內(nèi)船舶與冰層瞬時接觸力的平均值視為平均破冰阻力，并將其作為收斂性分析的考察參數(shù)。通過改變離散長度(相鄰節(jié)點間距離)，模擬不同離散長度下破冰船與冰層碰撞的平均破冰阻力，研究離散長度對平均破冰阻力的影響。

圖5給出了收斂性分析的結(jié)果。從圖5可知，當(dāng)冰體離散長度逐漸減小時，平均破冰阻力值逐漸下降，并且隨著離散長度的縮小而逐漸趨于收斂。綜合考慮破冰過程數(shù)值模擬所需時間和模擬精度等因素，選取冰層離散長度為125 mm進行后續(xù)數(shù)值模擬。

圖5　平均破冰阻力收斂性Fig.5　Convergence of mean icebreaking resistance

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

圖6給出了冰體離散長度為125 mm時，破冰船在破冰過程中所受到的破冰阻力時歷曲線。從圖6可知，在0～5 s階段，破冰船以3 m/s的恒定速度駛向冰層，由于在此過程中，船體不斷擠壓冰層單元，且船體與冰層擠壓接觸的面積不斷增大，破冰阻力曲線整體上隨時間的推移振蕩上升。大約5 s之后，破冰船與冰層接觸的表面積基本保持不變，破冰船進入連續(xù)破冰階段。在這一階段，破冰船單位時間內(nèi)擠壓、侵蝕的冰單元質(zhì)量趨于穩(wěn)定，使得船與冰的平均破冰阻力也趨于穩(wěn)定。此外，在船舶與冰層碰撞的整個時歷過程中，始終伴隨著冰體單元的失效，因此會產(chǎn)生接觸力卸載現(xiàn)象，導(dǎo)致破冰阻力時歷曲線始終處于振蕩狀態(tài)。

圖6　破冰阻力時歷曲線Fig.6　Time history of icebreaking resistance

圖7　破冰阻力時歷曲線與LINDQVIST破冰阻力比較(0.4 m厚冰層)Fig.7　Comparison of simulated time history of icebreaking resistance and LINDQVIST’s icebreaking resistance in 0.4 m thick ice

基于上述分析，選擇5 s之后的接觸力曲線，計算船舶破冰過程中的平均阻力，并將其與LINDQVIST提出的公式計算所得的平均破冰阻力進行對比。圖7給出了7～8 s時的船-冰接觸力和LINDQVIST提出的公式所得的平均破冰阻力，作為5 s之后破冰船在破冰過程中的受力例子。從圖7可知，在連續(xù)破冰階段，模擬所得的破冰阻力始終在LINDQVIST破冰阻力附近區(qū)間內(nèi)振蕩，數(shù)值模擬的平均破冰阻力與LINDQVIST破冰阻力二者基本吻合。

3.3冰層厚度對船舶破冰過程的影響

鑒于破冰船在航行過程中需在不同厚度的冰層中破冰，對破冰船與不同厚度冰層碰撞的情境進行了數(shù)值模擬。參考KUJALA[12]的長期實測統(tǒng)計數(shù)據(jù)，冰層的平均厚度約為0.4 m，標(biāo)準(zhǔn)差約為0.15 m；因此在模擬時考察冰層厚度為0.25～0.55 m時，破冰船所受的平均破冰阻力。

圖8給出了不同冰層厚度下數(shù)值模擬的平均破冰阻力，并將其與LINDQVIST經(jīng)驗公式結(jié)果進行對比。從圖8可知，模擬所得的平均破冰阻力總體上隨著冰厚的增加而上升，且與LINDQVIST破冰阻力曲線相比，數(shù)值模擬所得的平均破冰阻力基本在LINDQVIST破冰阻力曲線附近區(qū)間，二者之間差別較小。

圖8　A,B船不同冰層厚度下平均破冰阻力和LINDQVIST結(jié)果比較Fig.8　Comparison of the mean icebreaking resistance of Ship A and Ship B with LINDQVIST’s results in ice with different thickness

為進一步驗證冰材料模型及數(shù)值模擬方法的適用性，對另一艘同等船長、型寬、吃水，不同船首形狀的破冰船(B船)進行數(shù)值模擬，考察該船在連續(xù)破冰過程中的破冰阻力。圖9給出了該船首部的橫剖面型線圖。該船船首為勺型，首部橫剖線接近圓形，與楔形船首存在明顯差異。破冰船在連續(xù)破冰過程中與冰層發(fā)生碰撞的主要區(qū)域為船首，因此，不同船首形狀的破冰船可能產(chǎn)生不同的阻力響應(yīng)，對其進行模擬可進一步驗證所采用的數(shù)值模擬方法的適用性。

B船速度同樣取為3 m/s，分別在0.25～0.55 m冰層厚度情況下進行數(shù)值模擬。圖8給出了B船在不同冰層厚度下數(shù)值模擬的平均破冰阻力，并將其與所對應(yīng)的LINDQVIST破冰阻力曲線進行對比。從圖8可知，與A船相比，由于B船的首柱傾角較大，使得平均破冰阻力相對較小?？傮w上，模擬所得的A船、B船平均破冰阻力與LINDQVIST破冰阻力曲線基本吻合。因此，采用的冰材料模型和數(shù)值模擬方法可推廣應(yīng)用于不同船首形狀的破冰船的破冰過程模擬。

圖9　兩艘破冰船的首部橫剖線圖Fig.9　Bow profile of 2 icebreakers

4　結(jié)　論

采用有限元數(shù)值模擬方法，通過與實驗數(shù)據(jù)比對，驗證了一組冰材料模型，進而將其應(yīng)用于破冰船無限冰區(qū)破冰過程的數(shù)值模擬研究，并將模擬所得的破冰阻力值與經(jīng)驗公式計算的平均阻力值進行對比，以進一步驗證冰材料模型的準(zhǔn)確性和適用性。結(jié)果表明，該冰材料模型及數(shù)值模擬方法能夠比較準(zhǔn)確地模擬破冰船在不同厚度冰層中連續(xù)破冰的過程，且適用于具有不同形狀船首的破冰船。

[1] LINDQVIST G. A straightforward method for calculation of ice resistance of ships[C]//The 10th International Conference on Port and Ocean Engineering under Arctic Conditions. Lulea, Sweden, 1989, 2: 722-735.

[2] SU B,RISKA K, MOAN T. A numerical method for the prediction of ship performance in level ice[J]. Cold Regions Science and Technology, 2010, 60(3): 177-188.

[3] LIU Z, AMDAHL J, L?SET S. Integrated numerical analysis of an iceberg collision with a foreship structure[J]. Marine Structures, 2011, 24(4): 377-395.

[4] ZHOU L,RISKA K, MOAN T, et al. Numerical modeling of ice loads on an icebreaking tanker: comparing simulations with model tests[J]. Cold Regions Science and Technology, 2013, 87: 33-46.

[5] KIM H. Simulation of compressive ‘cone-shaped’ ice specimen experiments using LS-DYNA?[C]//13th International LS-DYNA Users Conference. Detroit, MI, 2014.

[6] TAN X,RISKA K, MOAN T. Performance simulation of a dual-direction ship in level ice[J]. Journal of Ship Research, 2014, 58(3): 168-181.

[7] 何菲菲. 破冰船破冰載荷與破冰能力計算方法研究[D].哈爾濱：哈爾濱工程大學(xué), 2010.

[8] 張健, 張淼溶, 萬正權(quán), 等. 冰材料模型在船一冰碰撞結(jié)構(gòu)響應(yīng)數(shù)值仿真中的應(yīng)用研究[J]. 中國造船,2013,54(4):100-108.

ZHANG Jian, ZHANG Miaorong, WAN Zhengquan, et al.Research on ice material model applied in numerical simulation of ship structure response under iceberg collision[J]. Shipbuilding of China, 2013, 54(4):100-108.

[9] 李紫麟, 劉煜, 孫珊珊, 等. 船舶在碎冰區(qū)航行的離散元模型及冰載荷分析[J]. 力學(xué)學(xué)報,2013,45(6):868-877.

LI Zilin, LIU Yu, SUN Shanshan, et al. Analysis of ship maneuvering performances and ice loads on ship hull with discrete element model in broken-ice fields[J]. Chinese Journal of Theoretical and Applied Mechanics,2013, 45(6):868-877.

[10] 王健偉, 鄒早建. 基于非線性有限元法的船舶-冰層碰撞結(jié)構(gòu)響應(yīng)研究[J]. 振動與沖擊, 2015,34(23):125-130.

WANG Jianwei, ZOU Zaojian. Study on structure response of ship collision with level ice based on nonlinear finite element method[J].Journal of Vibration and Shock, 2015, 34(23): 125-130.

[11] YAMAGUCHI H, SUZUKI Y, UEMURA O, et al. Influence of bow shape on icebreaking resistance in low speed range[C]//Proceedings Vol. IV of the 16th International Conference on Offshore Mechanics and Arctic Engineering/14th International Conference on Port and Ocean Engineering under Arctic Conditions (Joint OMAE/ POAC Conference). Yokohama, Japan, 1997: 51-62.

[12] KUJALA P. On the statistics of ice loads on ship hull in the baltic[R]. Mechanical Engineering Series No. 116. Finland: Ship Laboratory, Helsinki University of Technology,1994.

Numerical simulation of the continuous icebreaking by an icebreaker in level ice

REN Yizhou1, ZOU Zaojian1,2

(1. School of Naval Architecture, Ocean and Civil Engineering, Shanghai Jiao Tong University, Shanghai 200240, China;2. State Key Laboratory of Ocean Engineering, Shanghai Jiao Tong University, Shanghai 200240, China)

A numerical model of ice material was introduced for simulating the continuous icebreaking by an icebreaker in level ice. The numerical model was validated by comparing the numerical simulation results with the measured data in compressive cone-shaped ice experiments. The material model was applied to the numerical simulation of collision between the icebreaker and level ice in infinite ice field, and the icebreaking resistance was calculated. The results of numerical simulation and the icebreaking resistance calculated with different shape of bow and different ice thickness were compared with those by empirical formula.

ice material model; icebreaking resistance; level ice; numerical simulation

2015-08-20修改稿收到日期：2015-09-16

任奕舟 男，碩士生，1990年生

鄒早建 男，博士，教授，1956年生

U661.4

A DOI:10.13465/j.cnki.jvs.2016.14.034