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

基于FEM/SPH 耦合方法的极地
运输船舶冰阻力预报研究
（中国船舶科学研究中心，江苏无锡214082）
摘要：冰阻力是极地运输船舶动力需求分析和船型优化设计阶段的重要参考指标。

本文基于有限元理论，采用粘聚单元法构建海冰数值模型，开展冰样单轴压缩和三点弯曲试验过程模拟，验证了该模型满足海冰数值计算的要求。

在此基础上，以某极地运输船舶为研究对象，采用有限元/光滑粒子流体动力学耦合算法，通过罚函数接触算法传递有限元网格与水粒子之间接触状态和作用力，考虑船舶破冰过程中船/冰/水耦合作用，模拟破冰过程裂纹产生和扩展、冰块翻转与滑移等现象，并进行了冰阻力预报。

关键词：极地运输船舶；粘聚单元；海冰数值模型；破冰阻力
中图分类号：U661.1文献标识码：A doi:10.3969/j.issn.1007-7294.2021.01.002
Research on ice resistance prediction of a polar vessel based on FEM/SPH coupling method
WANG Zhi-peng ，HAO Zhai-liu ，TIAN Yu-kui ，WU Bao-shan (China Ship Scientific Research Center,Wuxi 214082,China)
Abstract ：Ice resistance is an important reference index for the marine power requirement analysis and ship form optimization design of polar vessels.In this paper,the cohesive element method is used to construct the ice numerical model based on the theories of finite element method,and the numerical simulations are car⁃ried out on uniaxial compression test and three-point bending test of the samples,the calculation results veri⁃fy that the model meets the requirement of the numerical simulation.On this basis,the FEM/SPH coupling al⁃gorithm is adopted to consider the coupling effect of ship/ice/water when a polar vessel is on an ice-break⁃ing voyage,and the contact states and forces between finite element meshes and particles of water are calcu⁃lated by the contact algorithm of penalty functions.The result shows the initiation and propagation of cracks on ice sheets,rolling over and sliding of the ice blocks while the ice resistance is predicted.Key words:polar vessel ；cohesive element ；numerical model of ice material ；icebreaking resistance 0引言
北极航道是联系亚欧美三大洲的最短航线，与通常的苏伊士运河和巴拿马运河航线相比，能够缩短数千公里航程，由此将大大降低商业运输成本，影响和意义深远。

随着气候变暖，海冰融化加速，北极航道通航条件不断改善，未来对极地运输船舶的需求会明显增大。

极地运输船舶不仅要满足开敞第25卷第1期船舶力学Vol.25No.12021年1月
Journal of Ship Mechanics Jan.2021
文章编号：1007-7294（2021）01-0009-07收稿日期：2020-07-19
基金项目：水动力学国防科技重点实验室基金（61422030209162203002）；工信部高技术船舶科研项目（[2017]614）作者简介：王志鹏（1987-），男，硕士，工程师，通讯作者，E-mail:********************。

水域性能和经济指标，而且需具备一定的破冰能力以满足冰区航行需求。

在极地运输船舶开发过程中，冰阻力预报是船舶破冰能力和动力需求分析的重要输入，也是总体设计的关键环节之一。

在船舶冰阻力预报方面，主要有模型试验、经验公式、数值模拟等方法，其中，模型试验是最可靠的方式，基于实船测量和模型试验数据总结的经验公式也被广泛应用，目前常用的经验方法有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数值模型
海冰采用弹塑性模型，塑性屈服前海冰的力学行为描述为
σ=KD del I +2GD′e (1)
式中：
σ为海冰应力；K =E /3()1-2υ为海冰体积弹性模量；G =E/3()1+2υ为海冰剪切弹性模量；
E 为海冰弹性模量；υ为海冰泊松比；D del 为膨胀应变；D′e =D e -D e I /2为弹性应变偏量；D e 为海冰弹性应变。

塑性模型屈服准则采用修正的Drucker-Prager 帽盖
模型。

船舶破冰航行过程中，冰层内发生应力集中产生裂纹，随
着裂纹的扩张，冰层出现断裂破坏现象。

为模拟海冰裂纹的产
生，引入粘聚单元损伤模型，即将冰层离散成孤立实体单元，在相
邻单元间插入粘聚单元作为潜在的断裂面，如图1所示，达到临界
强度时，粘聚单元失效，冰层内出现破坏，其中，实体单元采用三棱柱单元，粘聚单元采用零厚度六面体单元。

粘聚单元应力与相10船舶力学第25卷第1
期图1单元间粘聚方式示意图Fig.1Illustration of the ice mesh topology
对应变在达到黏结强度之前表现为线弹性：
{}t n t s t t =éëêêùûúúK nn K ns K nt K ns K ss K st K nt K st K tt {}
δn δs δt (2)式中，t 表示应力，δ表示位移，下标n 表示法向，下标s 和t 表示两个切线方向。

通过单元刚度矩阵K 定义应力和位移间耦合关系，如果仅对K nn 、K ss 和K tt 参数定义，设定其他刚度系数值为0，即表示法向和
切向为非耦合本构关系。

船舶破冰过程中，海水对海冰力学行为的影响不可忽视。

本文采用光滑粒子流体动力学方法模拟海水环境，实现破冰航行过程的船冰水动态耦合求解。

光滑粒子法是一种无网格的纯拉格朗日方法，采用一系列任意分布的粒子质点来代表连续介质流体，通过核函数W 来定义一定光滑长度h 范围内其他邻近粒子对目标粒子的影响程度，则流体域内的速度、压力及其梯度等变量均可通过一组无序质点的核函数插值集合表示。

f ()x ≃∑j m j ρj f i W ()||x -x j ,h (3)流体力学基本方程组离散为如下近似形式：
d ρi d t =∑j =1N m j ρi ρj
()U i -U j ∇i W ij (4)d U i d t =-∑j =1N m j ()
p i ρ2i +p j ρ2j ∇i W ij (5)式中，
d ρi /d t 、d U i /d t 分别为粒子i 相应的随体导数。

有限元网格与水粒子之间通过罚函数接触算法传递接触状态和作用力[12-13]，接触面的法向接触力表达式为
f s i =-f s k i l i n i (6)
式中，k i 是接触刚度，f s 是k i 的比例因子，l i 是粒子i 相对网格单元的穿透深度，n i 是网格单元的法向单位矢量。

接触面的切向接触力为f ti =μ||f s ，其中μ是粒子与网
格单元的摩擦系数。

计算得到耦合作用力后，将作用力以外力的形式加入到流体动力方程和有限元动力方程中：
d U i d t =-∑j =1N m j ()p i ρ2i +p j ρ2j ∇i W ij +F ext (7)Mx n +Cx n +Kx n =F +F ext (8)式中：
M 为质量矩阵；C 为阻尼矩阵；K 为刚度矩阵；x n 是各节点的加速度矩阵；x n 是各节点的速度矩阵；F 是节点的外载荷矩阵；F ext 是作用力。

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

第1期王志鹏等：基于FEM/SPH 耦合方法的极地 (11)
Fig.2The solution procedure of a ship-ice and water coupled problem
2海冰强度数值模拟
海冰压缩强度和弯曲强度是影响船舶冰阻力预报的主要参数，采用本文建立的海冰数值模拟方法，参考ITTC 推荐规程[14]中海冰物理强度测试方法，建立相应尺寸海冰试验数值模型。

采用相同加载方式，选取斯瓦尔巴群岛海域当年冰为目标海冰，根据实测统计结果[15]，选取冰样厚度为0.6m ，压缩强度为620kPa ，弯曲强度为470kPa 。

分别开展单轴压缩强度和弯曲强度计算，验证本文建立的海冰数值模拟方法，海冰数值模型主要参数见表1。

表1海冰材料参数表
Tab.1Material properties of ice
名称
符号
单位
数值海冰密度ρkg/m 3910弹性模量E Pa 109
泊松比υ-0.3粘聚单元刚度K Pa/m 1.18×107损伤起始应力T i Pa 2×105海冰-钢间摩擦系数μice -0.1
单轴压缩强度计算：参考物理模型试验方式，冰样尺寸为1h ×2h ×4h ，h 为冰样厚度，采用上下压盘对冰样进行固定和载荷施加，模拟计算的海冰破坏对比及载荷时历曲线如图3所示，可以得出：海冰破坏形式与试验相符，采用σc =F/()Wh 计算得到海冰压缩强度为637.5kPa ，与实测统计压缩强度（620kPa ）相比，数值模拟结果偏大2.8%。

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

12船舶力学第25卷第1期图3冰样单轴压缩压缩试验数值模拟
Fig 3Numerical simulation of the uniaxial compression test of ice
sample t /s
F /104
N
t /s
F /105
N
图4冰样三点弯曲试验数值模拟
Fig.4Numerical simulation of the three point bending test of ice sample
综上，采用本文建立的海冰数值模拟方法，建立海冰数值模型，开展冰样压缩破坏和弯曲破坏数值模拟，海冰破坏形式与试验相符，弯曲强度和压缩强度与实测统计结果相差分别为2.8%和4.5%，海冰数值模型可作为开展船舶破冰航行过程模拟和冰阻力计算的输入。

3船舶冰阻力计算分析
本文以某极地船舶为研究对象，船型参数见表2，船体外形如图5所示。

表2船舶主尺度
Tab.2Principal dimensions of the carrier
船模主尺度
符号
单位
数值总长L OA m 279垂线间长L pp m 259水线宽B WL m 44吃水T m 15.89排水体积▽m 3
144103
根据目标海域冰况特点，采用建立的海冰数值
模型，选取1.3m 和1.6m 两种厚度冰况，采用艏破
冰方式，对船舶破冰过程进行计算，计算得到破冰
过程和现象如图6~7所示。

船舶艏部采用挤压方
式破冰前进，船舶前进方向，冰层出现径向裂纹，随着航行前进，冰层断裂表现为环向裂纹，进而断裂，
破冰现象在下一次接触破冰过程中循环出现；冰层
破坏形成的冰块被排向船体两侧，少量冰块在浮力
和摩擦力作用下沿船体滑移，部分冰块在流场扰动
作用下漂浮至航道中；破冰和排冰接触作用构成了冰阻力的不同组成部分。

第1期王志鹏等：基于FEM/SPH 耦合方法的极地 (13)
图5船体外形图
Fig.5Outside view of the ship 图6船舶破冰过程
Fig.6Numerical result of the polar vessel crushing ice sheet 图7冰层局部破坏现象Fig.7Local damage of ice sheet (a)径向裂纹(b)环向裂纹
(c)冰块翻转
(d)碎冰滑移
船舶破冰阻力时历曲线如图8所示，初始破冰过程中，船舶靠近冰层利用艏部挤压进行破冰，船冰接触面积不断增加，冰阻力逐渐增大；而后冰层受到持续破坏，粘聚单元持续失效造成接触力卸载，冰阻力时历曲线表现为高频振荡（曲线a ），通过对阻力信号进行频谱分析，发现其能量谱峰值主要集中在100Hz 以
下低频范围内，高频区域对整体冰阻力
的影响较为有限，采用低通滤波方式滤
除高频噪声，获得较为直观的阻力（曲
线b ）。

同时，采用Lindqvist 公式[1,16]估算
船舶冰阻力，假定70%的船体湿表面积
被碎冰块覆盖，得到两种冰厚下阻力曲
线，如图9所示。

通过分析船舶破冰现象和阻力变化曲线可知：（1）随着航速和冰厚增加，船舶阻力明显增大；数值计算结果中，阻力点与阻力-航速曲线偏差较大。

数值计算考虑海冰的弯曲和剪切破坏，海冰破坏存在随机性，造成了阻力点随航速变化的离散性，使阻力点偏差增大。

（2）对比经验公式与数值计算结果，1.3m 厚度冰况阻力相差范围为7.32%~13.91%，1.6m 厚度冰况阻力相差范围为5.13%~8.07%，数值计算得到的冰阻力均大于经验公式预报结果。

Lindqvist 公式计算冰阻力时，航速与冰阻力是线性相关的，难以考虑冰水耦合作用及碎冰块运动
的影响，而数值计算对船冰水耦合作用进行模拟，
冰阻力结果较大。

（3）经验公式与数值计算的冰阻
力差别随着航速的增加而增加，随着层冰厚度的
增加而减小。

船舶破冰过程中，经验公式与数值
计算的冰阻力差别主要由冰水耦合及碎冰运动引
起；随着航速增加，冰层破坏后碎冰速度和碰撞程
度均明显增加，冰阻力差别随航速增加明显；随着
冰厚增加，冰块破碎后尺寸增大，由于惯性原因碎
冰块运动受到海水影响减弱，引起冰阻力差别随
厚度的增加而降低。

4结语
本文基于有限元理论，采用粘聚单元法构建海冰数值模型，开展了冰样单轴压缩强度和三点弯曲强度试验过程数值模拟，验证了该数值模型的合理性。

在此基础上，采用有限元/光滑粒子流体动力学耦合算法，建立船舶破冰过程中船/冰/水耦合作用数值计算方法，模拟破冰过程裂纹产生和扩展、冰块翻转与滑移等现象，针对某极地运输船舶，开展了冰阻力数值计算和分析，得到不同航速、冰厚状态下冰阻力变化情况。

该方法可为船舶连续破冰模式下阻力预报提供一种技术手段，可为浮/碎冰工况数值计算方法研究提供参考。

参考文献：
[1]Lindqvist G.A straightforward method for calculation of ice resistance of ships[C]//Proceedings of the Tenth International 14船舶力学第25卷第1期
0510
15F /103k N t /s 图8船舶破冰阻力时历曲线
Fig.8Time history of ice-breaking resistance 图9不同冰层厚度下破冰阻力和Lindqvist 结果比较Fig.9Comparison of the ice-breaking resistance of the polar ves-sel with Lindqvist ’s results in ice with different thickness
第1期王志鹏等：基于FEM/SPH耦合方法的极地 (15)
Conference on Port and Ocean Engineering under Arctic Conditions.Lulea,Sweden,1989:722-735.
[2]Spencer D,Jones S J.Model-scale/full-scale correlation in open water and ice for Canadian coast guard"R-Class"ice⁃breakers[J］.Journal of Ship Research,2001,45(4):249-261.
[3]Keinonen A,Browne R P.Icebreaker performance prediction[C]//Proceedings of the First International Offshore and Polar Engineering Conference.Edinburgh,Scotland.1991:562-570.
[4]Riska K.Performance of merchant vessels in ice in the Baltic[M].Sjöfartsverket,1997.
[5]Arne Gürtner,Morten Bjerkås.Numerical simulation of ice action to a lighthouse[C].Proceedings of the ASME28th Interna⁃tional conference on Ocean,Offshore and Arctic Engineering，2009.
[6]Wang J,Derradji-Aouat A.Numerical prediction for resistance of Canadian icebreaker CCGS Terry Fox in level ice[C]//IC⁃SOT2009,International Conference on Ship and Offshore Technology.Busan,Korea,2009:9-15.
[7]任奕舟,邹早建.破冰船在冰层中连续破冰过程的数值模拟[J],振动与冲击.2016，35(18)：210-223.
Ren Y Z,Zhou Z J.Numerical simulation of the continuous icebreaking by an icebreaker in level ice[J].Journal of Vibration and Shock,2016，35(18):210-223.
[8]Su B,Riska K,Moan T.A numerical method for the prediction of ship performance in level ice[J].Cold Regions Science and Technology,2010,60(3):177-188.
[9]季顺迎,李春花,刘煜.海冰离散元模型的研究回顾及展望[J].极地研究,2013,24(4):315-330.
Ji S Y,Li C H,Liu Y.A review of advance in sea-ice discrete element models[J].Chinese Journal of Polar Research,2013, 24(4):315-330.
[10]Tan X,Su B,Riska K,et al.A six-degrees-of-freedom numerical model for level ice-ship interaction[J].Cold Regions Sci⁃ence and Technology,2013,92:1-16.
[11]王钰涵,李辉,任慧龙,等.连续破冰模式下破冰船的冰力研究[J].海洋工程，2013,4:68-73.
Wang Y H,Li H,Ren H L,et al.Study of ice force about icebreaker based on continuous breaking parttern[J].The Ocean Engineering,2013,4:68-73.
[12]张志春,强洪夫,高巍然.SPH-FEM接触算法在冲击力学数值计算中的应用[J].固体力学学报,2011,32(3):319-324. Zhang Z C,Qiang H F,Gao W R.Application of SPH-FEM contact algorithm in impact dynamics simulation[J].Chinese Journal of Solid Mechanics,2011,32(3):319-324.
[13]Wang S,Soares C G.Numerical study on the water impact of3D bodies by an explicit element method[J].Ocean Engineer⁃ing,2014,78:73-88.
[14]ITTC recommended procedures and guidelines,testing and extrapolation methods,ice testing,test methods for model ice properties[S].7.5-02-04-02,2014.
[15]Aleksandr Yurevich Iko.Offshore ice-resistant fixed platform for the Olginskoye field in the Pechora Sea[D].University of Stavanger,2014.
[16]Daniela Myland,Soren Ehers.Influence of bow design on ice breaking resistance[J].Ocean Engineering,2016,119:217-232.。