基于SPH-FEM 算法的浮冰-水-船耦合作用数值模拟方法

毕东魁，王翔宇，孙树政，陈睿童，郁 峰，尤婷婷

(1. 中国船级社，山东 威海 264200；2. 哈尔滨工程大学船舶工程学院，黑龙江 哈尔滨 150000；3. 烟台哈尔滨工程大学研究院，山东 烟台 264000；4. 中国船舶集团有限公司，上海 200011)

0 引 言

我国处于北半球，船运走北极航线必须经过边缘冰带（MIZ），这个区域处于极地和温带气候系统的交汇处。该区域浮冰的大小从数米到数百米不等[1]。这些充满随机性的浮冰会对船舶的总体结构与人员生命安全造成严重威胁。所以，船舶在浮冰区域的冰阻力预测成为船舶安全保障的关键。

对于船舶浮冰阻力的研究，有实船试验法、水池模型船试验法、计算机仿真模拟等，综合对比水池试验法在耗费人力物力、试验周期、结果精确性等方面有较高优势，也是现代学者研究冰载荷的一种较为认可的方式。陈昭炀[2]在拖曳水池中进行了模型冰船阻力试验，国内首次应用试验方法研究浮冰的形状对冰阻力的影响，观察浮冰的堆积和贴体情况，分析工况、浮冰运动、阻力之间的影响规律。Jeong[3]在冰池中建立了一条碎冰航道，并根据船舶在该碎冰航道中航行时拖曳力与推进阻力的关系，推导出带有螺旋桨模型船的性能。随着计算机的发展，可以做到对水池试验进行准确仿真，不仅可以保证数据精准，而且有容易控制变量、灵活度高等优势。徐栋梁[4]在研究船舶舷侧与浮冰的碰撞过程中考虑了温度场对船体钢材料的影响，总结温度、速度、位置等因素对加强筋结构碰撞的影响，根据仿真结果拟合出相对应的经验公式并验证。

本文应用光滑粒子流体动力学方法(smoothed particle hydrodynamics, SPH)与有限单元法(finite element method，FEM)相耦合模拟冰水池中船与浮冰碰撞。SPH 模拟流体具有自适应的优势，本质上仍然是拉格朗日方法，但在处理高速碰撞中不会出现传统FEM 中网格严重变形导致结果不精确的问题。SPH 粒子在耦合中被视为特殊的节点，控制参数为节点编号、质量以及空间位置，粒子之间自行耦合，与有限单元点面接触耦合，可以较为准确做到冰-水-船耦合。

1 数值模型

SPH 的基本原理为拉格朗日法。将连续的海水与浮冰离散成诸多粒子质点，这些粒子拥有独立的质量、位置坐标、速度、加速度等物理量。随后求解整体的动力学方程，跟随每个质点的各项物理参数随时间的变化情况，最终综合所有质点得到整个流场与冰体的物理情况。

SPH 方程构造首先采用核函数逼近的方法，将描述场的函数转化成积分表达式，然后粒子近似，实现了对核函数近似积分表达式的粒子离散。

粒子近似式方程表达为：

式中，W为核函数。

核函数W由函 数 θ定义：

式中：d为空间维数；h为光滑长度。光滑长度随时间空间而变化。

SPH 方程最常用的光滑内核是三次B-样条函数，由下式定义：

式中：C为标准化常数，由空间维数决定。

边界采用刚性墙与虚粒子2 种方法结合处理。方法1：用有限或无限大小定义1 个平面刚性壁，保证内部粒子无法穿过壁面。方法2：创建虚粒子来定义对称面，虚粒子为靠近边界两倍初始光滑长度范围内实粒子的镜像，与对称实粒子具有相同的质量、速度等物理量。不仅提高了计算效率，并且仿真结果的精度得到了保证。

接触算法采用罚函数进行求解[5]，粒子的接触力分解为法向接触力fn和切向接触力，计算公式为：

2 计算模型

2.1 LNG 船体模型

本文研究一艘航行于北极航道的LNG 船，破冰能力为ARC-7，自主破冰厚度可达1.7 m，主要参数如表1所示，仿真模型缩尺比例为1:100, 材料设为刚体。忽略船体与浮冰碰撞时的变形，简化船体结构，保留船体外壳与内部主要骨架，用添加质量点的方式平衡船体的质量、重心位置和转动惯量。

表1 船体主要参数Tab. 1 Main parameters of the hull

2.2 浮冰域

由于北极航道环境复杂多变，浮冰的大小与形状也不尽相同，忽略不同大小与形状带来的结果偏差，只保证浮冰的密集度。将海面漂浮的浮冰统一设置成大小相同的长方体，为模拟真实环境，长方体的排布为随机布置。浮冰为SPH 建立，粒子直径为0.01 m，材料的本构模型选择各向同性弹塑性断裂模型，失效准则满足VON-MISES 屈服准则，以最大塑性应变模式作为材料的破坏模式，材料的分离模式为恒定最小压力模式[6]，具体参数见表2。

表2 冰材料模型参数Tab. 2 Ice material model parameters

由于船模宽度相较于拖曳水池宽度来说为小值，适当缩小浮冰域宽度对实验结果影响不大，用刚性墙模拟冰水池试验中浮式围栏的影响，浮冰域设置为长9 m，宽2 m。

海水模型同样使用SPH 建立，粒子半径为0.02 m。水的状态方程采用Gruneisen 状态方程，压力方程为[7]：

式中：C1为 水中声音的传播速度； α为 对系数 γ0的一阶修正；S1-S3为 μs-μp曲线斜率无量纲系数；E为初始材料内能； μ为体积变化率。

图2 为船在水中静置10 s 的浮力随时间变化的曲线图。前2 s 船与水开始接触，浮力数值发生较大变化，后8 s 趋于稳定，浮力稳定在1 275 N 附近，而船缩尺后满载排水量为127 kg，可以认为水的浮力模拟结果较为准确。

图2 浮力-时间曲线Fig. 2 Buoyancy-time curve

2.3 浮冰工况与结果对比

北极地区冬季最低温可至-50℃以下，夏季平均气温在-10℃至10℃之间，某些地区最高温度可达30℃以上，这就导致同一条航道中浮冰情况会有很大差别。综合实际情况并缩尺选取浮冰厚度均为0.01 m，密集度为80%，60%，40%，船舶释放六自由度，仿真时长9 s，前1 s 船舶静置在水中达到正浮，后8 s 以0.448 m/s2沿船长方向直线匀速运动。水域四周与底部设置刚性墙防止粒子穿透，船体纵剖面设置虚粒子边界。得出的阻力值与Du Brovin 在1950-1955 年通过4 次船模水池试验得出的结果推导出针对船模试验冰阻力预报的经验公式[8]对比。计算船模试验冰阻力Rice的经验公式为：

式中：Rice为船舶所受冰阻力；p1，p2为经验修正系数，与浮冰密集度和浮冰域宽度与船宽的比值有关；Fn为 傅汝德数；n为 功率系数。参数A与 φ的计算公式为：

式中：B为船型宽；L为船总长；r为浮冰密集度；hice为浮冰厚度； ρice为浮冰密度；f0为船体与浮冰的摩擦因数； αH为船舶前体棱形系数； α0为水线面入射角的1/2。

3 计算结果分析

3.1 浮冰应变云图分析

各时间浮冰应变云图如图3～图5 所示。

图3 40%浮冰密集度的应变云图Fig. 3 Strain cloud map of 40% ice floe density

图4 60%浮冰密集度的应变云图Fig. 4 Strain cloud map of 60% ice floe density

图5 80%浮冰密集度的应变云图Fig. 5 Strain cloud map of 80% ice floe density

从时间云图可以看到，船首先与浮冰碰撞，接触部位的浮冰应变达到0.001 失效应变后发生破碎，随后被船首挤压到两边位置，对周围浮冰产生挤压碰撞。碰撞主要发生在船体的首部与肩部，平行中体部位与尾部附近的浮冰应变几乎没有变化。但由于波浪和浮冰之间碰撞浮冰会向船首方向运动，密集度越低的浮冰之间运动的相对距离会远，密集度高的浮冰会聚在一起，整体向前。当船全部进入浮冰区域后，船尾被碰撞的碎冰会重新聚合。

3.2 浮冰运动状态分析

如图6 所示，在碰撞时浮冰的运动状态分为3 种情况。

图6 浮冰不同运动状态Fig. 6 Floating ice in different motion states

在仿真过程中可以明显看到模型浮冰有推积、贴体、平移等现象。浮冰在与船首接触时发生破碎且周围浮冰的挤压，导致部分碎冰无法从两边沿着船身向船尾方向移动，于是这些碎冰会堆积到船首的位置。船首的型线会导致在碰撞浮冰时，浮冰不仅会受到船长方向的碰撞力，也会受到来自船首向下的挤压力。这些挤压力将浮冰挤压成大小不一的碎冰，碎冰受到水中浮力作用紧紧贴附在水下船壳表面。平移现象在密集度小的工况中更为明显。因为浮冰在船头的碰撞与船身的摩擦之下向前运动，由于前方水域较为开阔，无浮冰阻拦，所以仅收到水阻力的情况下向前漂浮一段距离后静止。因为船速相对较低，所以整体工况中反转现象并不明显。这些浮冰的不同运动情况是使冰阻力增加的一个重要因素。

3.3 不同密集度下冰力随时间变化曲线

船体分别在40%，60%，80%密集度浮冰区域航行时受到的冰阻力随时间变化曲线如图7～图9 所示。

图7 40%浮冰密集度冰力-时间曲线Fig. 7 40% Ice density ice force time curve

图8 60%浮冰密集度冰力-时间曲线Fig. 8 60% Ice density ice force time curve

图9 80%浮冰密集度冰力-时间曲线Fig. 9 80% Ice density ice force time curve

实线为冰阻力-时间曲线，The mean 划线为4-6 s 船首进入浮冰区后数据平稳段阻力的平均值。Experience 划线为DuBrovin 经验公式值。各项详细数值如表3 所示。

表3 各项冰力数值表Tab. 3 Table of various ice force values

各浮冰密集度冰力-时间曲线具有明显的动态非线性特征，趋势分为3 部分，第1 部分0-2 s 时船开始向前运动，与浮冰未接触，冰力为0。在2-4 s 中船首部与浮冰接触，冰力开始增大。4-10 s 中船身与浮冰接触，直到整个船身完全进入浮冰区域。在这个阶段冰力大浮动波动，是因为船与浮冰接触时冰力增大。浮冰收到碰撞力向力的相反方向移动冰力卸载，浮冰又受到水阻力与其他浮冰的影响速度减小与船体发生二次碰撞冰力增大，直至破碎或沿着船身移动。但浮动范围变化很小，说明船身对冰力大小的影响较小。

冰力极值方面密集度为40%时峰值点较多且出现时间较晚，说明船舶航行中波浪对于浮冰的运动起重要作用。

由图10 可知，随着浮冰密集度的增加，浮冰阻力均值与经验公式值均程增加趋势，近似呈线性关系。仿真数据比经验公式在60%，80% 的工况中数值偏大，是因为浮冰随着密集度增加会出现破碎、堆积、贴体等不确定因素导致的阻力增加。

图10 浮冰阻力均值和经验公式值与浮冰密集度的散点图Fig. 10 Scatter plot of mean and empirical formula values of floating ice resistance and floating ice density

4 结 语

本文基于SPH-FEM 耦合算法建立模型船在水中与不同密集度浮冰域碰撞的仿真场景并进行模拟，分析船体航行冰阻力变化情况，主要结论如下：

1）基于FEM-SPH 耦合算法对LNG 船在水中与浮冰碰撞过程有很好的模拟效果。

2）在LNG 船航行过程中，浮冰与船体接触被破坏区域主要在船首与肩部，船身处浮冰应变变化不大。浮冰运动状态主要有推积、贴体和平移，密集度小的工况中推积与贴体较少，浮冰平移较多。

3）冰力-时间曲线可反映船体与浮冰碰撞全过程。浮冰密集度小的工况中船体所受平均冰力较小，与DuBrovin 经验公式对比，整体变化趋势一致，在高密集度工况下仿真结果偏大。