0 引　言

破冰船是北极航道开辟的重要载体，当破冰船与极地多年冰层接触破坏时，船体将不可避免地承受冰载荷的强冲击作用。该冲击力不仅给船载设备可靠性、船员舒适性造成危害，严重还将对船体结构安全性造成不可逆的破坏，加之北极天寒地冻，易造成船毁人亡的严重事故。因此，开展破冰船冰致冲击载荷研究意义重大。

在冰致冲击载荷研究方面，王键伟等^[1]基于对正应变测点布放位置敏感性的分析，合理制定光纤光栅传感器的安装方案，采用影响系数矩阵法对冰载荷进行了识别。Böhm等^[2]克服了冰载荷只有作用在普通仪表区域内才能准确获取的限制，进一步提出了一种确定负载位置的方法。Wang等^[3]基于典型极地船舷侧板架的模拟应变，确定了监测区域内识别误差的空间分布特征，提出了用最小二乘法（LSF）估计失效点处应变。Wu等^[4]对北极航行中船舶船体冰载荷进行统计分析，将最小二乘支持向量机（LS_SVM）算法应用于冰荷载识别，并通过实验验证了识别模型的可行性。曲悦含等^[5]采用数值模拟方法对某极地破冰船与大块浮冰碰撞产生的冰载荷展开了研究，利用ALE流固耦合方法实现船体-浮冰-海水的耦合作用，分析了浮冰尺寸、形状和船舶航速对船冰碰撞冰载荷的影响。何帅康等^[6]以内河破冰船为研究对象，开展了船体结构冰载荷的实船测量，发现不同厚度的河冰载荷作用形式和周期与海冰具有较高相似性。白雪寒^[7]考虑船-冰-水之间的耦合作用，提出了一种基于SPH（光滑粒子流体动力学方法）-FEM（有限单元法）算法的冰载荷数值预报方法，实现了考虑随机特性的冰载荷预报。徐张等^[8]基于Ls-dyna软件构建了船水冰全耦合模型，对极地运输船舶与规则层冰和不规则层冰发生碰撞分别开展了数值模拟，结果表明不规则层冰载荷下引起的船体结构响应更为剧烈，碰撞力均值也更大。吴刚等^[9]以“雪龙2”号为载体，在南极中山站完成了冰弯曲强度试验和艏艉双向破冰时的实船数据测量，摸索了一套极地冰区破冰船破冰试验方法。孔帅等^[10]基于CUDA-C并行处理技术，建立了可用于海冰模拟的高性能离散元方法（DEM），采用Voronoi分割算法生成平台作业时的碎冰域，研究发现极区浮式平台冰区运行时冰载荷受冰厚影响较大。庞福振等^[11]在松花江开展了内河破冰振动试验，测量了破冰及敞水工况对破冰船振动特性的影响，结果表明，船舶破冰作业时会引起船体总振动及局部振动响应，且船舶局部振动与破冰载荷作用位置及船体板架固有频率密切相关。

由以上分析可知，学者们针对破冰船破冰冲击载荷开展了大量且富有成效的研究，揭示了破冰速度、冰层厚度对冰载荷的影响机制，但实船受海冰作用下的冰载荷模拟仍有待进一步丰富，破冰冲击过程中海冰破坏模式尚有待进一步揭示。

1 破冰船破冰冲击载荷模拟方法

本文基于ALE（Arbitrary Lagrange-Euler）算法和侵蚀接触算法，建立船冰作用力学计算模型，并以“雪龙2号”破冰船为研究对象，数值模拟船-冰碰撞过程。

1.1 船冰碰撞有限元理论基础

船冰可以看作固-固相互碰撞，为模拟船舶破冰过程，本文采用动态接触分析算法，发生接触的2个表面被称为主面和从面，相应的接触算法包括动力约束法、分配参数法和对称罚函数法。本文采用对称罚函数法，罚函数原理示意图如图1所示。

以节点为例，如果当前时间步发生穿透，罚函数接触力f计算式为：

$ \left\{ \begin{gathered} {f} = {k}{\delta}，\\ {k} = p_f\frac{KA^2}{V}。\\ \end{gathered} \right. $

(1)

式中：k为接触刚度因子；δ为穿透量；p_f为接触刚度惩罚因子；K、A、V分别为接触单位体积弹性模量、接触段面积和主段体积。

船冰碰撞的时间离散采用显式时间积分方法，根据动力学问题的有限元方法，结构动力学方程经离散化变为：

$ {\boldsymbol{M}}\ddot {\boldsymbol x}(t) = {\boldsymbol P}(t) - {\boldsymbol F}(t) + {\boldsymbol H}(t) - {\boldsymbol{C}}\dot {\boldsymbol x}(t)。$

(2)

式中：$ {\boldsymbol{M}} $为质量矩阵；$ \ddot {\boldsymbol x}(t) $为节点加速度向量；$ {\boldsymbol P}(t) $为外力载荷向量；$ {\boldsymbol F}(t) $为内力向量；$ {\boldsymbol H}(t) $为沙漏阻力向量；$ {\boldsymbol{C}} $为阻尼矩阵；$ \dot {\boldsymbol x}(t) $为节点速度向量。在碰撞有限元计算中，当系统总沙漏能小于系统内能总和的5%，认为沙漏现象控制有效。

对离散化后的动力方程，采用显式中心差分方法求解，递推公式如下：

$ \left\{\begin{gathered}\ddot{\boldsymbol{x}}(t_n)=\boldsymbol{M}^{-1}\left[\boldsymbol{P}(t_n)-\boldsymbol{F}(t_n)+\boldsymbol{H}(t_n)-\boldsymbol{C}\dot{\boldsymbol{x}}(t_{n-1/2})\right]，\\ \boldsymbol{\dot{x}}(t_{n+1/2})=\boldsymbol{\dot{x}}(t_{n-1/2})+\boldsymbol{\ddot{x}}(t_n)(\Delta t_{n-1}+\Delta t_n)/2，\\ \boldsymbol{x}(t_{n+1})=\boldsymbol{x}(t_n)+\boldsymbol{\dot{\boldsymbol{x}}}(t_{n+1/2})\Delta t_n。\\ \end{gathered}\right. $

(3)

式中：${t_{n - 1/2}} = ({t_n} + {t_{n - 1}})/2$，${t_{n + 1/2}} = ({t_n} + {t_{n + 1}})/2$，$\Delta {t_{n - 1}} = {t_n} - {t_{n - 1}}$，$\Delta {t_n} = {t_{n + 1}} - {t_n}$；$ \boldsymbol{\ddot{\boldsymbol{x}}}(t_n) $为${t_n}$时刻节点加速度向量；$ \boldsymbol{\dot{x}}(t_{n+1/2}) $为${t_{n + 1/2}}$时刻节点速度向量；$ \boldsymbol{x}(t_{n+1}) $为${t_n}$时刻节点位置坐标向量。

对积分步长，一般情况下该参数主要取决于网格中最小单元尺寸，单元类型临界积分步长可表达为如下统一形式：

$ \Delta {t^e} = \alpha ({l^e}/c)。$

(4)

式中：$\Delta {t^e}$为临界时间步长；$\alpha $为时间步长因子；${l^e}$为特征尺寸；$c$为纵波速度。

1.2 海冰参数三点弯曲数值模拟

在海冰模型方面，首先采用有限元法对海冰参数进行三点弯曲反演模拟，算例模型如图2所示，试样为70 mm×70 mm×650 mm的标准三点弯曲冰梁。冰梁置于间距为600 mm的固定支座上，压头对准冰梁顶端中点，压头以0.30 mm/s的速度持续下压，直至冰梁发生弯曲破坏。

由图2(b)可知，当冰模型弯曲强度为0.7 MPa左右，弹塑性模型对海冰破坏模拟具有理想效果，此时可得如表1所示的海冰参数。

1.3 内河破冰冲击载荷试验验证

为说明ALE算法数值模拟的有效性，这里以内河破冰船为研究对象，以破冰船结构振动响应为评价指标，验证破冰冲击载荷数值模拟的有效性。内河破冰船总长29.7 m，排水量134 t，吃水1.15 m，本文采用实尺度模型进行破冰冲击载荷模拟，首先根据图纸资料建立内河破冰船几何模型，并将结构离散得到有限元模型，船体材料选用刚体材料，网格单元为壳单元，单元大小为0.2 m，单元总数为9519个。在船外侧设置水域及空气域，水域和空气域长度均为60 m，宽度为40 m，水域高度3.6 m，空气域高度2.6 m，水域和空气域采用六面体网格，单元总数为57400个；流场外边界均采用无反射边界条件，水和空气均采用Null材料模型，流体模型材料参数如表2所示。

冰模型长度为56 m，宽度35 m，厚度0.25 m，冰材料参数同表1。冰有限元模型采用六面体网格，大小为0.15 m，网格总数为175032，具体如图3所示。

为方便对比破冰冲击载荷的有效性，这里以破冰冲击下的结构振动为校核标准，破冰振动采用声固耦合方法，破冰船结构外侧设置无限大水域，网格尺度0.2～0.3 m，由近及远呈现逐渐增大的趋势，并在外流场外表面设置无反射边界条件，内河破冰船声固耦合模型如图4所示。考核点位置如表3所示。

试验破冰速度如图5所示，分析频率范围为0～100 Hz。

为验证方法的有效性，将计算得到的振动加速度响应与试验中记录的振动加速度响应进行对比分析，对比结果如图6所示。仿真结果与试验结果响应变化趋势吻合较好，最大总级偏差小于2.8 dB，验证了基于ALE算法进行冰载荷模拟的有效性。

2 典型破冰工况冲击载荷特性研究 2.1 破冰冲击模型

“雪龙2号”破冰船破冰冲击载荷模拟有限元模型如图7所示，船体所有节点在x方向以恒定航速运动，冰长约为垂线间长2倍，宽度约为船宽5倍，平整冰与水面初始间距为0.01 m。

破冰船主尺度参数如表4所示，船体材料选用刚体，弹性模量2.06×10¹¹ N/m²，泊松比0.3。船体采用2D网格划分，网格尺寸0.4 m；水域和空气域长宽分别为240 m、120 m，水域深度为16 m，空气域高度为8 m。流场外边界采用无反射边界条件，水和空气采用Null材料模型，有限元模型如图8所示。

2.2 连续式破冰模拟

本研究将数值计算中破冰船尾部进入冰场作为受力分析的初始时刻（T=0 s），如图9所示。图10对应给出了初始时刻冰层破坏模式。可以看出，由于破冰船持续破冰运动导致船首及舷侧受冰带区域对冰进行持续挤压，船首冰体发生变形后达到压缩强度临界值，冰体单元失效，从而发生破冰现象。与此同时，由于船首外飘型线，会导致首部斜侧方向冰层受斜向作用力，冰层在垂直于船体运动方向上的剪应力达到强度极限，冰体单元失效，冰层发生纵向剪切破坏并形成裂缝，从而发生冰层纵向剪切破坏。此外，由于冰船船首独特的外飘设计，破冰船运动过程中冰层受到正压力，使层冰发生弯曲，当应力达到冰材料极限后，冰层同样发生弯曲破坏，如图11所示。

图12为破冰船运动10 s后的破冰模式，此时破冰船在船首附近区域的破冰模式以剪切破坏为主，同时辅以挤压破坏等。图13为破冰船运动20 s时刻的破冰模式，可以看到，此时冰层破坏模式主要以弯曲破坏为主，在船首和舷侧发生碎冰脱落，并向水面以下运动，使破冰船水线以下区域也会受到船冰相互作用而产生碰撞力。总之，通过对不同时刻破冰模式分析可知，海冰与破冰船相互作用时主要有以下3种破坏模式：

1）挤压破坏。冰层被作用区域受挤压而逐渐破碎。

2）剪切破坏。冰层受到剪应力达到强度极限时产生的裂缝。

3）弯曲破坏。冰层沿斜面结构运动，随后受弯最终发生弯曲破坏。

图14为x方向破冰船破冰载荷合力时历曲线，由于破冰船在破冰过程中同时存在一种或多种破冰模式，导致船舶破冰载荷时历曲线成脉冲形式。如在T=6.2 s时刻载荷出现峰值，这是由于此刻冰层主要以弯曲破坏为主（见图15），此时破冰船会使冰层受到一个正压力，冰层发生弯曲，同时使破冰船受到一个较大支反力作用，从而导致载荷出现峰值。随着冰层局部失效脱落，破冰载荷出现降载现象，但随着破冰船持续运动，船首附近冰层区域会发生弯曲破坏，导致破冰船受到的冲击载荷又会增加，之后随着冰层失效再次出现卸载现象。

图16为T=9.4 s时刻和T=9.6 s时刻冰层失效局部放大图。可以看到，在T=9.4 s时刻，冰层破坏模式呈现剪切破坏、局部弯曲破坏及挤压破坏等多种复合模式，因此在T=9.4 s时刻载荷出现峰值。随着冰层的失效区域逐渐增多，大部分冰层出现破坏形成碎冰与船体脱离，造成载荷出现卸载现象，导致在T=9.6 s时刻载荷出现降低现象。

为了分析破冰船破冰冲击载荷特性，采用hann窗对船体表面节点受到的x方向冲击载荷进行傅里叶变换，图17为破冰船首部前侧附近冰致冲击载荷幅频曲线，其中x表示受力位置距艉封板的距离。可以看到，在船首附近区域，冰致冲击载荷主要集中在0～20 Hz频段内，船首破冰载荷以低频为主。这是由于破冰船低速运动时，船首通过缓慢压迫平整冰造成冰体呈现挤压破坏、弯曲破坏、剪切破坏等失效模式（见图17）。当然，由于平整冰失效模式出现的位置不同，会影响冰致冲击载荷的幅值，但主要在船首最先与冰作用区域的载荷幅值最大，如图17(a)所示。但当船冰作用出现局部弯曲时（见图18），也会对破冰船产生反向作用力，从而使得距离船首更远位置冰载荷幅值大于更近位置，图17(c)冰载荷幅值大于图17(b)。

2.3 冲撞式破冰模拟

为分析冲撞式破冰冲击载荷特性，设定数值模拟工况，破冰船初速度为8 kn，冰厚为3 m，破冰过程中速度逐渐减小为0。

图19为在航速8 kn、冰厚3 m工况下的冲撞破冰过程，本研究同样将船尾进入冰场作为受力分析初始时刻（T=0 s），终止时刻取破冰船运动15 s（T=15 s）。冲撞式破冰过程中，由于破冰船初动量较大，冰层破坏主要以弯曲破坏为主，剪切破坏及挤压破坏占比较小，如图19(a)所示。随着破冰船速度的逐渐减小，挤压破坏占比增大，如图19(b)所示。

破冰船冲撞破冰过程速度及x方向破冰载荷合力曲线如图20所示，破冰船冲撞破冰过程初始阶段，随着船体与冰层接触区域增大，船体x方向所受合力也逐渐增大（见图20(b)）；同时，随着破冰时间继续增加，破冰船速度逐渐减小，破冰船对冰层作用力也逐渐减小，导致船体x方向所受合力逐步减小（见图20(b)）。

3 结　语

1）完成了破冰船动态破冰数值模拟，并与已有内河破冰船破冰振动试验数据进行了对比，结果表明，破冰振动响应分析结果与实船测试结果总体吻合较好，最大总级偏差小于2.8 dB，验证了破冰冲击载荷ALE算法模拟的有效性。

2）对连续式破冰，破冰船持续破冰运动会对冰层产生挤压、剪切和弯曲破碎，剪切力和弯曲力主要来源于破冰船外飘型线设计，但不论何种破冰模式，冰致冲击载荷主要集中在0～20 Hz频段范围，船首破冰载荷以低频为主。

3）对冲撞式破冰同样存在层冰的挤压破坏、剪切破坏和弯曲破坏3种模式，破冰船初试运动将较大的冲量作用在平整冰，冰层破坏以弯曲破坏为主，剪切破坏及挤压破坏占比较小，随着破冰速度的减小，挤压破坏占比逐渐增大。