破冰船在冰层中连续破冰过程的数值模拟

任奕舟， 邹早建,2

(1. 上海交通大学 船舶海洋与建筑工程学院,上海　200240； 2.上海交通大学 海洋工程国家重点实验室,上海　200240)



破冰船在冰层中连续破冰过程的数值模拟

任奕舟1， 邹早建1,2

(1. 上海交通大学 船舶海洋与建筑工程学院,上海200240； 2.上海交通大学 海洋工程国家重点实验室,上海200240)

介绍了一种用于模拟破冰船在冰层中连续破冰的冰材料数值模型。通过与冰锥受压实验数据进行对比，对该数值模型进行了验证。将该模型应用于破冰船在无限冰区与冰层碰撞的数值模拟，对破冰船的破冰阻力进行了计算，并将不同船首、不同冰层厚度下计算所得的破冰阻力与经验公式计算结果进行了对比。

冰材料模型；破冰阻力；冰层；数值模拟

近十几年来，随着全球气候变暖，北极冰层逐年融化，各国对连接欧洲、东亚和北美的北极航线的探索不断深入，同时对北极地区油气和矿产资源的关注不断升温。极地冰区的通航和科学考察、资源勘探需依靠破冰船开辟航道和破冰航行。因此，准确模拟破冰船的破冰过程并预报该过程中破冰船所受的阻力，具有重要的理论和现实意义。

在国外，早期关于破冰阻力的研究主要采用经验公式方法。LINDQVIST[1]总结实船测量、模型实验的结果，归纳出直接计算船舶在冰区航行阻力的经验公式。近些年来，大量学者使用数值方法研究船冰相互作用问题。SU等[2]在经验数据的基础上，建立了研究船冰接触时船体所受冰载荷的数值模型，对全局和局部冰载荷进行了系统分析；LIU等[3]采用非线性有限元软件LS-DYNA分析了船舶和冰山碰撞的内部机理，用单元的侵蚀模拟裂缝的增长；ZHOU等[4]建立数值模型，模拟了破冰船破冰过程中的动态冰载荷，并与试验结果进行了对比；KIM等[5]进行了冰锥受压实验，并使用LS-DYNA进行数值仿真，验证了一组冰材料模型。TAN等[6]使用离散元数值方法，模拟船在冰层中连续破冰的受力过程。

在国内，数值模拟船舶破冰方面的研究方兴未艾。何菲菲等[7]应用商业软件MSC.DYTRAN研究了破冰船的破冰载荷和破冰能力；张健等[8]应用非线性有限元软件研究了船-冰碰撞下球鼻艏结构和船舶肩部的结构动态响应。李紫麟等[9]采用离散元模型对碎冰区浮冰与船舶的相互作用进行了数值研究。王健伟等[10]采用非线性有限元法建立散货船与冰层的三维有限元模型，对散货船与冰层的碰撞进行了数值模拟。上述方法对船舶破冰过程的数值模拟进行了一定尝试，但尚缺乏对具体冰材料模型的验证。

本文旨在使用非线性有限元软件LS-DYNA建立冰材料数值模型，并结合实测数据进行验证，进而采用该模型对破冰船在冰区航行时的连续破冰过程进行了数值模拟，对破冰阻力进行了预报。

1　冰锥受压实验的数值模拟

参考KIM的研究，冰材料选取LS-DYNA软件中的可破碎泡沫模型(MAT_63: Crushable Foam Model)，使用见图1的阶梯状正应力-体积应变曲线，并采用“截断应力”失效标准。材料属性和失效标准参数(见表1)。冰锥的直径为10 cm，锥角为120，有限元模型使用六面体单元，最大单元长度为5 mm。

图1　体积应变-正应力曲线图Fig.1　Volumetric strain-stress relation

表1　冰材料属性和失效标准

为了验证上述冰材料模型的准确性，通过数值模拟钢板挤压圆锥状冰体实验，将计算所得的相互作用力与实验结果进行比较。冰锥固定在实验台上，钢板位于其上方并以恒定的速度向下运动，与冰锥发生碰撞。图2给出了数值模拟该实验的有限元模型。冰锥数值模型和实验中所用样品的形状和尺寸保持一致，约束冰锥底部平面内所有节点的运动。在实验中，钢板与冰体接触的表面无明显变形，因此，数值模拟中将钢板设置为刚体，其材料参数见表2。

图2　冰锥和钢板的有限元模型Fig.2　FE model of ice cone and steel plate

表2　钢材料模型参数

为使验证所得结果适用性更全面，须考察冰体材料模型在不同应变率下的适用情况。因此，分别令钢板以1 mm/s和100 mm/s的速度向冰锥运动，模拟不同速度下碰撞的情境，并使用带有侵蚀算法的点-面接触模型(CONACT_ERODING_NODES_TO_SURFACE)模拟钢板与冰层的接触，使达到失效标准的冰体单元从模型中删除。

图3给出了数值模拟所得的不同速度下接触力-位移曲线。从图3可知，数值模拟所得的接触力在整体走势上与实验值相符，二者误差较小，不论是低速运动还是高速运动的计算值，都能够反映接触力上升的趋势，同时也伴随着一定程度的波动，这是由于钢板在挤压冰体时会产生冰体单元的破碎，导致瞬时接触力突然下降。由此可以确定，上述验证的冰材料与实际实验的冰体力学性能相符，可进一步推广。

另一方面，通过对比不同速度下的冰荷载过程可见，钢板速度为100 mm/s时，冰锥与钢板的接触力比1 mm/s时更大。产生这一差异的原因是，钢板高速运动时，对冰锥产生更强的冲击力，二者碰撞更为剧烈，因此产生的接触力更大。

图3　接触力-位移曲线Fig.3　Contact force-displacement curve

2　破冰船破冰的数值模型

基于上述冰锥受压实验与数值模拟，所选取的冰材料模型得到了验证。将该模型进一步应用于破冰船在冰层中连续破冰过程的数值模拟，并通过与经验公式计算值对比，进行进一步验证。

2.1船与冰层数值模型

选择某破冰船[11](A船)为计算对象，其主要船型参数见表3，横剖面型线图见图9。该船船首为传统的楔形船首，首部横剖线呈V型。

破冰船首部在破冰过程中变形很小，可以忽略，因此在数值模拟中视船体外板为刚体。船体使用厚度为35 mm的壳单元进行模拟，材料属性采用表2中的各项参数。对于船体的运动状态，不考虑破冰过程对船体运动的影响，设定船体以恒定速度向前运动，约束除X方向(定速运动方向)外的其他5个自由度的运动。

表3　 A船主要参数

在船舶破冰过程中，当冰层应力达到一定值时，冰体单元将发生破坏。因此，不同于船体使用壳单元进行模拟，考虑到冰层的厚度以及对单元失效、侵蚀过程的仿真，采用体单元模拟冰层。冰层材料采用表1中各项参数。

由于研究的是船舶在无限冰区中航行的破冰过程。不同于在水面可自由运动的浮冰与流冰，无限大的冰层可视为位置固定。为模拟无限大冰层条件，取冰层沿船长方向为60 m，垂直于船长方向为30 m(船体纵中剖面左右各15 m)，与船舶接触碰撞的边界为自由端，远场边界为刚性固定。

2.2船与冰层接触模型

在船舶破冰的过程中，冰层表面达到失效标准的单元须要不断被侵蚀，随之产生新的接触面，继续与船舶接触。为保证在数值模拟船舶破冰过程中，当单元的最大主应力超过设定值后，单元失效，失效后的单元从模型中删除，研究中使用带有侵蚀算法的点-面接触模型(CONTACT_ERODING_NODES_TO_SURFACE)模拟破冰船与冰层的接触。在模拟两种力学属性差异较大的材料时，使用这一接触模型可获得较为精确可靠的模拟结果。冰层水线面与船体满载吃水水线面的高度保持一致(见图4)。

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

3　船舶破冰过程数值模拟

3.1收敛性分析

为了获得较可靠的数值模拟结果，选择合适的冰体单元尺寸是至关重要的。为此首先进行收敛性分析，以为合理选择冰体单元尺寸提供依据。

设置船速为3 m/s，冰层厚度为0.4 m，碰撞前船体与冰层之间的距离取为0.1 m；计算时间设为8.0 s。

破冰船与冰层碰撞过程中，二者之间的接触力始终处于周期性的振荡状态中。因此，将一定周期内船舶与冰层瞬时接触力的平均值视为平均破冰阻力，并将其作为收敛性分析的考察参数。通过改变离散长度(相邻节点间距离)，模拟不同离散长度下破冰船与冰层碰撞的平均破冰阻力，研究离散长度对平均破冰阻力的影响。

图5给出了收敛性分析的结果。从图5可知，当冰体离散长度逐渐减小时，平均破冰阻力值逐渐下降，并且随着离散长度的缩小而逐渐趋于收敛。综合考虑破冰过程数值模拟所需时间和模拟精度等因素，选取冰层离散长度为125 mm进行后续数值模拟。

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

3.2数值模拟结果

图6给出了冰体离散长度为125 mm时，破冰船在破冰过程中所受到的破冰阻力时历曲线。从图6可知，在0～5 s阶段，破冰船以3 m/s的恒定速度驶向冰层，由于在此过程中，船体不断挤压冰层单元，且船体与冰层挤压接触的面积不断增大，破冰阻力曲线整体上随时间的推移振荡上升。大约5 s之后，破冰船与冰层接触的表面积基本保持不变，破冰船进入连续破冰阶段。在这一阶段，破冰船单位时间内挤压、侵蚀的冰单元质量趋于稳定，使得船与冰的平均破冰阻力也趋于稳定。此外，在船舶与冰层碰撞的整个时历过程中，始终伴随着冰体单元的失效，因此会产生接触力卸载现象，导致破冰阻力时历曲线始终处于振荡状态。

图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可知，在连续破冰阶段，模拟所得的破冰阻力始终在LINDQVIST破冰阻力附近区间内振荡，数值模拟的平均破冰阻力与LINDQVIST破冰阻力二者基本吻合。

3.3冰层厚度对船舶破冰过程的影响

鉴于破冰船在航行过程中需在不同厚度的冰层中破冰，对破冰船与不同厚度冰层碰撞的情境进行了数值模拟。参考KUJALA[12]的长期实测统计数据，冰层的平均厚度约为0.4 m，标准差约为0.15 m；因此在模拟时考察冰层厚度为0.25～0.55 m时，破冰船所受的平均破冰阻力。

图8给出了不同冰层厚度下数值模拟的平均破冰阻力，并将其与LINDQVIST经验公式结果进行对比。从图8可知，模拟所得的平均破冰阻力总体上随着冰厚的增加而上升，且与LINDQVIST破冰阻力曲线相比，数值模拟所得的平均破冰阻力基本在LINDQVIST破冰阻力曲线附近区间，二者之间差别较小。

图8　A,B船不同冰层厚度下平均破冰阻力和LINDQVIST结果比较Fig.8　Comparison of the mean icebreaking resistance of Ship A and Ship B with LINDQVIST’s results in ice with different thickness

为进一步验证冰材料模型及数值模拟方法的适用性，对另一艘同等船长、型宽、吃水，不同船首形状的破冰船(B船)进行数值模拟，考察该船在连续破冰过程中的破冰阻力。图9给出了该船首部的横剖面型线图。该船船首为勺型，首部横剖线接近圆形，与楔形船首存在明显差异。破冰船在连续破冰过程中与冰层发生碰撞的主要区域为船首，因此，不同船首形状的破冰船可能产生不同的阻力响应，对其进行模拟可进一步验证所采用的数值模拟方法的适用性。

B船速度同样取为3 m/s，分别在0.25～0.55 m冰层厚度情况下进行数值模拟。图8给出了B船在不同冰层厚度下数值模拟的平均破冰阻力，并将其与所对应的LINDQVIST破冰阻力曲线进行对比。从图8可知，与A船相比，由于B船的首柱倾角较大，使得平均破冰阻力相对较小。总体上，模拟所得的A船、B船平均破冰阻力与LINDQVIST破冰阻力曲线基本吻合。因此，采用的冰材料模型和数值模拟方法可推广应用于不同船首形状的破冰船的破冰过程模拟。

图9　两艘破冰船的首部横剖线图Fig.9　Bow profile of 2 icebreakers

4　结　论

采用有限元数值模拟方法，通过与实验数据比对，验证了一组冰材料模型，进而将其应用于破冰船无限冰区破冰过程的数值模拟研究，并将模拟所得的破冰阻力值与经验公式计算的平均阻力值进行对比，以进一步验证冰材料模型的准确性和适用性。结果表明，该冰材料模型及数值模拟方法能够比较准确地模拟破冰船在不同厚度冰层中连续破冰的过程，且适用于具有不同形状船首的破冰船。

[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].哈尔滨：哈尔滨工程大学, 2010.

[8] 张健, 张淼溶, 万正权, 等. 冰材料模型在船一冰碰撞结构响应数值仿真中的应用研究[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] 李紫麟, 刘煜, 孙珊珊, 等. 船舶在碎冰区航行的离散元模型及冰载荷分析[J]. 力学学报,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] 王健伟, 邹早建. 基于非线性有限元法的船舶-冰层碰撞结构响应研究[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