冰区加强型集装箱船碎冰航道航行阻力数值模拟

张远双，齐江辉，郑亚雄，吴述庆

(1.武汉船舶职业技术学院，湖北 武汉 430050；2.武汉第二船舶设计研究所，湖北 武汉 430064)

0 引 言

冰区加强船舶不需要自主破冰，通常需要在破冰船已经开辟的航道中航行，即在碎冰区域航行。船舶在冰水混合的碎冰区航行时，航行阻力与开敞水域有明显不同。关于船舶在冰水混合流动中的阻力性能预报，通常有理论计算、经验公式估算、数值模拟和船模试验等方法。理论计算及经验公式估算受船型尺度等限制，通常只针对特定船型且精度较差。

在数值模拟方面，国内外开展的相关研究并不多。目前应用最多的是离散元法（DEM）、有限元法（FEM）以及光滑粒子流体动力学方法（SPH）等。Jungyong Wang 等[1]基于Ls-dyna 软件模拟了碎冰条件下Terrry Fox 号破冰船的阻力性能，分析了碎冰密集度对阻力性能的影响。Moon-Chan Kim 等[2]对一艘冰区加强型散货船的碎冰区阻力性能进行了数值模拟，通过与试验结果的对比验证了数值结果的准确性，其数值结果与试验值误差在20%左右。郭春雨等[3]同样采用Ls-dyna 软件计算各种工况下的船舶航行阻力特性，并与试验值进行了对比。王超等[4]基于离散元模型结合欧拉多相流对碎冰区船舶的冰阻力进行了数值模拟，得到了碎冰阻力随航速等的变化规律。在自编程计算方面，李紫麟[5]基于离散元理论，采用三维圆盘模拟海冰单元，对不同航速碎冰作用于船体的力进行了数值模拟。J.Tuhkuri 等[6]对DEM 方法应用于冰水船相互作用的研究进行了系统的评述。总的来说，碎冰区域航行船舶的阻力性能数值预报仍需要提高精度。

在船模试验方面，冰水池船舶试验通常分为冻结模型冰试验和非冻结模型冰试验。随着冰水池技术的发展，国内外冰池船模试验技术迅速发展，并取得了较多研究成果。首次破冰船阻力试验是在1964 年，Corlett 和Snaith[7]以石蜡代替冰，对1 艘小型破冰船进行了阻力试验。随后德国、芬兰、韩国、加拿大、俄罗斯等冰水池相继在船舶冰区阻力、操纵性等方面开展了大量的模型试验。在国内，郭春雨等[8–9]依托哈尔滨工程大学常温拖曳水池采用非冻结模型冰进行了一系列的船-冰相互作用试验。黄焱等[10–11]在天津大学冰力学实验室，开展了碎冰和平整冰阻力试验。

本文采用STARCCM+商用软件，基于离散元方法，结合拉格朗日方法的颗粒接触动力学模型，建立船-冰-水的CFD-DEM 耦合数值计算模型，对一艘冰区加强型集装箱船在碎冰区航行的阻力特性进行研究。参照冰水池试验结果，选取不同的颗粒形状进行数值模拟，将数值计算现象与试验现象进行对比，对不同航速、不同冰厚下的阻力特性进行分析。

1 理论基础

1.1 数值计算方法

计算模型中，流体满足不可压缩及连续性方程，同时忽略各相之间的热传递过程[12]。基于VOF 方法捕捉自由液面，湍流模型为SST-K-omege 模型。为防止出流边界处的流体反射对船体周围流场的影响，对出流边界附近采用Choi 方法进行数值消波[13]。数值离散方法采用离散元方法。

离散元方法在20 世纪70 年代初由Cundalland Strack[14]首次提出。离散元方法将研究对象离散成一系列的单元，单元之间通过节点进行连接。其求解过程一般为：1）将研究对象离散为单元颗粒，选定合理的单元连接形式，给定初始条件；2）根据接触条件判断单元之间是否发生接触，若发生接触则根据单元间的本构模型，计算单元之间相互作用力，从而得到各单元之间的相对速度和位移，更新单元信息。通过对计算域内所有单元的物理量的求解，可以得到整体的研究对象的变形及运动情况。

1.2 离散项接触力模型

两单元之间接触作用模型如图1 所示，在离散元中接触力公式实际是弹簧-阻尼器模型的一种变形。弹簧产生排斥力将颗粒分开，而阻尼器表示粘性阻尼并允许模拟除完全弹性外的碰撞类型。作用在接触点处的接触力可以看做是一对弹簧-阻尼器振子。

图1 接触力模型Fig.1 Model of contact force

Hertz-Mindlin 无滑移接触模型是一种非线性弹簧-阻尼器接触模型的变形[15]。作用在2 个单元A 和B 之间的力可表示为：

其中： Fn为法向力分量； Ft为切向力分量。

法向力分量可以写成：

切向力分量可以写成：

其中：Nndamp和 Ntdamp分 别为法向阻尼系数和切向阻尼系数，Cnrest和Ctrest分 别为法向和切向的恢复系数，若Cnrest=0 则 Nndamp=1，若 Ctrest=0 则 Nndamp=1； Req， Meq和Geq分别对应等效半径、等效杨氏模量和等效剪切模量。

对于单元-壁面碰撞，上述公式保持不变，将壁面半径和质量设定为 Rwall=∞ 和 Mwall=∞，因此等效半径为 Req=Rparticle， 等效质量 Mwall=Mparticle。

2 数值计算流程

2.1 离散元模型

离散项冰粒子的建模参照冰水池试验的碎冰图像，碎冰粒子表现为一系列的圆球组成的复合粒子，碎冰粒子形状选取为金字塔形和不规则棱柱形，如图2所示。冰的密度为917 kg/m3。冰的弹性模量为E=9 GPa，根据冰属性的缩尺比关系[16]，冰粒子模型的弹性模量为300 MPa，泊松比为0.3。

图2 碎冰几何形状Fig.2 Geipetric shape of crushed ice

图3 数值航道中的碎冰粒子Fig.3 Ice particles in numerical channels

根据北极地区碎冰尺寸统计数据，碎冰尺寸大致服从对数正态分布规律，其分布形式为：

本文数值模拟中，颗粒（2 种碎冰模型）、流体（水）和壁面（船体壁面、航道边界壁面）多相之间会产生碰撞等相互作用，在计算中考虑7 种相互作用类型，如图4 所示。其中，船-冰之间的摩擦系数为0.1，冰-冰之间的摩擦系数为0.1，动摩擦系数为0.3。冰粒子采用喷射器的形式进入流场，喷射器定义冰粒子进入的速度、位置及方式等。

2.2 计算模型设置

本文计算模型为1 艘典型冰区加强型集装箱船，模型缩尺比为30，实船及模型主尺度参数如表1 所示。

图4 多相相互作用示意图Fig.4 Schematic diagram of mutiphase interaction

表1 船体主尺度参数Tab.1 The main parameters of ship

碎冰航道计算域以船体为中心，船体前后距计算域边界距离均为8 L，船体两侧距计算域边界距离均为1.5 Lpp，在垂直方向上计算域尺寸为−1.5 Lpp～2.0 Lpp。船体、平整冰及计算域的4 个壁面均为无滑移壁面边界条件，速度入口及压力出口边界如图5 所示。本文计算域分为空气域和水域，计算过程中船体保持不动，通过改变速度入口处的速度调节船舶的航速。为准确捕捉自由液面位置、船体周围流场等，在相应的区域进行网格加密处理，边界层网格形式为棱柱形边界层，边界层采用y+壁面处理，y+值范围为30～60。船体表面网格尺度为5 mm，船体表面第1 层网格厚度为2 mm，计算域中等网格总数约为212 万，图6 为网格划分结果。

图5 计算域及边界条件Fig.5 Calculation domain and boundary conditions

图6 网格划分结果Fig.6 Result of meshing

3 数值计算结果及分析

3.1 船冰接触验证

图7 为汉堡冰水池进行的某冰区加强型散货船的碎冰阻力试验图像[17]。船舶在碎冰区域航行时，碎冰在与船舶发生碰撞时会发生翻转现象，而随着船舶的航行，在船体首部区域会发生碎冰堆积现象，同时一部分碎冰运动到船体底部或舷侧随船体滑移。在这个过程中，碎冰之间以及碎冰与船体之间会发生碰撞。从图7 可以看出数值计算现象与试验现象吻合度较好，证明本文采用的碎冰粒子模型及数值方法可以有效模拟碎冰区船舶与碎冰接触时的现象。

图7 船-冰碰撞现象模拟Fig.7 Simulation of ship-ice collision

3.2 碎冰运动状态分析

船舶在碎冰区域航行时，碎冰与船体之间的碰撞、摩擦等作用产生的阻力是船体阻力的重要成分。分析船体周围碎冰的运动状态对阻力特性的预报十分重要。计算V=0.563 m/s 时船体周围碎冰分布如图8 所示。可以看出碎冰在船体首部发生碰撞后堆积，碎冰与船体间的相对速度骤降。而在船体尾部，由于船体尾流的作用使得碎冰分布较为分散，聚集在螺旋桨周围的碎冰较少，这可以有效减少冰-桨接触，有利于螺旋桨的推进效率。从图9 也可以看出，随着航速的增加，船体周围的兴波已经较为明显，兴波可以有效地将船体周围的碎冰排开。同时可以看到，在船体尾部形成了一条浮冰很少的航道，该航道宽度大致相当于船体的宽度，这会有效降低船体-冰之间的接触力。

3.3 航速对碎冰阻力的影响

图8 船体周围碎冰的运动状态Fig.8 Movement position of crushed ice around hull

图9 不同航速碎冰状态对比Fig.9 Movement position of crushed ice in different speed

船舶在敞水区域航行时，航速对船体阻力的影响很大。计算航速分别为0.187 m/s，0.376 m/s，0.563 m/s，0.751 m/s 和0.939 m/s（分别对应实船航速2 kn，4 kn，6 kn，8 kn 和10 kn）时，船体受到的冰阻力时历曲线如图10 所示。时历曲线中阻力分为3 部分，contact force 为船-冰接触力，resistance 为水阻力，net force 为船体受到的总阻力。可以看出船-冰接触力随着航速的增大而增大，同时接触力有强烈的震荡特性。接触力的3 个分量中，纵向力以某一基准值为下限震荡；侧向力则在0 值附近上下震荡，这是由于船体左右两舷对称的原因；而垂向力基本为0，这表明船-冰碰撞基本为水平方向碰撞。低速航行时，纵向接触力震荡下限远大于水阻力；航速较高时，纵向接触力的震荡下限小于水阻力。这说明随着航速的增大，水阻力在船体总阻力中所占的比重越来越高。

对船体阻力的各个分量取稳定段时间平均即得到各航速下的船体阻力分量如图11 所示。可以看出，船体总阻力及各分量均随航速的增大而增大，这与敞水航行阻力的规律类似。在航速从0.187 m/s 增大到0.939 m/s 时，航速变为原来的5 倍，总阻力增大为原来的3.5 倍，纵向接触力增大为原来的2.2 倍，水阻力增大为原来的20.4 倍，说明随着航速的增大，水阻力的增大速度远高于船-冰接触力的增大速度。但在航速增大的过程中，纵向接触力占总阻力的比例越来越小，航速为0.939 m/s 时比重约为59%；水阻力占总阻力的比例越来越大，航速为0.939 m/s 时比重约为41%。

图10 低速和高速时船体阻力分量时历曲线Fig.10 Time history curve of hull resistance at low speed and high speed

图11 船体总阻力及分量随航速变化曲线Fig.11 Comparision of resistance components

图12 不同航速船体表面接触力分布Fig.12 Contace force distribution of hull surface in different speed

为进一步解释航速增大而船-冰接触力增大不明显的现象，图12 给出2 种航速下船体表面与碎冰接触力的分布。可以看出，在低航速时，碎冰在船体周围聚集使得船体与碎冰的接触面积较大；而高航速时，高航速及船体的兴波等使得船体排开碎冰的效果更明显，船体周围聚集的碎冰数量明显减少，因此船体与碎冰的接触面积变小，因此航速增大后接触力的增长不明显。

4 结 语

本文以1 艘典型冰区加强型集装箱船为研究目标，建立DEM 和CFD 相结合的数值模拟方法对碎冰区航行船舶的船-冰-水耦合作用进行研究。讨论了船体周围碎冰的运动状态、航行阻力特性分析以及航速对船体航行阻力的影响等，结论如下：

1）基于CFD-DEM 结合的数值方法可以有效模拟船舶在碎冰区航行的过程，对碎冰在船体附近的堆积、翻转等运动状态可以有效模拟，采用该方法进行船-冰-水耦合作用分析是有效可行的。

2）碎冰在与船体接触作用后会在首部形成一个明显的减速区，使得碎冰在首部堆积。首部的碎冰一部分沿着船体侧边或底部滑行。在尾部区域，由于尾流的作用会在船体尾部形成一条碎冰数量较少的航道，而随着航速的增大，该航道碎冰变得更少，同时船体左右两侧的兴波使得船体排开碎冰的作用更加明显。

3）航速对碎冰区航行船舶阻力有重要影响，碎冰区航行阻力与敞水航行阻力曲线特性一致。随着航速的增大，水阻力的增大速度远高于船-冰接触力的增大速度，水阻力占总阻力的比重随着航速的增大而增大明显。