基于FEM/SPH耦合方法的极地运输船舶冰阻力预报研究

王志鹏，郝寨柳，田于逵，

（中国船舶科学研究中心，江苏无锡214082）

0 引 言

北极航道是联系亚欧美三大洲的最短航线，与通常的苏伊士运河和巴拿马运河航线相比，能够缩短数千公里航程，由此将大大降低商业运输成本，影响和意义深远。随着气候变暖，海冰融化加速，北极航道通航条件不断改善，未来对极地运输船舶的需求会明显增大。极地运输船舶不仅要满足开敞水域性能和经济指标，而且需具备一定的破冰能力以满足冰区航行需求。在极地运输船舶开发过程中，冰阻力预报是船舶破冰能力和动力需求分析的重要输入，也是总体设计的关键环节之一。在船舶冰阻力预报方面，主要有模型试验、经验公式、数值模拟等方法，其中，模型试验是最可靠的方式，基于实船测量和模型试验数据总结的经验公式也被广泛应用，目前常用的经验方法有Lindqvist 方法、Kei⁃nonen 方法、Riska 方法和Spencer 方法等[1-4]。相比于经验公式中模型过分简化，数值计算能够考虑海冰破坏模式和更多船型信息对冰阻力的影响。大量学者对海冰数值模型和破冰航行过程开展了研究，主要包括有限元、离散元、半经验数值模型等方法，使数值计算的精度和适用范围得到了持续改进和提高。Arne Gürtner 等[5]基于有限元方法，利用LS-DYNA 软件基于黏结单元建立海冰数值模型，对冰层与灯塔作用过程进行了数值模拟，获得了作用区域载荷分布和变化过程，并与实测结果进行了对比，验证了黏结单元法模拟海冰破坏的可行性；Wang 等[6]利用LS-DYNA 软件，通过用户自定义冰材料，将海冰破坏前表述为线弹性，采用多表面失效准则进行破坏判断，并将数值计算结果和破冰船实测结果进行了对比，验证了数值方法的可行性；任奕舟等[7]基于有限元理论，采用可破碎泡沫模型模拟冰材料，并通过冰锥试验结果验证其可行性，同时采用经验公式计算附连水质量并叠加在船体质量上，通过对破冰船不同航速下破冰过程进行数值模拟进而预报了船舶冰阻力；Su 等[8]采用离散元数值积分方法，耦合求解冰力和船体三自由度运动方程，建立船-冰作用时船体载荷数值模型，分析了船舶在破冰航行过程中冰阻力和船体局部载荷情况；季顺迎等[9]研究了不同形态海冰离散单元模型特点，对极区海冰动力特性、海冰重叠堆积、海冰-海洋结构物相互作用开展了数值模拟及分析；Tan 等[10]提出了基于半经验方法的数值模型，考虑压力-面积关系计算船舶与冰层间接触力，以此研究了船舶六自由度运动对破冰模式的影响，并通过实尺度破冰船数据进行了对比验证；王钰涵等[11]采用船冰接触过程理想化、挤压-弯曲过程理想化和海冰破损几何形状理想化等理想化破冰假设，对直航情况下连续破冰过程进行模拟，获得了冰力时历曲线和船舶运动响应，分析了相关破冰参数对连续破冰作用下破冰形状和平均冰力的影响程度。

在上述冰阻力数值预报方法中，计算模型通常考虑船冰间相互作用，采用船冰水耦合求解的研究相对较少。开展船冰水耦合作用数值计算，能够获得破冰过程中冰层破坏模式、冰块运动行为、船舶破冰排冰现象等更多信息，可为船舶性能预报和船型设计优化提供更为直接有效的支撑，具有重要的理论意义和实际工程价值。本文基于有限元/光滑粒子流体动力学耦合算法，初步实现了破冰过程中的船/冰/水耦合过程模拟，并对船舶冰阻力进行了预报。

1 数值模型

海冰采用弹塑性模型，塑性屈服前海冰的力学行为描述为

式中：σ 为海冰应力；K = E/3( 1- 2υ )为海冰体积弹性模量；G =E/3( 1+ 2υ )为海冰剪切弹性模量；E 为海冰弹性模量；υ 为海冰泊松比；Ddel为膨胀应变；D′e= De- DeI/2 为弹性应变偏量；De为海冰弹性应变。塑性模型屈服准则采用修正的Drucker-Prager 帽盖模型。船舶破冰航行过程中，冰层内发生应力集中产生裂纹，随着裂纹的扩张，冰层出现断裂破坏现象。为模拟海冰裂纹的产生，引入粘聚单元损伤模型，即将冰层离散成孤立实体单元，在相邻单元间插入粘聚单元作为潜在的断裂面，如图1所示，达到临界强度时，粘聚单元失效，冰层内出现破坏，其中，实体单元采用三棱柱单元，粘聚单元采用零厚度六面体单元。粘聚单元应力与相对应变在达到黏结强度之前表现为线弹性：

图1 单元间粘聚方式示意图Fig.1 Illustration of the ice mesh topology

式中，t表示应力，δ表示位移，下标n表示法向，下标s和t表示两个切线方向。通过单元刚度矩阵K定义应力和位移间耦合关系，如果仅对Knn、Kss和Ktt参数定义，设定其他刚度系数值为0，即表示法向和切向为非耦合本构关系。

船舶破冰过程中，海水对海冰力学行为的影响不可忽视。本文采用光滑粒子流体动力学方法模拟海水环境，实现破冰航行过程的船冰水动态耦合求解。光滑粒子法是一种无网格的纯拉格朗日方法，采用一系列任意分布的粒子质点来代表连续介质流体，通过核函数W 来定义一定光滑长度h范围内其他邻近粒子对目标粒子的影响程度，则流体域内的速度、压力及其梯度等变量均可通过一组无序质点的核函数插值集合表示。

流体力学基本方程组离散为如下近似形式：

式中，dρi/dt、dUi/dt分别为粒子i相应的随体导数。

有限元网格与水粒子之间通过罚函数接触算法传递接触状态和作用力[12-13]，接触面的法向接触力表达式为

式中，ki是接触刚度，fs是ki的比例因子，li是粒子i 相对网格单元的穿透深度，ni是网格单元的法向单位矢量。接触面的切向接触力为fti= μ ||fs，其中μ是粒子与网格单元的摩擦系数。

计算得到耦合作用力后，将作用力以外力的形式加入到流体动力方程和有限元动力方程中：

图2 船/冰/水耦合求解过程Fig.2 The solution procedure of a ship-ice and water coupled problem

船冰水耦合求解方法如图2 所示，船体为刚体模型，船冰间进行接触判断并计算相互作用力，同时根据水粒子的分布和位置计算水载荷，得到冰单元应力，当应力达到损伤起始值时采用损伤模型进行单元破坏判定，当达到破坏极限时，粘聚单元失效，更新碎冰单元接触面使其能够参与新的接触判断进而计算碎冰单元运动，同时更新冰单元边界，如此迭代实现船冰水耦合求解。

2 海冰强度数值模拟

海冰压缩强度和弯曲强度是影响船舶冰阻力预报的主要参数，采用本文建立的海冰数值模拟方法，参考ITTC推荐规程[14]中海冰物理强度测试方法，建立相应尺寸海冰试验数值模型。采用相同加载方式，选取斯瓦尔巴群岛海域当年冰为目标海冰，根据实测统计结果[15]，选取冰样厚度为0.6 m，压缩强度为620 kPa，弯曲强度为470 kPa。分别开展单轴压缩强度和弯曲强度计算，验证本文建立的海冰数值模拟方法，海冰数值模型主要参数见表1。

表1 海冰材料参数表Tab.1 Material properties of ice

单轴压缩强度计算：参考物理模型试验方式，冰样尺寸为1h × 2h × 4h，h为冰样厚度，采用上下压盘对冰样进行固定和载荷施加，模拟计算的海冰破坏对比及载荷时历曲线如图3 所示，可以得出：海冰破坏形式与试验相符，采用σc= F/( Wh )计算得到海冰压缩强度为637.5 kPa，与实测统计压缩强度（620 kPa）相比，数值模拟结果偏大2.8%。

图3冰样单轴压缩压缩试验数值模拟Fig 3 Numerical simulation of the uniaxial compression test of ice sample

弯曲强度计算：采用三点弯曲方式验证弯曲强度，梁的尺寸为1h × 2h × 6h，采用三个圆柱体对冰样进行固定和载荷施加，模拟计算的海冰破坏对比及载荷时历曲线如图4 所示，可以得出：海冰破坏形式与试验相符，采用σf= M/W = 3Fl/( )2bh2计算得到海冰弯曲强度为491.3 kPa，与实测统计弯曲强度（470 kPa）相比，数值模拟结果偏大4.5%。

图4 冰样三点弯曲试验数值模拟Fig.4 Numerical simulation of the three point bending test of ice sample

综上，采用本文建立的海冰数值模拟方法，建立海冰数值模型，开展冰样压缩破坏和弯曲破坏数值模拟，海冰破坏形式与试验相符，弯曲强度和压缩强度与实测统计结果相差分别为2.8%和4.5%，海冰数值模型可作为开展船舶破冰航行过程模拟和冰阻力计算的输入。

3 船舶冰阻力计算分析

本文以某极地船舶为研究对象，船型参数见表2，船体外形如图5所示。

表2 船舶主尺度Tab.2 Principal dimensions of the carrier

图5 船体外形图Fig.5 Outside view of the ship

根据目标海域冰况特点，采用建立的海冰数值模型，选取1.3 m 和1.6 m 两种厚度冰况，采用艏破冰方式，对船舶破冰过程进行计算，计算得到破冰过程和现象如图6～7 所示。船舶艏部采用挤压方式破冰前进，船舶前进方向，冰层出现径向裂纹，随着航行前进，冰层断裂表现为环向裂纹，进而断裂，破冰现象在下一次接触破冰过程中循环出现；冰层破坏形成的冰块被排向船体两侧，少量冰块在浮力和摩擦力作用下沿船体滑移，部分冰块在流场扰动作用下漂浮至航道中；破冰和排冰接触作用构成了冰阻力的不同组成部分。

图7 冰层局部破坏现象Fig.7 Local damage of ice sheet

船舶破冰阻力时历曲线如图8 所示，初始破冰过程中，船舶靠近冰层利用艏部挤压进行破冰，船冰接触面积不断增加，冰阻力逐渐增大；而后冰层受到持续破坏，粘聚单元持续失效造成接触力卸载，冰阻力时历曲线表现为高频振荡（曲线a），通过对阻力信号进行频谱分析，发现其能量谱峰值主要集中在100 Hz 以下低频范围内，高频区域对整体冰阻力的影响较为有限，采用低通滤波方式滤除高频噪声，获得较为直观的阻力（曲线b）。同时，采用Lindqvist公式[1,16]估算船舶冰阻力，假定70%的船体湿表面积被碎冰块覆盖，得到两种冰厚下阻力曲线，如图9所示。

通过分析船舶破冰现象和阻力变化曲线可知：（1）随着航速和冰厚增加，船舶阻力明显增大；数值计算结果中，阻力点与阻力-航速曲线偏差较大。数值计算考虑海冰的弯曲和剪切破坏，海冰破坏存在随机性，造成了阻力点随航速变化的离散性，使阻力点偏差增大。（2）对比经验公式与数值计算结果，1.3 m 厚度冰况阻力相差范围为7.32%～13.91%，1.6 m 厚度冰况阻力相差范围为5.13%～8.07%，数值计算得到的冰阻力均大于经验公式预报结果。Lindqvist 公式计算冰阻力时，航速与冰阻力是线性相关的，难以考虑冰水耦合作用及碎冰块运动的影响，而数值计算对船冰水耦合作用进行模拟，冰阻力结果较大。（3）经验公式与数值计算的冰阻力差别随着航速的增加而增加，随着层冰厚度的增加而减小。船舶破冰过程中，经验公式与数值计算的冰阻力差别主要由冰水耦合及碎冰运动引起；随着航速增加，冰层破坏后碎冰速度和碰撞程度均明显增加，冰阻力差别随航速增加明显；随着冰厚增加，冰块破碎后尺寸增大，由于惯性原因碎冰块运动受到海水影响减弱，引起冰阻力差别随厚度的增加而降低。

图8 船舶破冰阻力时历曲线Fig.8 Time history of ice-breaking resistance

图9 不同冰层厚度下破冰阻力和Lindqvist结果比较Fig.9 Comparison of the ice-breaking resistance of the polar vessel with Lindqvist’s results in ice with different thickness

4 结 语

本文基于有限元理论，采用粘聚单元法构建海冰数值模型，开展了冰样单轴压缩强度和三点弯曲强度试验过程数值模拟，验证了该数值模型的合理性。在此基础上，采用有限元/光滑粒子流体动力学耦合算法，建立船舶破冰过程中船/冰/水耦合作用数值计算方法，模拟破冰过程裂纹产生和扩展、冰块翻转与滑移等现象，针对某极地运输船舶，开展了冰阻力数值计算和分析，得到不同航速、冰厚状态下冰阻力变化情况。该方法可为船舶连续破冰模式下阻力预报提供一种技术手段，可为浮/碎冰工况数值计算方法研究提供参考。