船-冰碰撞下冰材料有限元数值方法研究进展

蔡 伟，朱 凌，3，毕璐泽

（1.高性能船舶技术教育部重点实验室（武汉理工大学），武汉 430063；2.武汉理工大学交通学院，武汉 430063；3.高新船舶与深海开发装备协同创新中心，武汉 430063）

0 引 言

极地资源的开发、航道运输以及科考成为当今世界各国关注的热点话题。随着全球气候变暖，北极冰川加速融化，极地生产运输活动逐渐增多，致使浮冰，特别是冰山逐渐增多，这给航行于冰区的船舶的结构安全带来了巨大的挑战。极地地区海面上漂浮的冰体主要分为淡水冰和咸水冰，其中淡水冰主要是冰山冰，或者是一些从冰川分离而出的小型冰山。极地船舶在冰区航行时难免会遭受到浮冰的碰撞作用，在极端冰碰载荷作用下船体结构难免会出现塑性变形，甚至会出现破坏断裂，如图1和图2所示，这将会导致人员伤亡、货物泄露以及环境污染等问题，甚至会造成沉船事故。例如，2015年，俄罗斯大型拖网渔船“远东”号与浮冰相撞，导致船舶沉没，56人死亡以及13人失踪；2019年1月，“雪龙”号破冰船在南极阿蒙森海与冰山相撞，船艏桅杆被撞倒，部分舷墙受损，对考察任务造成了严重影响。因此，亟需开展冰碰载荷下船舶结构动力响应评估方法研究，以探讨冰体碰撞下结构的塑性变形损伤。

船舶与海洋工程结构与浮冰碰撞是结构在短时间内遭受巨大冲击载荷作用下的一种复杂非线性动态响应过程，存在着材料非线性、几何非线性、接触非线性和运动非线性等问题。针对船-冰碰撞结构动力响应问题，目前尚未有较为成熟全面的理论方法，尤其是在冰材料力学模型方面，目前国内外学者主要是采用经验公式法、试验方法和数值方法来研究冰体碰撞下结构动力响应问题。随着计算机技术的不断提高，计算力学方法在不断发展，数值模拟方法已成为一种最为经济高效的方法[1]，主要包括有限元方法（FEM）、离散元方法（DEM）[2]、光滑粒子动力学方法（SPH）[3]以及近场动力学方法（PD）[4]等。有限元法在模拟介质间的间断和离散性方面能力有限，相比之下，离散元方法在模拟冰体离散特性方面具有一定的优势[5]。除此之外，有限元方法在模拟冰体断裂以及裂纹扩展方面有一定的局限性，而光滑粒子动力学方法和近场动力学方法的粒子离散形式使其在海冰破碎以及裂纹扩展方面具有更突出的优势[4]。然而，在船舶结构动力响应研究方面，尤其是针对结构塑性大变形以及断裂失效问题，有限元方法比离散元、光滑粒子动力学以及近场动力学方法显得更有优势，更加成熟。

由于冰体是一种力学性质较为复杂的天然复合材料，其材料力学特性的影响因素众多，尤其是对于海冰材料，主要受盐分、杂质等成分含量、温度、加载速率和形成年份等影响[6]。船-冰碰撞过程是一种瞬态动力响应问题，由于冰材料的固有特性，冰随着应变率的增加会由韧性破坏到脆性破坏[7-8]。由于冰体材料的复杂力学特性，不同条件下冰的力学性质差异较大，其材料数值模型仍然处于研究探索阶段。目前，国内外学者提出冰体理想化材料模型，主要为各向同性弹性失效模型[9-14]、各向同性弹塑性失效应变率相关材料模型[15-18]、可压碎泡沫型材料模型[19-25]、损伤材料模型[26-27]以及一些其它自定义弹塑性材料数值模型[28-35]等。在船-冰碰撞过程中，冰体的破坏包括了挤压、弯曲破坏、堆积等动态过程，其破坏力学特性不仅与应变率、温度等因素有关，还与其自身形成条件有密切关系，如年份、盐分、杂质等成分含量，利用数值模型方法来研究冰体破坏失效过程具有较大的挑战性。除此之外，通过一些冰体压缩、弯曲等试验测得冰体材料参数，如杨氏模量、屈服应力、剪切模量和泊松比等，这些结果具有较大的离散性，开发一种涵盖所有条件下的冰体材料模型是不现实的。因此，目前尚未有良好的数值模型来描述所有条件下的冰体材料力学特性，上述所提出的材料模型都有各自的局限性。

因此，本文主要总结了现有几种常见的冰体有限元数值材料模型，并与现有的一些试验结果、ISO 推荐的压力-面积曲线以及单位体积能量吸收经验值进行对比分析，讨论现有一些冰体材料模型的适用性和优缺点。

1 冰材料有限元数值模型研究进展

有限元数值计算方法是船-冰碰撞下结构动力响应的重要评估手段，但是碰撞过程中冰材料的数值模拟是其最为关键的难点问题，下面将介绍几种常见的冰材料数值模型。

在大多数船-冰碰撞情况过程中，船体与冰在接触过程中高应变率使冰体较明显地体现出脆性性质，故国内外许多学者将冰体材料模型简化为各向同性弹性失效材料模型。最早在针对冰雹与复合材料结构碰撞响应进行研究时，Kim 等[9-10]以及Anghileri 等[11]利用ANSYS/LS-DYNA 软件中的各向同性弹性断裂模型来模拟在碰撞过程中冰体的力学特性。该模型采用Von-mises 失效准则，将最大塑性应变定义为其破坏模式，把恒定最小压力模式定义为其分离模式，来模拟冰体失效行为。虽然冰体材料模型比较简单，可以节省计算时间，但是仅仅依据最大塑性应变和截断压力两个恒定值来判断冰体失效模式有所不足。结合相关冰体力学参数，一些学者[12-14]逐渐将此模型运用于船-冰碰撞分析之中。考虑到冰体材料力学特性受应变率影响较大，Carney 等[15]采用一种各向同性弹塑性失效应变率相关材料模型来模拟冰体失效，此模型广泛运用于冰雹的高速冲击数值模拟[16-17]。胡志宽[18]将各向同性弹塑性失效应变率相关材料模型运用于冰与螺旋桨的碰撞数值模拟研究中。此种材料模型虽然考虑到了应变率变化对冰体力学特性的影响，但是在船-冰碰撞结构动态响应研究中，碰撞速度有限，其应变率变化范围并不明显，没有冰雹高速冲击的应变率大，故该材料模型的适用性具有一定的局限性。

由于冰体材料参数试验测量数据离散性较大，具体的数值大小难以确定，一种可压缩泡沫型材料模型[19-21]被广泛运用于船-冰碰撞研究之中，通过运用自定义应力-体积应变关系来定义可压碎泡沫模型的力学行为。考虑到在碰撞过程中冰体与结构物接触区域存在着高压区，在其接触区域周围存在着低压区，如图3所示，但是上述的可压缩泡沫型材料模型还不能体现出这一区别。为此，Gagnon[22]对此材料模型所定义应力-体积应变关系进行了改进，通过对冰体模型进行分层，在每层赋予不同材料特性的锯齿形的应力-体积应变关系曲线，如图4所示，来对冰体的高压区域和低压区域进行模拟。类似于Gagnon 的分层方法，Kim 等[23-25]通过对冰体压缩试验现象与结果的观察和分析，定义了两条锯齿形的应力-体积应变关系曲线，并赋予冰体不同区域处，以表达冰体与结构物接触时高压区域和低压区域。此外，Kim 还运用“最大主应力”失效准则用来定义冰体破坏，即当冰体网格单元大于最大失效应力时，单元自动失效删除，采用此种材料模型的数值计算结果与试验结果吻合较好。相比于其它冰材料模型，可压碎泡沫材料模型在描述冰体力学特性时在一定程度上有较高的准确性，它可以宏观地体现出冰体与结构物相互接触时冰体的变化，可以比较真实地反映实验过程中的冰体与结构物的接触力大小，但是其缺乏强有力的物理解释，无法模拟出冰体脆性失效特性以及冰体在碰撞过程中出现的内部缺陷，该模型只反映特定条件下的冰体材料力学行为。

冰体在碰撞过程中主要体现出脆性失效性质，考虑到岩土、混凝土等材料本构模型与冰体材料在某种程度上有一定的相似性，而混凝土材料模型研究相对于冰体材料数值模型研究较为成熟，一些学者逐渐将一些成熟的混凝土损伤材料模型用于冰体材料模型数值模拟当中。Polach 和Ehlers[26]使用ANSYS/LS-DYNA 软件中损伤材料模型（*MAT_DAMAGE_3）来模拟冰体材料，使用Lemaitre 损伤模型来模拟冰体失效，当冰体损伤值大于临界损伤值时，冰材料就会失效。该模型的数值计算结果与实验结果吻合较好[37]。此外，Cai 等[27]基于混凝土材料本构模型，结合冰体材料力学模型，建立了表述应力-应变关系的本构方程的屈服面公式，以及失效应变与截断压力之间关系等式，提出了基于混泥土材料本构模型的冰体数值材料模型，数值计算结果与模型试验结果吻合较好。由于冰体材料特性的影响因素众多，如温度、应变率和盐度，开发出一种涵盖所有条件下的冰体材料模型是不现实的。相对于海冰来说，冰山冰为淡水冰，材料属性较为简单，可将其视作各向同性材料[38]。Liu 等[28-29]通过开发LS-DYNA 软件中的自定义各向同性的弹塑性材料模型来模拟冰山的力学行为，在所开发的弹塑性模型中，使用屈服面公式和流动法则来表征冰材料的应力-应变关系，建立失效应变与压力之间的关系曲线来定义冰体的失效准则。基于Liu 等人的自定义冰体弹塑性材料模型，Gao 等[31]提出了一种新的失效准则来模拟冰体失效，Song 等[33]利用Liu 的冰材料模型开展了考虑附连水质量的船-冰碰撞有限元计算研究。此外，Shi 等[34]提出了一种与温度梯度相关的冰体弹塑性材料模型，用于研究温度对船舶-冰山碰撞冲击力的影响。上述的自定义材料模型虽然在冰体失效准则和冰体三向应力状态模拟上更加准确，但是依然很难反映冰体真实裂纹延伸和断裂形态，相比于可压缩材料模型，无法准确模拟计算出冰体在碰撞过程中接触面压力大小分布特征。除此之外，Ince等[35]基于KOSORI冰模型的本构关系[36]，考虑应变速率、温度和盐度等参数的影响，对Johnson-Cook 材料模型进行了修正，通过ABAQUS软件用户自定义材料模型（UMAT）来实现冰材料模型的数值模拟。

在自然界中，无论是海冰还是淡水冰，不同环境下冰体的形成机理各不相同[6,36,38-39]，如在极地地区有1年冰和多年冰，不同年份的冰层的材料力学特性相差较大。冰的应力-应变特性依赖于应变率和温度，与冰的形成成分如盐分等也有直接关系，不同条件下测得的冰体材料参数也各不相同，如冰体密度、屈服应力、杨氏模量、剪切强度、拉伸强度、压缩强度以及断裂韧度等参数。开发一种涵盖所有条件下的冰体材料模型是难以实现的，每种冰体数值材料模型都具有其局限性。

2 模型试验与数值模拟结果对比分析

船-冰相互作用试验研究方法主要包括压缩和拉伸试验、弯曲试验、贯穿与剪切实验、压痕试验以及碰撞试验[5]。其中，压缩与拉伸试验、弯曲试验、贯穿与剪切实验以及压痕试验主要用于冰体力学特性研究。在船-冰相互作用过程中，有些工况是低速大质量碰撞情形，通常采用准静态加载实验方法来研究结构的弹塑性动力响应；对于碰撞速度较大时，便采用动态冲击实验方法来研究结构的弹塑性动力响应。下面针对三种模型试验工况来进行数值模拟分析，并对不同冰材料模型的计算结果进行分析。值得说明的是，由于所分析的模型实验工况有各自的特点，所以对下面所选取的冰材料模型参数与所参考的原文中的冰材料模型进行了微小的改动，以提高数值计算结果与模型实验结果的可比性。

2.1 冰体压缩实验数值模拟

Kim 等[24]开展了不同加载速率圆锥冰体压缩实验，实验装置如图5 所示，还采用可压缩泡沫型冰体材料模型来进行相应的数值仿真研究。在数值模拟中通过定义两条锯齿形的应力-体积应变关系曲线，并赋予冰体不同区域处，采用“最大主应力”失效准则来定义冰体破坏失效准则，数值结果与实验结果吻合较好。将Kim 等人所开展的冰体压缩实验其中一组实验工况来作为有限元数值模拟对象，其主要实验参数如表1所示，有限元数值模型如图6所示。图7对比了模型实验和6组数值模型计算结果，可以发现：各向同性弹性失效模型（MAT_13）和混凝土损伤模型（MAT_78）的计算结果会出现多个峰值，并出现多次加卸载过程。这是由于这两种材料模型在加载过程中冰模型会出现网格失效和删除，导致压块与冰体存在多次接触的过程。因此，在模拟冰体压缩实验过程中，采用网格单元失效删除的方式来模拟冰体失效是不准确的，很难对冰体压缩力-位移变化机制进行预测，会导致最终所预测的载荷值偏小。

表1 冰体压缩实验和数值模拟主要参数Tab.1 Main parameters of ice compression test and numerical simulation

图7 圆锥冰压缩实验与数值的力-位移曲线结果对比Fig.7 Comparison of experimental and numerical force-displacement curves for conical ice compression

2.2 冰体水平碰撞实验数值模拟

Cai 等[27]利用水平冲击试验机开展了楔形冰与船体板模型碰撞实验，实验装置如图8 所示。本文选取的碰撞实验工况主要参数如表2 所示，对应的有限元数值模型如图9 所示，冰体前端长度为200 mm，前端夹角为90°。下面选取了5种材料模型来进行对比，模型实验与有限元数值计算结果对比如图10所示。图10对比了各向同性弹性失效模型（MAT_13）、可压碎泡沫型材料模型（MAT_63）以及损伤模型（|MAT_78，MAT_153），从对比结果可以发现：采用低压区的低应力-体积应变曲线（M2）的可压碎泡沫型材料模型的计算结果与实验结果相比误差比较大；在碰撞力预测方面，采用各向同性弹性失效模型和损伤材料模型（MAT_78）的计算结果与实验结果比较接近；在船体板结构变形预测方面，采用损伤材料中的MAT_78 和MAT_153 单元材料模型的计算结果与实验结果比较接近。在船-冰碰撞数值模拟过程中，由于可压碎泡沫型材料模型很难模拟出冰体破坏失效特性，在模拟冰体破坏几何特征方面存在很大的不足之处。相比较而言，船-冰碰撞过程中冰体与船体局部结构相互作用特性研究方面，由于冰体在碰撞挤压过程中会出现以颗粒形式剥落的现象，采用如各向同性弹性失效模型（MAT_13）和损伤材料模型（MAT_78、MAT_153）等具有单元失效删除功能的冰材料模型来进行数值模拟具有较大的优势。

表2 冰体水平碰撞实验主要参数Tab.2 Main parameters of ice horizontal impact test

图10 楔形冰碰撞实验和数值的力-位移曲线结果对比[27]Fig.10 Comparison of experimental and numerical force-displacement curves for wedge ice collision

2.3 冰体落体碰撞实验数值模拟

Ince等[35]开展了圆锥冰体落体碰撞实验，实验装置如图11所示，此外还进行相对应的有限元数值模拟，基于KOSORI冰模型本构关系，采用了改进的Johnson-cook材料模型来模拟冰体力学行为，数值计算结果与实验结果吻合较好。本文选取的碰撞实验工况主要参数如表3 所示，对应的有限元数值模型如图12 所示。从图13 的对比结果可以发现：不同的冰材料模型数值计算结果差距较大，在板的变形为60 mm之前时，MAT_63（M2）材料单元模型、MAT_78单元模型与KOSORI冰模型的计算结果比较接近，与实验结果吻合较好：当板的变形超过60 mm时，上述的冰材料模型计算结果差异较大。

表3 冰体落体碰撞实验主要参数Tab.3 Main parameters of ice drop impact test

图13 圆锥冰落体碰撞实验和数值模拟的力-位移曲线结果对比Fig.13 Comparison of the experimental and numerical force-displacement curves for conical ice drop impact

3 冰材料数值模型可行性评估方法

在船-冰碰撞下结构动力响应数值模拟中，冰材料数值模型的准确性与可行性评估是船-冰碰撞数值研究中的一个重要环节。目前国内外学者主要采用的评估方法有模型实验验证方法[20,24,26-27,35]、ISO 规范等推荐的压力面积曲线（P-A）对比方法[13,27-28,31,40]以及比吸能（SEA，Specific Energy Absorption）判定方法[41-43]。模型实验方法是评估数值模拟方法可行性的最好方法之一，可以通过实验直接观察和分析冰体失效过程，并获得直观的实验数据，但是模型实验成本较高，耗时较长，且其提供的验证数据有实验条件限制，具有较大的不确定性。为此，许多学者基于大量实验测量得到的数据，对冰体碰撞过程中接触压力和接触面积进行分析，通过拟合得到压力-面积关系，从而为数值模型验证提供参考标准。比如，Masterson 等[44]总结了一系列的试验数据，提出了一个由许多研究人员根据现场测量数据拟合得到的P-A关系式，并将其纳入ISO 规范[45]。Cai等[27]开展了冰体冲击刚性墙的有限元数值计算，参考一些现有的有限元计算结果，总结了不同冰体材料模型下压力-面积值，并与ISO/CD 19906[43]以及API/CSA[46]推荐的压力-面积曲线进行对比，如图14 所示。利用现有的一些基于实验数据得到的P-A曲线，来对所提出的冰材料数值计算模型进行标定和校对，是目前比较常用的方法。

图14 不同冰体材料模型P-A值与ISO规范推荐曲线对比[27]Fig.14 Comparison between P-A values calculated by different ice material models and P-A curve recommended by ISO rule[27]

对于一种破坏失效吸能材料，单位的质量材料破坏失效所吸收的能量值（SEA）常用来衡量材料的吸能特性。船-冰碰撞过程中冰体体现出一种较强的挤压破坏失效特性，因此，SEA也用来表征冰体材料的吸能特性。图15展示了不同缩尺比实验测定得到的SEA 值，根据冰体破坏体积Vc，主要划分为三个区域：区域1，当1.0E-8 m3≤Vc≤1.0E-7 m3时，SEA值的范围大致在1.0E4 J/kg～1.0E5J/kg之间；区域2，当1.0E-7 m3≤Vc≤1.0E-5m3时，SEA值的范围大致在1.0E3 J/kg～1.0E4J/kg之间；区域3，当1.0E-5m3≤Vc≤1.0E-3m3时，SEA值主要集中在1.0E4 J/kg 左右。在各自区域内，随着冰体破碎体积的增大，SEA 值在逐渐减小。由于图15 中总结的实验数据所对应的模型实验冰体尺度不同，所以展现出了SEA 值的三个划分区域[42]。根据Cai等[27]文中的楔形冰与板模型碰撞实验工况，主要碰撞参数如表2 工况1 所示，开展不同冰材料模型下有限元数值计算，得到不同破碎体积下的SEA 值，如图16所示。由于受到模型实验尺寸的限制，随着碰撞时间的变化，冰体的破坏体积范围主要集中在区域3。图16 给出了不同冰材料模型下不同冰体破坏体积所对应的SEA数值计算值，可以发现：图中的冰材料模型所计算得到的SEA值都小于实验参考平均值10 000 J/kg，各向同性弹性失效模型（MAT_13）[13]和混凝土损伤模型（MAT_78）[27]计算得到的SEA 值随着破坏体积的增大有下降的趋势，而模拟高压区的可压碎泡沫型材料模型（MAT_63，M1）[22]计算得到的SEA 值随着破坏体积的增大有着明显的上升趋势，对于模拟低压区的可压碎泡沫型材料 模 型（MAT_63，M2）[22]和 损 伤 模 型（MAT_153）[26]计算得到的SEA 值没有很明显的变化趋势，相对比较稳定。除此之外，从表4 可以看出，各向同性弹性失效模型（MAT_13）和损伤模型（MAT_153）所计算得到的SEA 平均值相对于其它模型的SEA 平均值和实验参考值较小。每种材料模型所定义的材料参数都有着各自对应的适用条件，与图15 中的实验所对应的条件有所区别，因此模型实验的SEA值和数值模拟的SEA值有所差别。

图15 不同缩尺比实验测定的SEA值[42]Fig.15 SEA values measured by model tests with different scaled ratios[42]

图16 不同冰材料模型下不同冰体破坏体积所对应的SEA数值计算值Fig.16 Numerical calculation values of SEA under different ice material models versus different ice crushed volumes

表4 不同冰材料模型下的SEA数值计算值Tab.4 Numerical calculation values of SEA under different ice material models(Experimental reference value:10 000 J/kg)

4 结 语

在自然界中，无论是海冰还是淡水冰，不同环境下冰体的形成机理各不相同，如在极地地区有1年冰和多年冰，不同年份的冰层的材料力学特性相差较大，此外，冰体一些主要的力学参数和力学性质还受应变率、温度以及盐度、杂质等冰的形成成分因素的影响。因此，对冰船作用有限元数值模拟方法来说，开发一种涵盖所有条件下的冰体材料数值模型是目前冰工程领域里一项极具挑战性的研究工作，其中建立准确的冰材料本构关系是其关键技术和难点问题，仍处于研究探索阶段；除此之外，在船-冰碰撞数值模拟过程中，冰体变形破坏失效模式和裂纹扩展的真实模拟也是今后亟需解决的有限元数值模拟难点问题。

本文总结了一些现有冰材料有限元数值计算模型，分析了各自的特点，并选取了几种模型实验工况，开展了有限元仿真计算，对比并分析了不同冰材料模型计算结果的差异性；此外，还总结了一些现有的冰材料数值模型计算可靠性验证方法，针对一些算例，讨论并分析了用于冰材料模型验证的P-A曲线验证方法和比吸能（SEA）判定方法。从实验结果与数值模型计算结果对比可知，每种材料模型都有各自的适用范围及前提条件，目前所开发的冰材料数值模型很难适用于各种实验工况，利用模型实验验证方法来对冰材料模型准确性和可行性进行判断具有一定的局限性。为此，考虑温度、应变率以及冰体形成条件等影响因素，建立不同实验工况下的实验结果数据库，如SEA 值，为冰材料数值模型提供验证数据，是目前船-冰碰撞数值方法研究的关键之处。