极地某构型船首与浮冰碰撞的数值模拟研究

刘 昕 傅昱晓 蒋武杰 倪宝玉＊

（1. 中国人民解放军92942 部队 北京 100161；2. 哈尔滨工程大学 船舶工程学院 哈尔滨 150001；3. 上海交通大学 船舶海洋与建筑工程学院 上海 200030）

0 引 言

在地球上，冰川的覆盖面积约为1 600 万km2,约占地球表面积的3.1%[1]。随着全球气候变暖以及冰川逐渐消融，北极的战略意义和商业价值正被逐步发掘。与其他海域相比，极地地区由冰川融化而生成的浮冰由于面积更大、冰层更厚，因此与船舶发生碰撞时的危险性也大大提高。这些大尺寸浮冰融化、相互碰撞，又可能会生成更多数量的小块浮冰。这些浮冰主要集中在冰盖区与开阔水域中间的边缘冰区内，这就使船舶在极地航行时，与浮冰碰撞的概率大大增加。因此，模拟冰材料、构建船体模型以及分析冰船碰撞的动态过程，逐渐成为极地领域热门的研究方向。

基于冰材料物理与力学特性，学者们提出了各种可用于数值模拟的冰材料模型。SCHNEIDER 等[2]依据海冰的失效模式建立了1 种几何网格方法，用于计算海冰与船体的相互作用。SAWAMURA 等[3]基于瞬态载荷与浮冰作用的过程，得到了浮冰的弯曲特性。EHLERS 等[4]通过把数值模拟结果与早期进行的一系列海冰四点弯曲实验进行对比，发现数值模拟结果在受力、破坏时间和位移方面与试验结果吻合较好，验证了塑性材料可用于模拟冰船碰撞时的海冰材料。马国宝[5]通过讨论海冰性能的影响因素，确定海冰材料在进行数值模拟时应采用弹塑性材料模型，同时模拟圆台冰撞击刚性墙的过程，将得到的压力-面积曲线与ISO 给定曲线进行对比，验证了模拟冰材料的可行性。

针对船舶与浮冰的作用，LUBBAD 等[6]通过数值模拟方法与实验方法，研究了冰船碰撞问题，分析船舷承载能力、结构变形、结构稳定性以及板架的承载能力。何菲菲[7]采用理论推导和数值仿真方法，基于非线性有限元法，对破冰船船体运动响应以及破冰能力进行了探索。韩月[8]通过模拟破冰船与海冰的相互作用，讨论了不同形状船首的破冰性能，并基于数值模拟方法预报了船体的运动响应和破冰载荷。张健等[9]提出了连续破冰模式下的船首部冰载荷分析方法，并与数值仿真结果进行对比。

本文首先通过整理海冰的物理力学特性，确定模拟海冰材料的参数，运用LS-DYNA 软件中的各向同性材料来模拟冰材料；并设置了冰锥与钢板撞击过程，通过与ISO 规范提出的P-A 曲线[10]对比，验证模拟冰材料的可行性；然后基于某极地船型原型构建船首模型，通过模拟船首以不同速度与浮冰碰撞的工况，分析船首速度对于冰船碰撞的影响。

1 模拟海冰材料及其有效性验证

在研究船舶与海冰碰撞的过程中，LS-DYNA是较常用的软件。本章将根据海冰材料的物理力学特性，利用LS-DYNA 设置海冰材料的相关参数，并进行冰材料模型的验证。

1.1 模拟海冰材料

通过对海冰材料的整理，得到其物理力学特性如表1[1]所示。

表1 海冰的物理力学特性统计值

海冰是具有各项异性的弹塑性材料，但考虑到模拟各向异性材料耗费的计算资源较大，且对于研究船舶与浮冰碰撞过程来说，各向异性与各向同性材料的差异性可忽略。因此，在研究浮冰与海上大型船舶相互作用时，可采用LS-DYNA 中的13 号材料（*MAT_ISOTROPIC_ELASTIC_FAILURE）来模拟海冰材料，并输入参数值：密度为910 kg/m3、剪切模量为2.2 GPa、屈服应力为2.12 MPa、塑性硬化模量为4.26 GPa、体积模量为5.26 GPa、塑性失效应变值为0.35。

该材料是带有塑性应变失效准则的各项同性弹性断裂失效模型，即当材料单元的应力或应变达到最大许用值时，单元会失效且无法再承载应力。

1.2 海冰材料的有效性验证

为验证模拟冰材料的合理性，本文通过对刚性板撞击锥形冰体进行模拟，将模拟结果与ISO 收录的P-A 理论标准曲线[9]进行对比验证。其中，锥形冰体和刚性板的模型如图1 所示。

图1 锥形冰体与刚性板模型

在图1 中，上部圆锥半径为0.5 m、高为0.289 m，下部圆柱半径为0.5 m、高为0.1 m。刚性板为1.2 m×1.2 m 的正方形板，厚度为0.1 m。锥形冰的网格尺寸为0.04 ~ 0.05 m，刚性板网格尺寸为0.03 m。锥形冰前端距刚性板表面距离为0.02 m，刚性板以0.1 m/s 的速度向锥形冰撞去。刚性板材料选用了LS-DYNA 材料库中第20 号材料“020-RIGID”。该材料不会变形且没有应力、应变等，相当于固定的刚体，因此在计算中，仅需采集冰锥的受力情况。为保证碰撞的连续发生，锥形冰底部所有节点的6个自由度运动均被约束。

在刚性板向锥形冰撞击的过程中，定义刚性板与锥形冰的接触为侵蚀接触，即冰单元体失效时,程序会自动搜索下一个接触节点。其中，从面定义为浮冰模型，主面定义为刚性板模型。计算过程中,仅包含从面上的作用力。将静摩擦系数设置为1.5,动摩擦系数设置为1.0。在对称平面的选择上，不包括有正常边界约束的面；在内部节点的选择上,储存器的分配使腐蚀接触可以发生，否则，单元格被侵蚀后将没有接触。在对固体元素进行相邻材料的处理上，设定实体元素面位于材料子集的边界上。图2 为刚性板撞击锥形冰体的应力云图。

图2 刚性板撞击锥形冰体的应力云图

由图2 可见，在刚性板与锥形冰碰撞的过程中，随着刚性板的运动，锥形冰逐渐断裂，冰单元失效且失效面不平整。从时间分别为0.125 s、0.15 s、0.17 s 和0.19 s 的锥形冰体应力云图可以看出，接触面的失效应力分布不均匀，存在明显的高应力区和低应力区。其中，最大应力分别为12.950 MPa、38.785 MPa、36.400 MPa 和39.701 MPa。

为验证冰材料的合理性，需要将刚性板-锥形冰体碰撞的压强-面积曲线与ISO 标准理论曲线进行对比。国际标准化委员会在ISO/CD 19906（2010）中，推荐了冰山碰撞载荷下造成严重事故时的冰体材料P-A 曲线[9]，这是1 条反比例曲线。

本文所模拟的刚性板与锥形冰碰撞中，接触面的压强与面积关系式为：

式中：P为接触面的压强，MPa；A为碰撞接触面积，m2。

刚性板与锥形冰碰撞的压强-面积曲线和ISO标准理论曲线的对比情况如图3 所示。

图3 压强-面积曲线对比情况

由图3 可见，初始阶段压强值都比较大，随着碰撞接触面积的增加，压强值迅速减小并趋于平稳。当面积达到0.4 m2时，本文模拟冰材料对应的曲线产生小幅波动。但整体来看，两者吻合度较高，所以本文可以选取该材料模型进行极地海冰与船舶相互作用的数值模拟。

2 船冰碰撞动态响应的工况分析

2.1 典型工况分析与主要船首部件受力分析

本文所用的船体模型是基于某极地船型建立的，且该船型非破冰船型，在极区航道航行时，会尽量避免与大块浮冰、冰脊碰撞，浮冰厚度主要参考较薄的一年冰或多年冰来选取。本文设置的冰厚为0.6 m，浮冰形状为正方形，边长为10 m。遇到较大尺寸的浮冰时，建议绕行。综合考虑计算精度与计算性能，冰单元的大小取200 mm×200 mm×100 mm。为模拟浮冰在海上漂浮的状态，设置浮冰模型初始状态为静止，边界条件为自由边界，可在6 个自由度上进行自由运动。

考虑到整船结构建模计算量巨大，且船冰碰撞主要发生在船首区域，因此只选择船首进行建模，如图4 所示，网格大小为300 mm×400 mm。船首模型除外部壳体，还包括甲板、纵骨、纵桁、横舱壁、肋骨、肋板和舭部等24 个壳部件和47 个梁部件，船首与浮冰距离为0.5 m。

图4 船首有限元模型

模型船首材料采用高强度钢，当应力超出材料的弹性极限时，它便由弹性状态变为塑性状态。本文用LS-DYNA 材料库中的“001_ELASTIC”来模拟钢材料。该材料是一种各向同性的弹性材料，当所受应变超出其承受范围时会导致结构变形失效，对应参数选择为杨氏模量2 060 GPa、密度为7 850 kg/m3、泊松比为0.3。

选择船首速度6 kn 为典型工况，整个碰撞过程持续约5 s，各时间点船首与浮冰碰撞的应力云图及浮冰模型的应力云图如下页图5 所示。

可以看出，在船首与浮冰发生第1 次碰撞时，船首的应力开始由首端向后传递，但船首端部与浮冰接触的区域仍然是应力最大区域。浮冰也因与船首碰撞而发生破坏，但此时碰撞速度较慢（即加载速率较小），浮冰模型表面虽有明显裂纹，但裂纹并未扩张，且碰撞过程中，船首的一部分动能转化为浮冰的动能和内能。所以在第1 次碰撞后，浮冰有了速度，并发生位移和翻转。

在t=3 s 附近，船首与浮冰发生第2 次碰撞。浮冰由于第1 次碰撞时已具备一定速度，所以在这次碰撞中，船首和浮冰的碰撞应力相对较小，浮冰模型只发生了局部破坏，并没有明显的裂纹扩张，船首所受应力也有所减少。

对于壳部件而言，船首底部甲板是碰撞区域主要的支撑部件，船首端部结构是与浮冰直接接触的部件。在船首区域，船底板距离碰撞区域很近，而且也是承受总纵弯曲的主要部件之一。所以，壳部件主要研究这3 个部件。目标船型为纵骨架式船，总纵弯曲主要由纵向梁分担，所以选取的梁部件均为纵向梁，包括上层建筑主要纵桁、主甲板纵骨和船底甲板纵骨，如下页图6 所示。

图6 船首主要部件选取

图7 给出了船首上述各主要受力部件所包含单元，在每一时间步所受应力的最大值组成的曲线，表2 则给出了此过程中上述各部件的最大应力及对应时间。

图7 船首主要部件上单元的最大应力

表2 各部件最大应力及对应时间

可以看出：纵向梁部件中承受应力较大的是上层建筑主要纵桁和船底甲板纵骨，其应力最大值分别达到了26.95 MPa 和32.82 MPa；对于壳部件而言，应力在达到最大值之后波动减小，其中承受应力最大的是直接与浮冰接触的船首端部构件，在0.4 s 达到最大值105.75 MPa。

由于各部件相对碰撞区域的位置不同，因此应力的峰值不同，并且距碰撞区域越远的构件，应力峰值出现的时间也相对较晚。各部件结构不同，纵向梁部件碰撞峰值后的应力波动情况相对较平缓，而壳部件则呈现间接脉冲式分布。

2.2 不同速度下船首受力情况对比

破冰船在浮冰水域航行时，既要完成破冰，还会受到碎冰和小型浮冰的阻力，因此航行速度较慢。为研究不同速度对船首与浮冰碰撞模拟结果的影响，本文选取船首航行方向（x轴方向）速度分别为5 kn（工况一）、6 kn（工况二）和7 kn（工况三），其余方向的运动都被约束。浮冰的速度设置为0 m/s，6 个自由度的运动均为自由运动，尺寸保持10 m×10 m×0.6 m 不变。

图8 所示为工况一、工况二和工况三对应的船首所受合力-时间图像。可以看出，上述各工况下均出现2 个较明显的波峰，说明船首与浮冰发生了2 次较剧烈的碰撞（参见图5 工况）。

图8 所选工况下的船首所受合力-时间图

随着速度增大，波峰的出现逐渐提前，第2 次碰撞的波峰情况尤其明显。通过比较2 个峰值，可以看出船首速度对于第1 次碰撞的冰力影响更为明显。

由于冰船相距较近，浮冰初始状态是静止的，因此第1 次波峰出现的时间比较接近，且随着船速增大，各工况下对应的船首所受合力的最大值随之增大，如表3 所示。对于冰船第2 次碰撞，随着船速增大，第2 次的峰值也随之增大，且出现时间随之提前。

表3 所选工况下的合力极大值

图9 给出了各工况下的船首与浮冰碰撞的结果。

图9 工况一、工况二和工况三的碰撞结果

由图9 可见，速度对于浮冰失效方式的影响，在工况一（5 kn）和工况二（6 kn）时，海冰的加载速率较小，浮冰模型表面的裂纹没有明显扩张；在工况三（7 kn）时，随着船速增大，海冰的加载速率变大，浮冰模型的裂纹明显扩张，最终变为两半，出现了劈裂破坏。

3 结 语

本文根据冰材料常见的物理及力学特性，基于LS-DYNA 中的各向同性材料设置了冰材料模型，模拟锥形冰撞击钢板的过程，通过与ISO 规范提出的压力-面积曲线对比，验证了模拟冰材料的可行性。

文中设置3 组速度分别为5 kn、6 kn 和7 kn时的工况，分析船首与浮冰碰撞过程。船首速度决定了外力对浮冰的加载速度。船速越大、冰载荷也越大，速度为7 kn 时的合力达到2.1 MN，且碰撞后浮冰会发生翻转。随着船首前进，冰与船会发生二次碰撞。由于本船型船首瘦削，所以浮冰容易产生劈裂破坏。在碰撞过程中，船首端部构件承受的应力最大，建议根据实际航行工况决定是否对该部分构件进行加固。