不同船艏类型破冰模式差异的数值预报分析

2021-09-02叶礼裕

船舶力学 2021年8期

张媛，王超，叶礼裕，杨波

（哈尔滨工程大学，哈尔滨 150001）

0 引言

极地日益增长的航运活动、科学考察工作、资源开采需求以及军事领域的应用价值极大提升了极地结构物的需求量，并对其极地航行船舶的性能提出了更高的要求。破冰船是极地航行必需的开辟航道、引领冰区航行、保证极地船舶安全航行的重要交通工具[1-2]。破冰能力是衡量破冰船性能的关键标准之一，而破冰船的船艏形状特征参数在很大程度上决定了破冰效率、破冰模式，以及破冰第二阶段的冰翻转过程和浸没过程的碎冰块运动轨迹、运动模式和冰的浸没阻力。为了深度掌握船艏型线特征对破冰船破冰能力的影响，探索各关键船艏特征参数与破冰过程的相关性机理，以帮助设计出型线更优、性能更佳的破冰船艏，提高国家极地航行事业的竞争力，开展破冰船艏形状特征的破冰模式差异预报分析研究是十分必要的。

破冰船艏部形状可通过5 个特征参数来表达：外飘角、水线进流角、纵剖线倾角、首柱倾角、艏部长度[3]。其中，外飘角主要影响破冰船的破冰效率及破碎冰块的下沉效率，水线进流角的影响主要体现在清除船艏前端及其两侧形成的堆积冰的速率，即清冰效率，纵剖线倾角和首柱倾角也会一定程度地影响到破冰和冰块下沉效率。破冰船艏的设计主要是围绕上述特征开展的，目前常见的破冰船艏主要有7种：具有平行纵剖线的直线艏、凹形艏柱艏（怀特艏）、高外飘角艏（米勒艏）、具有车刀的勺形艏、具有脊骨的半勺形艏、平板系艏和狄森-瓦斯艏[4]。其中，前三种破冰艏均保持了平滑的船体，敞水域阻力性能较好，通常被称为传统型破冰艏；而后四种因艏部型线非平滑过渡，通常被称为非传统型破冰船艏。

数值预报方法被认为是研究冰-船接触过程中破冰模式以及破冰载荷的有效手段。在平整冰、冰脊和冰山等冰场中，常采用有限元方法（FEM）模拟局部冰载荷和全部冰载荷[5-6]。然而，在描述冰-船相互作用过程中的破冰特性时，由于离散元方法（DEM）既可以在微观尺度上描述冰结构的离散性，又可在宏观尺度上模拟船-冰相互作用过程中的冰破裂过程，DEM 表现出了其在冰荷载预报方面的固有优势[7-10]。文献[11]采用SPH方法模拟了层冰工况下，冰-船接触过程中的冰破坏过程并预报了破冰船的破冰阻力。近场动力学方法是另外一种具有模拟材料断裂特性优势的无网格方法，现已成功应用于冰的模拟预报中[12-14]。

本文的主要目的是应用数值预报的手段来分析传统型破冰船艏和非传统型破冰船艏的破冰模式差异。为了更好地模拟破冰模式，文中采用了一种新型非局部无网格方法——常规状态近场动力学理论（后文简称OSB-PD）来建立冰的弹脆性本构模型；建立了一种快速接触检测算法（后文简称FCDA）来处理基于粒子方法的固固耦合接触判断问题；嵌入了消息传递接口（MPI）的并行技术来提高无网格方法的计算效率；在FORTRAN 语言环境下编写了破冰船破冰过程的数值预报程序；预报了两种典型船艏形状破冰船的破冰过程；最后总结分析了两种船艏的破冰模式和破冰阻力的差异。

1 冰-船作用数值模型

1.1 冰的本构模型

文中冰的本构模型是基于OSB-PD 理论建立的各向同性、均质分布的标准弹脆性材料。在OSBPD 理论中，冰材料离散为无限个材料粒子点，通过关注每个粒子点的运动信息和物理信息，进而建立整个冰体的变形和破坏的数学模型。由于PD方法是非局部方法，因此每个材料点与其相邻一定范围Hx内的粒子存在相互作用关系，超出这个范围的粒子与该粒子不存在任何作用关系。在笛卡尔坐标系下，粒子空间位置为x，其占据一定的空间体积Vx，其密度用ρ(x)表示，邻域粒子的坐标为x'。当变形发生时，粒子发生位移u，且具有新的坐标位置y，同样地，邻域粒子的位置发生变化为y'。粒子间的相互作用用力密度T来表示，两个粒子间的力密度方向相反，大小不同，方向分别指向对方粒子。由此可以得知，两个粒子间的作用力是两个不同的力密度，分别为[x,t]和[x',t]。OSB-PD 方法的最终控制方程如式（1）所示[15]：

通过应变能密度和力密度关系与经典介质力学中对应关系的对照可推导出OSB-PD 方法的相关常量参数，代入公式（1）中，即可得到PD的详细积分表达式：

式中，a、d和b为PD 常量，Λ为辅助参数，δ为近场域的半径，θ和θ'分别为当前粒子点和当前作用邻域粒子点的体积膨胀，s为粒子间相互作用的伸长量，b为外部作用力。s和θ的表达式分别为

材料破坏的标准通过伸长量s表示，当s超过材料的极限伸长量时，则材料发生破坏，这种破坏过程是不可逆的，即粒子间变形伸长量超过极限之后，相互作用会永久消失。由此可以简单引入一个历史变形状态标量Ω来表示粒子间的作用关系：当Ω取值为0 代表作用关系消失；当Ω取值为1 时代表作用关系仍旧存在。

1.2 接触检测算法（FCDA）

在FCDA 方法中，船体离散为一系列足够表述船体表面形状的四边形面元，则冰粒子与船体的碰撞检测过程可以简化为空间中点和面的相对位置与距离判断的数学问题。检测步骤如下：

（1）冰材料离散为粒子形式，船体划分为接近平面的四边形的面元；

（2）在t时刻，任何一个粒子都有与曲面接触的可能性，为了实现高效的搜索，可以先排除该时刻完全不可能与曲面接触的粒子以减少粒子搜索量，因此建立一个包含整个不规则目标碰撞结构的规则矩形体，只有进入这个矩形体的粒子才有可能与目标面发生碰撞，由此排除了大量不可能发生接触的粒子，减少了不必要的粒子搜索过程；

（3）在t+1时刻冰粒子穿透矩形体和目标面，接下来只需关注进入矩形体的所有粒子；

（4）在这些穿透矩形体的粒子点里面找到每个粒子点接触或者穿透的唯一面元；

（5）假设粒子的坐标为(x0,y0,z0)，对于目标撞击体表面的所有四边形面元，可以找到其四个点在x方向上的最小值xmin和最大值xmax，在y方向上的最小值ymin和最大值ymax，以及在z方向上的最小值zmin和最大值zmax。若面元与点关系为xmin＜x0＜xmax，且zmin＜z0＜zmax，或者xmin＜x0＜xmax，且ymin＜y0＜ymax，则可认为这些面元可能与该物质点发生碰撞，该面元为可能与粒子碰撞的元素；

（6）在上一步找到的可能碰撞面元中，开始判断粒子是否穿透或者接触面元，此时接触检测过程已经简化为空间点面位置关系的数学问题，应用空间点面距离公式即可，建立面元的空间方程Ax+By+Cz+D=0，则最终判断粒子是否与面接触的准则如下：

冰粒子与船体的接触检测过程完成后，对进入船体的粒子重新分配位置，并更新粒子的运动信息和计算冰载荷，该过程可直接参考文献[15]第十章的处理过程，本文不再解释。

1.3 信息传递接口并行技术（MPI）

本文编译了基于MPI手段的OSB-PD 并行程序，可以在小内存的硬件上实现大数据量的运算，也可以为以后的超大数据量的并行扩展提供基础技术保障。

如图1（a）所示，模型冰层在长宽方向上远远大于厚度方向，因此针对冰层所在平面，将计算域平均分解为两个维度的9 个子域并行进程，所占用线程从0 编号至8。该平整冰层x方向上的粒子数量为nx，y方向上的粒子数量为ny，z方向的粒子数量为nz，划分在各个处理器后三个方向上的粒子数分别为、和，且=nz，由此每个处理器中的总粒子数为n=··nz。PD 方法中每个粒子与其邻域粒子存在作用关系，邻域的尺寸为δ=m·dx，dx为粒子间距，m为正整数，代表倍数关系。由此并行中每个处理器边界处的粒子会与相邻处理器中的粒子存在信息交互传递，要传递信息的粒子占据的宽度为m个dx，假设m=3，则可以设定不同相邻处理器之间的信息交换的计算域宽度为3dx，如图1（b）所示，则每个处理器最大粒子总容量为ntot=n+nz·2 +·nz·2 +nz·3·3·4。

图1 并行计算域的划分示意图Fig.1 Schematic diagram of the parallel computing domain partition

通过上述的并行策略和计算域的划分方法，使用MPI 常用指令，即可编译基于MPI 方法的PD 并行程序。

2 破冰过程数值预报分析

基于上述建立的数值计算方法，本章分别建立了传统型破冰船艏和非传统型破冰船艏的破冰船连续破冰的数值模型。为了探究不同船艏类型在破冰过程中的差别，本文预报了相同工况下两种破冰船船艏的破冰过程，并对比分析了破冰模式和破冰载荷。

2.1 两种船艏形状描述

传统型的破冰船艏是具有平行纵剖线的直线艏型，该船型起源于1950 年Souiet和Finland 的破冰船，呈尖瘦型的船艏形状，具有优秀的破冰能力，是发展至今仍旧被广泛应用和改进的破冰船艏型。非传统型破冰船艏选用简化的狄森-瓦斯型艏（Thyssen-Waas Bow），这种船艏和传统型船艏相比，具有明显的区别，该船型主要靠船艏最大宽度位置来造成冰面的剪切破坏进而依靠重力造成冰面的弯曲破坏，该船型具有极好的冰清除能力，可以开辟较宽且碎冰较少的破冰航道。两种船型的模型图如图2所示，详细的船型参数如表1所示。本文中两种船型的计算缩尺比均为25。

表1 两种船型的主尺度参数Tab.1 Principal dimensions of two kinds of icebreakers

图2 两种不同船艏形状破冰船模型图Fig.2 Icebreaker models with two different bow shapes

2.2 计算输入

冰参数、离散信息和工况设定参见表2，该输入信息对应于模型的缩尺比25。例如冰的实尺度厚度为1.0 m，则计算中输入为0.04 m。参量的换算关系参考ITTC冰试验的换算规则[16]。

表2 计算模型参数Tab.2 Calculation parameters of the model

2.3 计算结果

2.3.1 数值方法验证

本小节首先对比了数值计算结果和试验结果的破冰模式、数值计算结果的破冰载荷与试验[17]和经验公式[18]的破冰载荷来验证数值方法的可靠性，该实验中的破冰船采用本文中传统型船艏破冰船。破冰模式的结果对比如图3 所示。为了和经验公式的结果做对比，破冰载荷依据ITTC 中冰区船舶换算公式换算为实尺度的冰载荷数据，对比图如图4所示。

图3 破冰模式验证示意图（传统破冰船艏；实船4节航速）Fig.3 Icebreaking pattern verification comparison snapshot(Conventional bow with prototype;speed:4 kn)

图4 破冰阻力验证对比图（传统破冰船艏；实船4 kn航速）Fig.4 Icebreaking resistance verification comparison(Conventional bow with prototype;speed:4 kn)

从破冰模式来看，环形裂纹的生成和扩展引发了冰层的弯曲破坏模式，整个破冰过程的破坏模式主要都体现为弯曲破坏模式。每一次冰层破坏模型的循环基本为：首先产生沿船长方向的环形裂纹并向船艏两侧不断扩展，该类型裂纹总是从船艏的半船宽处开始成型；随着船舶的前进，环形裂纹向部和船艉扩展，同时在船艏半船宽处产生放射状的裂纹；紧接着，基本平行于船艏侧边的二次环状裂纹产生并且扩展，伴随着局部区域的冰挤压破碎；船艏产生较短的环状裂纹。数值模拟结果很好地捕捉到了上述在试验中观测到的破冰现象。

从破冰阻力均值来看，数值计算破冰力的均值为0.969 MN，试验破冰力的均值为1.124 MN，经验公式计算得到的均值力为1.117 5 MN。由此可以得到结论：数值计算结果与试验结果和经验公式的结果具有很高的吻合度，本文建立的数值计算方法能够高效地模拟破冰过程和准确地获取破冰阻力。

2.3.2 两种船艏的破冰模式

本小节计算了两种破冰船艏在实船航速为3 kn 的工况下的破冰模式和破冰载荷，分别如图5和图6所示。

从图5 可以发现，两种不同船艏的破冰船破冰过程表现出明显不同的破冰模式。在裂纹初始阶段，传统尖瘦型的破冰船艏首先对冰层造成局部挤压，形成船艏形状的破冰缺口，进而在船艏两侧半宽处产生沿船长方向的环状裂纹，向艉部和部演变发展，此时的破冰模式体现为船艏轮廓处的挤压破坏和平行于船艏外侧轮廓线的弯曲破坏，船艏两侧形成少许碎裂的冰块；而非传统型船艏在冰船接触位置发生极小区域的挤压破坏，形成轻微的类似于船艏形状的弧状冰层缺口，当破冰船进一步穿过冰层时，由于船艏没有明显的艏柱冲破冰层，而是平缓的过渡型船艏压向冰面，形成了典型的放射状裂纹和穿越放射裂纹端部的圆弧状裂纹。在裂纹扩展演变阶段，传统型船艏破冰模式和冰层破坏的模式循环过程与2.3.1小节介绍的相同；而非传统型船艏的破冰模式简单地由放射状裂纹和穿越放射裂纹的圆弧状裂纹组成，冰层破坏模式的循环过程为放射状裂纹和圆弧状裂纹有规律地层层扩展，且每次破坏循环周期内的圆弧状裂纹半径都比上一个循环大。总体来看，两种破冰模式中，弯曲破坏主导了整个破冰过程，冰层破坏的成型形状基本与船艏轮廓相一致。

图5 两种类型船艏的破冰模式对照Fig.5 Comparison of icebreaking patterns of two bow types

从实尺度破冰航道来看，传统尖瘦型破冰船的型宽为23.0 m，开辟的航道宽度为28.56 m。非传统型船的型宽为20.0 m，破冰宽度为27.3 m。两种船型的破冰航道宽度与船宽的关系分别为1.25B和1.36B。因此，对比这两种破冰船艏可发现，非传统型船艏破开航道宽度的能力更显著。另外，传统型船艏的航道周边有明显的冰尖形状，而非传统型的破冰航道边缘更平滑。最后还可以发现，传统尖瘦型船艏的破冰过程形成的碎冰块普遍偏大，会对后续船舶的航道航行造成一定的影响，非传统平缓船艏的破冰块的形状和大小分布较均匀，且尺寸偏小。

2.3.3 两种船艏的破冰载荷

两种破冰船艏破冰力和经验公式的对比如图6 所示。从破冰载荷来看，两种船艏的载荷趋势和破冰模式相对应，传统型艏的冰载荷具有连续性的特征，这是由于传统型船艏破冰过程没有明确的周期性分界线，往往下一个破冰循环过程和本次破冰循环过程的裂纹扩展同时进行；非传统型艏的冰载荷具有明显的周期循环特征，每个周期的冰载荷呈现先增加后减小的趋势，增加过程对应于弧状裂纹的产生和扩展，减小的过程对应于放射状裂纹的扩展进而该区域冰全部破碎为小冰块。