船-冰碰撞中考虑温度影响的冰体材料本构模型研究

2023-03-01俞同强刘俊杰王自力

船舶力学 2023年2期

俞同强，刘昆，刘俊杰，王自力

（1.江苏科技大学船舶与海洋工程学院，江苏镇江 212003；2.中国船舶科学研究中心，江苏无锡 214082）

0 引言

随着全球气候变暖，极地地区的冰雪消融日益加剧，北极海冰尤其是夏季海冰正以每十年10%的速度减少[1]，据估计，北冰洋十年内将出现夏季无冰年，北极航道的全面开通成为可能。对于我国来说，利用北极东北航道，从上海以北港口到欧洲西部港口的航程比传统航程缩短25%～55%[2]，不仅可以大大减轻贸易成本，而且可以摆脱对马六甲海峡的依赖，同时对于我国渔业、旅游业以及油气开采等行业具有深远影响，其重要性不言而喻。

北极航道的通航对极地船舶的安全性提出了更高的要求，船-冰碰撞所带来的结构强度与冰载荷问题日益受到人们的重视。由于极地海域复杂的环境特征，即使在夏天也有浮冰存在，可能导致船舶破损、进水甚至沉没，造成重大的生命财产损失。即使辅以现代设备，船-冰碰撞问题仍然难以避免。2019年“雪龙”号在执行我国第35次南极考察任务期间，在阿蒙森海密集冰区航行中，因受浓雾影响，在南纬69°59.9'、西经94°04.2'位置与冰山碰撞，船艏桅杆及部分舷墙受损，所幸无人员受伤。可见，开展极地船舶船-冰碰撞性能研究，对于确定冰载荷大小、评估船体结构性能、减少船-冰碰撞的事故损失以及对未来极地船舶的结构设计等均具有非常重要的意义。

研究船-冰碰撞问题首先需要确定冰载荷的大小，目前主要有理论计算、试验以及仿真三种方法。理论研究方法主要是在一定尺度下，通过能量法等力学理论结合单轴、压痕试验数据，得到不同工况下冰载荷的经验计算公式[3-5]，方法简单但是精度较低，适用范围不够广泛。随着计算机技术快速发展，数值算法不断进步，商业有限元软件被越来越多地用于解决海冰与结构作用时冰材料的模拟和冰载荷的计算问题，此时船-冰碰撞问题的研究已发展到一个新的阶段。但是，由于海冰材料复杂的物理和力学特性，目前对海冰以及冰载荷的研究仍然面临诸多困难，其中对其本构关系的合理准确描述一直是影响相关问题计算分析的关键因素。国内外学者相继提出粘塑性、粘弹塑性、弹塑性、弹脆性和泡沫材料等多种材料模型[6-10]，具有一定的实用性但也有各自的局限性。目前由于商业有限元软件中缺乏海冰材料相关定义，使用自定义材料本构来模拟海冰力学特性和计算冰载荷成为一种切实、有效的手段[11-12]。

本文建立一种带脆性失效特征的弹塑性冰材料本构模型，考虑温度对屈服函数及失效准则的影响，编写相应材料子程序并通过二次开发嵌入到LS-DYNA 材料库，以此为基础开展船-冰碰撞研究，分析船体与海冰相互作用时冰体和结构的损伤和变形特点。本文研究成果可为极地船舶碰撞损伤评估和结构设计制造提供参考。

1 考虑温度影响的海冰材料本构模型

1.1 考虑温度影响的海冰材料本构理论

自然界中的海冰是由结晶冰、卤水、空气和杂质等组成的复杂混合物，其成分比例随外界条件、时间和空间的变化而变化，其中最重要的控制因素是温度，所以从材料学角度来看，海冰可以视为温度敏感性复合材料[13]。随着海冰力学试验尤其是多轴试验的广泛实施，海冰三维应力状态和材料力学特点得到全面直观反映[14-17]，为理论模型建立和数值仿真提供了有效的数据支撑。其中，Derradji-Aouat[18]在总结归纳三轴试验结果的基础上给出了各向同性淡水冰、冰山冰以及柱状冰的多曲面椭圆屈服方程，该屈服方程得到了广泛接受和使用。

即

式（3）可以写为

该屈服方程即为Tsai-Wu 屈服准则，其中，f( )p,J2为屈服函数，J2为第二偏应力不变量，a0、a1、a2为系数。Tsai-Wu 准则能够体现出复杂应力状态下静水压力能影响材料屈服的性质，同时十分适用于描述对于拉压强度相差很大的海冰材料。该屈服方程在主应力空间下的屈服曲面及在坐标轴上的投影如图1所示，可以看到在三维空间下，屈服面呈椭球型，在坐标面上的投影为椭圆形，且由于未考虑各向异性，各坐标面下屈服面投影形状保持一致，对于同一类型冰，屈服面在温度等因素影响下膨胀或收缩，同时保持椭球焦心不变，即λ、pc、η的数值保持不变。

图1 主应力空间下的屈服曲面及其在坐标轴上的投影Fig.1 Yield surface in principal stress space and its projection on coordinate axis

由以上分析可知，参数a0、a1、a2的值仅与参数qmax有关，在三轴试验中，qmax=σ1-σ3，因此可以较为系统地总结出qmax与试验中各变量的关系，Derradji-Aoua 给出式（5）～（6）来表征qmax与应变率和温度的关系。

船-冰碰撞下冰体通常处于高应变率（10-3以上）范围，在此种场景下海冰表现为弹脆性，应变率的影响可以不计。根据Gagnon 等[15-17]的三轴试验数据，针对该应变率附近不同温度下的试验数据进行统计、拟合，得到的屈服曲线如图2所示。

因此，参数a0、a1、a2可以表示为温度的函数ai(T)，将该屈服方程发展为受温度T影响的多重屈服面准则作为本文海冰本构的屈服准则：

图2 不同温度下屈服方程曲线Fig.2 Yield equation curves at different temperatures

该本构模型在弹性阶段，满足各向同性广义胡克定律，其增量形式为

式中，剪切弹性模量G=E/2( 1 +μ， )E为弹性模量，μ为泊松比，γj为工程切应变，γj= 2εj(j = xy,yz,zx )，Δε0为静水应变增量，Δσ0为静水应力增量。在塑性阶段，由于屈服面的外凸性，塑性应变增量为

式中，dγ为塑性一致性参数增量。采用关联流动法则，总应变为弹性应变与塑性应变之和，

式中：dσij为应力增量；δij为Kronecker符号，其值为0（i≠j）或1（i=j）。

此外，还必须建立合理的单元失效准则来模拟海冰材料的脆性破坏特征，海冰存在多种破坏形式，在受压情况下剪切和挤压破坏是两种主要形式。但在实际碰撞问题的处理中，当采用剪切破坏准则时，实际上过高估计了材料的抗剪强度，往往会导致单元提前失效。这是因为冰在经历塑性阶段的挤压过程中，单元静水应力增加，此时单元的强度明显增加，此种失效标准显然更难达到。基于以上分析，本文采用以率无关和累积塑性应变为表征的动态失效准则来定义单元的失效，该准则起源于Jordaan的经验失效方程[19]，如下式所示：

式中，a、p1为常数，p为静水应力。

本文进一步发展该经验失效准则，具体为：考虑温度的影响，将p1设置为屈服方程较大的一个根，同时，给定初始失效应变ε0= 0.01，当＞εf时材料失效，其中为等效塑性应变=，εf为失效应变，εf=ε0+(p/p1- 0.5)2。此外，考虑到拉压强度的不同，给定拉伸强度极限，当p＜pc时，单元删除。pc为截断应力，设为拉伸极限2 MPa。该失效准则特点是注重由静水应力带来的高应力区和低应力区的承载能力差别，同时抑制摩擦滑移，在碰撞问题中弱化剪切破坏的影响。

1.2 子程序算法实现

积分率本构方程的数值算法被称为应力更新算法或本构积分算法，对于率无关塑性，常用的本构积分算法有显式算法、隐式算法以及半隐式算法。一个广义的率无关的弹塑性本构模型[20]可以写为

式中：εij分别为总应变、弹性应变以及塑性应变；q为内变量；hα为塑性模量；rij为屈服函数的偏导数为塑性应变增量；γ为塑性参数，由加卸载准则确定。本文采用最近点投影法来求解，其为一种半隐式的返回映射算法。计算中对塑性参数采用隐式而对塑性流动方向和塑性模量采用显式算法，在步骤结束时计算塑性参数的增量并同时强化屈服条件。在某时刻t的下一个时间步长Δt中，试验应力由胡克定律求得，若在弹性阶段内，则直接更新应力，若超过屈服极限，则进入迭代步，其各参数计算公式为

式中，Δγn+1即为塑性参数γ在n+1 步的增量，C为弹性系数矩阵。在每一个积分步内，首先计算方程塑性参数值与应力分量，

每个迭代步内检查收敛条件，若收敛，退出迭代算法。由此可得n+1 步的塑性应变增量为

图3 本构算法流程Fig.3 Flow chart of constitutive algorithm

1.3 单元测试和材料验证

建立单元测试模型分别施加恒定压缩速度，测试温度为-10 ℃，考察单元变形中相关参数的变化情况，在子程序计算中，输出第二偏应力不变量J2、静水应力p以及失效应变的值，得到J2-p、εf-p的曲线分别如图4～5 所示。可以看到，偏应力不变量J2与静水应力P呈二次曲线关系，其中标准曲线为理论曲线。图4 中在a～b范围内，输出值低于标准值，此时单元在弹性范围内，即屈服应力f(p,J2,T)＜0。当到达b点之后，两线重合，单元达到屈服状态，输出值与理论值一致，初始屈服点在10 MPa 左右，与理论解符合。失效应变与静水应力也呈二次曲线关系，静水应力输出值与理论值一致，当静水应力为55 MPa左右时，失效应变达到最小值为0.01，与理论解相符，且能合理表现冰材料的脆性特征，子程序算法的执行过程是正确的。

图4 J2-p曲线图Fig.4 Relationship curve of J2-p

图5 εf-p曲线图Fig.5 Relationship curve of εf-p

同时，建立刚性板-球形冰体挤压模型，通过压力-面积曲线验证材料的准确性。球体直径为1 m，施加1 m/s 的恒定速度，0.5 s 时球体的损伤变形如图6 所示。可以看到球形冰体表面发生破坏失效，失效面不平整，接触边缘有冰体脱落，冰体接触面内应力分布不均，有明显的高应力区域和低应力区域，最大应力为23.5 MPa左右，表明在三维应力下，材料的强度得到了一定程度的增强。

图6 0.5 s时球形冰体的破坏情况Fig.6 Damage of ice blocks at 0.5 s

图7给出了碰撞力-时间曲线。可以看到，在整个作用过程中，碰撞力呈现典型的非线性特征，刚开始接触时，由于球体接触区域面积变化较快，碰撞力上升明显且上升速度较快，其后，面积变化速率逐渐减小，碰撞力上升趋势变缓，总体趋于平稳。计算得到的压力-面积曲线如图8 所示，可以看到，压力-面积曲线总体趋势与ISO 标准[21]以及其他学者实际测得的变化趋势相同，初始阶段由于接触面积较小，其值小于标准值，且此时误差较大，随着面积增大，平均压力值逐渐减小，与标准值逐渐接近，最后阶段，压力值趋于平缓，且与标准曲线的一致性较好，实际情况中，有围压下海冰强度最大约为30 MPa 左右，因此在使用中应以面积较大（≥0.25 m2）时为准。由于ISO 标准的压力-面积曲线是基于极限状态设计（ULS）以及偶然极限状态（ALS）方法，压力测量取极限值，有一定的富裕度。本文给出Timco 和Molikpaq[22-23]通过试验测得的压力-面积曲线（试验温度-5 ℃～10 ℃），相比之下，本文所得压力-面积曲线在面积较小时与试验得到的曲线吻合得更好。因此，该材料模型能够反映海冰强度的真实变化情况。

图7 碰撞力-时程曲线Fig.7 Force-time curve

图8 压力面积曲线Fig.8 Pressure-area curve

图9 和图10 给出了碰撞过程中刚性板局部（1 m×1 m）对角线节点的位置在0.2 s、0.4 s、0.6 s、0.8 s四个时间点各个节点的节点力。可以看到，随着作用过程的进行，存在着高压区与低压区的相互转化，越接近中心的节点力变化幅度较大，沿着四周的逐渐减小。

图9 节点位置图Fig.9 Nodal position

图10 节点力Fig.10 Nodal force

考虑温度的影响，约束单元底端的四个节点，顶端施加恒定的压缩速度，直至单元失效，输出不同温度下失效时单元的等效应力和等效应变。不同温度下单元失效时的等效应力和应变以及计算所得压力-面积曲线如图11～13所示。可以看到，随着温度的不断降低，单元失效时的应力增加，表明材料强度不断增加，材料逐渐“变硬”，同时失效应变逐渐减小，材料逐渐“变脆”，与实际情况相符。不同温度下的压力-面积曲线变化趋势相同，压力趋于平稳时，在较低温度下的压力值较高。

图11 单元失效时应力-温度曲线Fig.11 Failure stress-temperature curve in single unit test

图12 单元失效时应变-温度曲线图Fig.12 Failure stain-temperature curve in single unit test

通过以上验证，可以表明本文提出的自定义冰材料模型子程序执行正确且能够较好地模拟海冰材料特性和冰载荷，合理体现其温度敏感性特性，可运用于实际的碰撞场景中。

2 船-冰碰撞场景仿真

图13 不同温度的压力-面积曲线Fig.13 Pressure-area curve at different temperatures

2.1 场景介绍

通过调研分析确定典型计算工况，对于船艏部位的撞击，选取船舶排水区域内船艏典型结构进行研究。本文选取艏部结构中的球鼻艏结构与块状海冰的碰撞场景，如图14 所示。船体艏部选取球鼻艏至艏部舱壁向后6.5 m处，尺寸为20 m×48 m×26 m，建模采用4 积分点、5 节点壳单元，网格尺寸为200 mm，船体结构所用材料为DH32型低温钢，材料属性由低温准静态拉伸试验所得，参数如表1所示，其真实应力-应变曲线如图15 所示。块状冰体尺寸为10 m×10 m×5 m，采用实体单元，网格尺寸为100 mm。目前极地商船的航速一般不超过16 kn[24]，本文中碰撞相对速度为8 m/s，温度为-10 ℃。

图14 碰撞场景Fig.14 Collision scenario

表1 DH32钢的力学参数Tab.1 DH32 steel mechanical parameters

图15 真实应力-应变曲线Fig.15 Real curves of stress-strain

2.2 结果分析

冰体、结构的损伤变形如图16 所示。可以看到，碰撞过程中冰体在与船艏接触区域发生较大程度的破坏，艏部结构也发生较大程度变形。冰体最大应力达到15 MPa 左右，破碎表面不规则并伴有碎片产生，冰体有明显的高低应力区且高应力区分布具有显著的局部性。由于碰撞区域较小，速度较大，艏部结构最大应力达到764 MPa 左右，外板变形呈现连续凹陷的形式，内部构件交叉区域产生一定程度的应力集中现象，应力分布区域较为集中，可见球鼻艏内部加强结构在碰撞过程中起到重要作用。

图16 冰体与结构损伤变形Fig.16 Damage and deformation of the ship structure and ice block

图17 给出了碰撞过程中碰撞力时程曲线。可以看到，总体上碰撞力呈先迅速上升后逐渐趋缓的特征，在整个作用过程中碰撞力曲线表现出显著的非线性特征。在碰撞初始阶段由于接触面积迅速增大，碰撞力上升速度较快，随着碰撞的进行，冰体单元存在同时大量失效的情况，碰撞力有瞬间减小的情况，此后接触面积增长速度逐渐变缓，碰撞力增长速度减慢，在最后阶段，面积趋于最大值，碰撞力变化趋于平稳略有上升。