船舶数字碎冰场构建与冰阻力计算方法研究

金 强，张佳宁，张 雷，田于逵

（1.大连海事大学船舶与海洋工程学院，辽宁大连 116026；2.中国船舶科学研究中心，江苏无锡 214082）

0 引 言

近年来，随着全球气候变暖，北极航线受到更广泛的关注，北极资源的开采也成为众多环北极国家争相开展的目标，从而推动了极地船舶的建造与研究。对于极地船舶的研究开展最多的是对航行于层冰冰况下的船舶冰阻力性能预报，而对航行于极地碎冰区船舶冰阻力性能的探究却并不多见。在现已开通的北极航线中，近70%航段为碎冰所覆盖，船体所受冰阻力是冰区船设计要考虑的主要因素。

对于极地船舶冰阻力的研究，Jones 等[1]采用“Healy”号破冰船模拟层冰冰况下船舶破冰阻力性能试验，并将试验结果与实船航行数据进行对比。DuBrovin 等[2]在1950—1955 年间四次船舶冰池试验数据基础上总结推导出了针对预报船模冰池试验中所受冰阻力的计算公式。Kim 等[3]使用LS-DYNA软件模拟了具有破冰能力货船在碎冰航道下航行的冰阻力性能，并将结果与冰池试验结果进行了对比。Wang 等[4]研究了层冰冰况下“Terry Fox”号破冰船的破冰阻力以及自航试验结果，并将试验结果与实船航行资料进行了对比。Selvadurai等[5]运用二维离散元模型模拟移动层冰与冰区锚泊结构物间的相互影响。国内对于冰区船舶研究起步较晚，郭春雨等[6]采用非冻结模型冰试验方法研究了碎冰条件下船舶冰阻力随速度、密集度变化的规律，并采用Colbourne 实船转化方法预报了不同冰厚下的冰阻力。李紫麟[7]采用离散元方法，用三维圆盘模型模拟碎冰，将数值模拟结果与非冻结模型冰试验结果进行对比，分析影响冰阻力大小因素。本文根据冰池试验规程及相应的冰级规范要求，构建浮/碎冰的离散元数值模型及用于模拟计算的数值冰水池，采用离散元方法（discrete element method，DEM）数值计算航行于碎冰航道中船舶所受的冰阻力，并将数值模拟结果与HSVA冰池试验数据及Dubrovin经验公式方法的结果进行对比分析，进而验证数值模拟方法的精度。

1 碎冰航道船舶冰阻力数值计算

1.1 数值计算理论基础

离散元法方法由Cundall 提出，其核心为牛顿第二运动定律基本思想，将离散分布的介质看作是一系列离散单元的集合[8-9]。极地碎冰区的海冰具有较强的离散分布特性，可以将碎冰区的海冰视为一系列离散分布的单元集合。对于单个碎冰的离散元模型采用若干个离散粒子填充的方式创建出具有一定物理形状的碎冰模型，填充的基本粒子为球形颗粒。构成碎冰的粒子单元i和单元j间简化受力如图1 所示。碎冰-船体间接触时，假定与组成碎冰的球形单元间发生接触的船体表面视为半径R趋向于无穷大的球形单元。每个离散单元均满足相应的运动方程，一个单元求解的输出作为下一单元求解的输入，通过时间步迭代求解各个单元的运动方程进而求得碎冰的整体运动状态。

图1 碎冰单元i，j间接触示意Fig.1 Schematic diagram of contact between brash ice units of i and j

离散元算法的关键在于对单元间接触力的计算，构成碎冰的球形单元i和单元j间接触的等效物理模型如图2所示。

图2 碎冰单元i，j间接触的等效物理模型Fig.2 Equivalent physical model of contact between brash ice units of i and j

式中，Kt和Nt分别表示切向弹簧刚度和阻尼，Cfs为摩擦系数，dn和dt分别表示在单元接触点处法向和切向上的重叠量，vn和vt分别表示在接触点处的速度沿法向和切向上的速度分量。

式中，Nn,damp和Nt,damp分别表示法向和切向阻尼系数，定义如下：

式中，Cn,rest和Ct,rest分别为法向和切向恢复系数[9]，Req、Meq、Eeq及Geq分别表示相对于两接触单元i和j的半径、质量、杨氏模量以及剪切模量。

式中，vi和vj分别表示接触的海冰粒子单元i和j的泊松比；Ei和Ej分别表示相接触的海冰粒子单元i和j的杨氏模量。

1.2 船-冰数值模型及数字碎冰场构建

用于船模试验的某冰区油船数值模型如图3所示。

图3 某冰区油船数值模型Fig.3 Numerical model of a tanker

根据芬兰-瑞典冰级规范（FSICR）[14]中Class IA 冰级对碎冰冰厚的规定及试验中生成的海冰厚度确定数值计算中创建的碎冰厚度。创建碎冰离散元数值模型，首先绘制碎冰模型草图如图4所示，初步设定碎冰的尺寸形状为图4（a）中0.08 m×0.08 m的方形，图4（b）为基于图4（a）中草图绘制的不规则多边形。根据所要创建的模型冰冰厚尺寸，以冰厚为37.1 mm 为例对草图4（a）和（b）拉伸相应的厚度尺寸得到碎冰模型的物理实体，如图5所示。为在碎冰域中生成尺寸大小不同的碎冰模型，依据冰池试验中碎冰尺寸大小，对于生成的离散元碎冰模型采用罗辛-拉姆勒随机分布方法将模型冰尺寸在0.08 m×0.08 m 的基础尺寸上在[0.04，0.1]m 内随机分布，采用不同大小的球形基础颗粒对图5 中生成的物理模型进行粒子填充，填充粒子数为50。生成的碎冰离散元模型如图6所示，冰池试验与数值计算中碎冰的分布如图7和图8所示。

图4 创建的碎冰物理模型草图Fig.4 Sketch of greated physical model of brash ice

图5 生成的碎冰物理模型实体Fig.5 Generated brash ice physical model entity

图6 碎冰的离散元模型Fig.6 Discrete element model of the brash ice

图7 冰池试验中碎冰形状及分布Fig.7 Shape and distribution of brash ice in the ice tank

图8 数字碎冰场中碎冰分布Fig.8 Distribution of brash ice in the numerical brash ice field

网格的划分基础尺寸选为0.215 m。为提高计算精度，在数值碎冰场自由液面处叠加三层网格加密体，并且只在垂直方向（Z方向）上对网格尺寸进行细化。三层网格控制体网格尺寸依次取基础尺寸的10%、20%和40%，即三层网格控制体区域在垂直方向尺寸分别为0.021 5 m、0.043 m 和0.086 m。在船首尾部曲率变化较大部位以及船体两侧与碎冰相互作用的部位同样添加控制体对网格进行加密处理。船首尾部及船体表面水线面处的控制体设置如图9 所示，控制体区域内网格加密方式设置为各向同性，对X、Y、Z三个方向进行相同比例的加密处理，网格尺寸取为基础尺寸的10%，即0.021 5 m，同样对于船体两侧水线面处网格的加密尺寸取基础尺寸的10%。生成网格如图10和图11所示。

图9 船首尾及水线面处网格加密控制体Fig.9 Encryption control body of the grid at the bow,stern and waterline surface

图10 计算域中纵剖面处网格Fig.10 Grid at longitudinal section in the computational domain

图11 船体局部网格加密处理Fig.11 Encrypted grid of the corresponding parts on the ship

根据试验条件及ITTC冰池试验相关规定[15]，对于数字碎冰场的构造：船首前方为1.5Lpp，船尾后方为2.5Lpp，Lpp为船舶的垂线间长，碎冰航道宽度为两倍船宽。航道两侧创建壁面区域模拟试验中的层冰区域，整个冰场长、宽、高分别为46.0 m、16.0 m和10.0 m，如图12所示，初始化的数字碎冰航道见图13。

图12 数字碎冰场的三维尺寸Fig.12 Three-dimensional size of numerical brash ice field

图13 计算域初始化后的碎冰航道Fig.13 Brash ice channel after initialization of calculation domain

2 试验模型与方案设计

2.1 相似定律

冰池试验中为确保船模与实船的相似性，船模试验相似数（傅汝德数、雷诺数及柯西数）均需考虑，

式中，V表示速度，g表示重力加速度，L表示长度，ν表示运动粘度，ρ表示水密度，E表示杨氏模量。

由于在模型试验中傅汝德数与雷诺数不能够同时满足，且对于冰池试验而言流体粘性影响相对较小，因此冰模试验中通常忽略雷诺数相似。根据傅汝德数与柯西数相似规律，船体、碎冰及流体主要的几何及物理量相似关系满足表1所示，其中λ为缩尺比。

表1 各物理量相似规律Tab.1 The similarity laws of physical quantities

2.2 试验模型

试验采用某冰区油船为研究对象，测试冰况依据芬兰-瑞典冰级规范（FSICR）[14]选取Ice Class1A冰级条件，碎冰航道宽度为两倍船宽尺寸。船模首尾部如图14 所示（本章节所用船模试验说明图示均引用自文献[16]），船体主要参数见表2。

表2 某冰区油船主参数表Tab.2 Main parameters of a certain oil tanker

图14 试验模型船首尾部示意图Fig.14 Schematic diagram of the bow and stern of the model ship

2.3 冰池试验方案设计

2.3.1 冻结模型冰碎冰航道构造

冰池试验相对于一般敞水条件下的船模试验的技术更加复杂，碎冰航道的生成是试验前期准备最为耗时且技术要求复杂的工作。根据HSVA 标准模型冰制造程序生成厚度满足要求的层冰，在层冰冰厚达到设定要求后冰池室内温度开始上升，当温度接近-3℃时，切割出宽度等于两倍船宽大小的航道并采用冰凿等设备将航道内的层冰破碎构造初始碎冰航道。碎冰制备及生成的碎冰航道如图15和图16所示。

图15 碎冰制备Fig.15 Generation of the brash ice

图16 生成的碎冰航道Fig.16 Generated brash ice channel

2.3.2 试验工况

试验过程中流体设为静止，即不考虑实际冰况下碎冰漂移速度及水流速对船舶所受阻力的影响。将模型船航速设为定值0.464 m/s，将碎冰冰厚设置为37.1 mm、39.8 mm、46.3 mm 和49.4 mm，探究船舶所受冰阻力随冰厚的变化规律。试验中采用拖曳小车的方式使模型船匀速航行于碎冰航道中，如图17所示。

图17 碎冰航道中的模型船Fig.17 Model ship in the brash ice channel

3 数据分析与处理

3.1 数值模拟过程与冰池试验对比分析

数值模拟中采用相对坐标系的方式将船体视为静止，设置碎冰及水流以恒定速度V=0.464 m/s（实船航速5 kn）向船舶运动。以模型冰厚hice=39.8 mm 计算为例，碎冰航道中船舶航行如图18 和图19 所示。

图18 试验中碎冰航道中船舶航行状态Fig.18 Navigation status of the ship in brash ice channel in test

图19 数字碎冰场中船模航行状态Fig.19 Navigation status of the ship in numerical brash ice field

船首及船体两侧水线部位为船冰间的主要接触区域，航行于碎冰航道时船冰间主要发生三种相互作用形式：

（1）船首部位与碎冰间的碰撞、堆积与挤压，这也是船体所受冰阻力产生的主要因素。

（2）船首前方碎冰被不断推离至船体两侧，与船体两侧表面发生摩擦与挤压作用并在船尾后方留下等同于船宽大小的无碎冰分布航道，试验与数值计算中相应的现象分别如图20和图21所示。

图20 试验中船体后方碎冰航道Fig.20 Brash ice channel behind the ship in the test

图21 数字碎冰场中船尾后方碎冰航道Fig.21 Brash ice channel behind the ship in the numerical brash ice field

（3）碎冰自船首部位被压入船体以下，由于浮力作用碎冰紧贴船体底部表面运动直至船尾处浮出水面，试验与数值计算中相应的现象分别如图22和图23所示。

图22 冰池试验中压入船底的碎冰Fig.22 Brash ice pressed into the bottom of the ship in the test

图23 数字碎冰场中压入船底碎冰Fig.23 Brash ice pressed under the bottom of the ship in the numerical brash ice field

在数值模拟过程中，通过将船冰间三种典型的相互作用形式与冰池试验过程中相应的现象进行对比可知，采用粒子喷射的方式创建的数字碎冰场能够有效地模拟出船冰间的各种相互作用，并且与试验现象相吻合。

3.2 数值计算及经验公式结果与试验数据对比分析

数值模拟计算中为得到稳定结果，船模应航行至少两倍船长距离。选取t=200 s时刻不同冰厚条件下船舶所受阻力时历曲线如图24所示。

图24 不同冰厚条件下船体所受冰阻力时历曲线Fig.24 Time history curve of ice resistance on hull under different ice thicknesses

从图24中可以明显得出：

（1）不论在何种冰厚条件下，船舶所受冰阻力特点均是幅值随时间波动较大。这是因为作用于船体表面的冰阻力并不是一个恒定的外力作用，随着船舶在碎冰航道中的行进，船首部位与碎冰间的作用形式不断变化，发生碰撞的碎冰数量也随时间而改变。

（2）四种冰厚条件下，船舶所受冰阻力变化的总体趋势均随时间增长而增大，随后又趋于平稳。在船首初次接触碎冰到整个船体进入碎冰航道的过程中，船冰间发生相互作用的船体表面面积与碎冰数量不断增加，因而在这一过程中船舶所受冰阻力呈现不断增长趋势；而随着船体整个进入碎冰分布均匀的航道后，与船体表面发生碰撞的碎冰数及船体表面积均为定值，从而船舶所受冰阻力数值也保持稳定。

1970年DuBrovin提出针对碎冰试验中模型尺度下船舶所受阻力的经验公式[2]：

式中，R表示模型船所受冰阻力；p1与p2分别为系数，取决于碎冰分布密度及航道宽度，Fr为傅汝德数，n为功率系数，取决于船型，参数A和Φ定义[17]如下：

式中，r为碎冰分布密度，ρice为碎冰密度，f0为船冰间摩擦系数，α0和αH分别表示模型船水线面入射角的一半与船舶前体菱形系数。具体参数选取见表3。

表3 DuBrovin经验公式各参数选取Tab.3 The parameters for DuBrovin empirical formula

根据式（14）～（16）及表3 确定的参数可计算得出，在冰厚分别为37.1 mm、39.8 mm、46.3 mm 和49.4 mm条件下时，模型船相应地所受阻力大小为49.3 N、51.51 N、56.57 N和58.91 N。

对数值计算所得冰阻力数值进行处理，得出四种不同冰厚条件下船舶所受冰阻力均值，并将试验数据及经验公式所得结果整理，如表4所示。

表4 三种方法所得冰阻力数值Tab.4 Ice resistance values obtained by three methods

根据表4数据可知三种方法所得船舶冰阻力变化趋势如图25所示。

根据图25及表4数据分析，数值计算结果与DuBrovin经验公式所得的结果更加接近，但后者所得冰阻力数值普遍高于前者，且在冰厚为37.1 mm、39.8 mm、46.3 mm 和49.4 mm 四种条件下与冰池试验数据误差分别为9.24%、8.85%、24.33% 和25.38%，而数值计算所得结果在相应冰厚条件下与试验结果误差分别为12.60%、11.33%、28.64%和29.45%。从数据中可以看出，在冰厚尺寸较小条件下，经验公式方法与离散元算法计算所得阻力数值与试验数据误差较小，有较高的精度；随着冰厚增大，这两种计算方法的误差开始增大。这是由于对于经验公式方法而言，式（14）～（16）为DuBrovin 总结1950—1955 年间的四次船模冰池试验数据推导而得，原始数据相对较少并且从公式中可以看出冰厚因素hice对阻力计算并未起到主要影响作用，因此当冰厚增大时阻力计算的误差也相应增大；就离散元算法而言，船冰间的接触力主要取决于接触力Fc的分量，即法向力Fn与切向力Ft的计算，而从理论部分的式（1）可知，对于碎冰单元间接触力计算影响最大的为单元间法向与切向的刚度与阻尼系数，与碎冰厚度有关的质量M对接触力计算影响不大，因此随着冰厚增加所得的阻力值变化并不明显。

图25 三种计算方法所得阻力变化趋势Fig.25 Variation trend of resistance obtained by three methods

4 结 论

本文目的在于依据试验中冰池构造，采用粒子喷射的方式创建数字碎冰场，并且使用基于离散元算法的数值计算方法模拟船舶在冰厚为37.1 mm、39.8 mm、46.3 mm 和49.4 mm 四种条件下所受冰阻力，将计算结果及经验公式数据分别与模型船冰池试验数据进行对比分析，验证数值计算方法精度及DuBrovin经验公式方法估算船舶在碎冰区域所受冰阻力的可靠性，得出如下结论：

（1）对比冰池试验，在数值模拟中船冰间相互作用形式与试验中现象吻合一致。表明基于粒子喷射方式创建的数字碎冰场是合理的，能够准确地模拟出船舶航行于碎冰区船冰间的物理现象。

（2）使用DuBrovin 经验公式计算四种冰厚条件下的冰阻力，与试验结果相比较的误差分别为9.24%、8.85%、24.33%和25.38%，数据误差在一定的合理范围之内，从而表明该经验公式在冰阻力初步估算与预报中具有一定的准确度与参考意义。

（3）基于离散元算法的数值模拟方法与经验公式计算结果吻合较好，最大误差仅为5.70%；在冰厚尺寸较小条件下与冰池试验数据误差也在合理范围之内，即在模型冰厚为37.1 mm 和39.8 mm 时误差仅为12.60%和11.33%；而当模型冰厚增大为46.3 mm 和49.4 mm 时，由于离散元算法机理导致数值计算误差分别增至28.64%和29.45%。从而，采用基于离散元算法的数值计算方法，在冰厚较薄情况下计算船体所受冰阻力能够达到较高的精度，而在碎冰厚度较厚时该方法的误差偏大。