0 引　言

船舶在向极地航行的过程中，经过了开敞水域、极地冰缘区至中央浮冰区等过程，此过程冰块尺寸逐渐增加。国内外学者研究船舶阻力性能的2种方法。一是基于CFD数值模拟，运用势流或粘流方法；二是基于模型试验法。张金龙等^[1]采用离散元法，研究了船舶在开敞水域、沿岸壁碎冰水域航行时的冰载荷特性，并总结了船舶沿岸壁航行，冰阻力的变化规律。Dogrul等^[2]分别对KCS和KVLCC2数值模拟，对比了不同尺度下水动力及船体周围流体特性。Kazeroni等^[3]采用模型试验，研究了高速模型船阻力性能，发现低速航行主要是摩擦阻力。Marco等^[4]采用模型试验法、RANS和LES法分析了相同流动条件下，发现LES法更适合模拟船体周围的水流现象。Corlett等^[5]以某破冰船为船模，第一次采用非冻结模型冰石蜡完成了冰区船阻力模型试验。Luo等^[6]采用CFD-DEM方法计算了船-冰-流间的相互作用，对碎冰航道中船舶所受冰载荷进行了预报。McGovern等^[7]采用模型试验法，研究了波浪对浮冰作用的影响，得出波长对浮冰运动的影响最大，而傅汝德数、浮冰尺寸和粗糙度的影响较小。Cundall等^[8]引入了离散元法，并应用于冰力学研究。Løset等^[9]使用离散元模拟了受外力作用碎冰场的响应特性，采用线性粘弹性模型研究了船冰碰撞。Huang等^[10]采用离散元法，基于STAR-CCM+，研究了不同航速、密集度、海冰尺寸下的冰阻力。李正等^[11]基于STAR-CCM+，研究了碎冰翻转过程中的船冰摩擦系数、碎冰厚度、航速对冰载荷的影响。Xue等^[12]基于CFD-DEM方法，考虑了流体速度和船舶波浪的变化，并预测了破冰中的船舶阻力。Tang^[13]采用耦合方法对船舶所受冰载荷进行了预测，发现船首浮冰的碰撞、堆积和挤压是造成船冰阻力的主要因素。Song等^[14]使用Ls-Dyna中流固相互作用技术进行验证，确定了球体和一个矩形块的附加质量系数。

综上所述，发现开敞水域中船舶阻力研究全面且细致，冰载荷特性的研究相对较少。但国内外关于冰载荷的影响分析较少，考虑到极地船冰作用的复杂性，不同影响因素对极地船舶航行安全有很大的影响，本文考虑碎冰的影响，对船冰耦合水动力特性进行讨论分析。

1 数学模型 1.1 流体控制方程

在CFD研究领域，牛顿第二定律，质量及能量守恒定律构成了流体的三大定律。任何流动问题都要满足质量守恒定律和动量守恒定律。

$ \frac{{\partial \rho }}{{\partial t}} + \frac{{\partial \left( {\rho {{\boldsymbol{u}}_i}} \right)}}{{\partial {x_i}}} = 0 ，$

(1)

$ \frac{{\partial \left( {\rho {{\boldsymbol{u}}_i}} \right)}}{{\partial t}} + {\rm{div}} \left( {\rho {\boldsymbol{u}}{{\boldsymbol{u}}_i}} \right) = - \frac{{\partial p}}{{\partial {x_i}}} + \frac{{\partial {{\boldsymbol{\tau}} _{xi}}}}{{\partial x}} + \frac{{\partial {{\boldsymbol{\tau}} _{yi}}}}{{\partial y}} + \frac{{\partial {{\boldsymbol{\tau}} _{zi}}}}{{\partial z}} + {{\boldsymbol{F}}_i}。$

(2)

式中：$t$为时间；$\rho $为流体密度；$p$为压强；${\boldsymbol{u}}$为速度矢量；${{\boldsymbol{u}}_i}$为速度分量；${x_i}$为速度分量位移；${{\boldsymbol{F}}_i}$为质量力；${{\boldsymbol{\tau}} _{xi}}$，${{\boldsymbol{\tau}} _{yi}}$，${{\boldsymbol{\tau}} _{zi}}$为黏性力${\boldsymbol{\tau}} $的分量。

1.2 离散元接触模型

碎冰区航行的船舶，存在碎冰与船舶的碰撞。选择弹簧接触模型进行碰撞力计算，如图1所示。

无滑移接触模型Hertz-Mindlin，用来模拟粒子间接触。A和B两球体间作用力为：

$ {{\boldsymbol{F}}_{{\text{contact }}}} = {{\boldsymbol{F}}_n} + {{\boldsymbol{F}}_t} ，$

(3)

$ {{\boldsymbol{F}}_n} = - {K_n}{d_n} - {N_n}{v_n} 。$

(4)

式中：$ {K_n} $为法向弹簧刚度；$ {d_n} $为离散单元接触点处的法向重叠；$ {N_n} $为法向弹簧阻尼；$ {v_n} $为接触点物体表面速度的法向分量。

切向力表达式为：

$ {F_t} = \left\{ {\begin{array}{*{20}{l}} { - {K_t}{d_t} - {N_t}{v_t},}\;{\left| {{K_t}{d_t}} \right| \lt \left| {{K_n}{d_n}} \right|{C_{fs}}}，\\ {\left| {{K_n}{d_n}} \right|{C_{fs}}{d_t}/\left| {{d_t}} \right|,}\;{\left| {{K_t}{d_t}} \right| \geqslant \left| {{K_n}{d_n}} \right|{C_{fs}}} 。\end{array}} \right. $

(5)

式中：$ {K_n} $、$ {K_t} $分别为法向弹簧刚度和切向弹簧刚度；$ {N_n} $、$ {N_t} $分别为法向弹簧阻尼和切向弹簧阻尼；$ {C_{fs}} $为静摩擦系数；$ {d_n} $、$ {d_t} $为法向和切向相对位移；$ {v_n} $和$ {v_t} $分别为接触点处碎冰单元间法向和切向速度。

2 数值仿真模型设置 2.1 KCS船体模型

选取ITTC公布的标准KCS船模，该船拥有完善的实验数据。KCS的主要尺寸见表1，KCS模型见图2。

2.2 计算域及边界条件

计算域的设置满足ITTC标准，三维计算域长宽高分别为4.5L×4L×3L（L为船长）。计算域及边界条件如图3所示。

2.3 碎冰模型及海冰参数

表2为海冰参数。本文选择圆柱状碎冰，采用多相间的相互作用。

2.4 网格设置

数值计算时，采用棱柱层网格，对船体附近和水线面进行加密，图4为水线面、船体表面及棱柱层网格。

3 数值计算模拟 3.1 开敞水域船舶阻力精度预报

以航速v=2.196 m/s为例，靠近壁面第一层网格满足$ 30\leqslant {y}^{+}\leqslant 300 $，按照$ \sqrt[4]{2}$的比例进行阻力验证，图5为不同网格下的船舶阻力计算结果。

计算发现，7套网格计算结果与试验值误差在0.5%～2.6%之间，292万网格下数值模拟与试验值相差最小。

3.2 冰水阻力网格收敛性验证

考虑网格尺度对阻力及冰载荷计算结果的影响。按照$1:\sqrt 2 :2$的尺度比由密到疏建立了3套网格，参数如表3所示。

碎冰水域中，选取0.915、1.281、1.647 m/s等3组航速，碎冰厚度为0.02 m，半径为0.12 m，密集度为80%。不同航速下，水阻力数值结果与试验值对比、3套网格下冰载荷如图6所示。

图6(a)中水阻力与试验值误差在10%以内，图6(b)～图6(d)对应3套网格下的总冰载荷，发现总冰载荷均值相差较小，水阻力及冰载荷的计算误差在M₂套网格下最小，因此选择M₂网格。M₂网格下，不同航速的冰载荷均值分别为7.01、9.51、12.04 N。

3.3 不同水域兴波对比

选择2个航速1.647、1.922 m/s以船体的中纵剖面为参考，研究了船模右舷距离中纵剖面B/2、L/2、L 3处，不同航速下船模波高曲线见图7，对比并得出了无冰和密集度为80%时碎冰条件下波高随航速变化的规律。计算发现，船模首尾出现明显的兴波。由于船模首尾存在水流流速停滞，伯努利方程中流速减小，水压增大形成波峰。产生低水压区是由于船肩部流速增大，形成波谷。

开敞水域中，随着航速的增大，船模周围兴波作用增强，曲线变化趋于一致，在$Y = B/2$处，船模位置向两侧，兴波产生的波高不断减小。从船尾向船首方向，波峰与波谷交替出现。且船模尾部的波高整体大于艏部波高。在$Y = B/2$处，由于水面处碎冰的缓冲作用，碎冰水域处波面起伏较开敞水域减小。在$Y = L$处，此处的波面起伏很小，在0附近处振荡。

4 结果讨论与分析

基于以上条件，本文研究航速、碎冰厚度、碎冰半径、密集度等变化参数下，冰载荷及冲量大小的变化规律。

冰阻力产生的主要来源为船-冰之间的碰撞、挤压和堆积。船首处的碎冰由于船冰撞击被推离至船体的两侧。在流的作用下，部分碎冰沿着船体两舷侧发生滑动，船后方形成了约船体宽度的无碎冰水域，自由液面处形成了清晰的开尔文波。

4.1 航速对冰水阻力及冲量影响

选取密集度为$80\% $的碎冰区，碎冰粒子厚度为0.04 m，直径为0.24 m，此时满足冰厚相似定律，选取的航速分别为0.915、1.281、1.647 m/s，对应的弗汝德数为0.108、0.152、0.195。计算结果如图8所示。

计算发现，船体纵向冰载荷与船舶运动方向相反，船体冰载荷在y、z方向上始终在0附近处振荡，冰载荷不断地加载与卸载。随着航速的增加，冰载荷在3个方向处的峰值点不断增大。总阻力中，船体所受冰载荷和水阻力占比不同。图9选择对应航速分别为0.915、1.647 m/s下的碎冰粒子速度云图和冲量大小。

船舶航速的增加，冲量峰值点呈增大的趋势。航速较小时，开敞水域及碎冰水域船舶横向水阻力吻合较好，当航速大于1.281 m/s时，两环境下的纵向水阻力差值逐渐增大。

4.2 碎冰厚度对冰载荷影响

为了研究碎冰厚度对冰载荷大小的影响，选择航速为1.281 m/s和密集度为80%的工况，船体所受的冰载荷时历曲线如图10所示。表4为数值模拟结果。

碎冰厚度由0.01 m增大至0.04 m时，船体纵向、横向及垂向冰载荷呈非线性增加趋势，碎冰厚度的增加使碎冰与船体的接触面积增大，且碎冰粒子质量增加，碎冰与船体间的碰撞力和摩擦阻力增大，碎冰粒子动能增大。图11为不同厚度下波面流场及碎冰粒子速度云图。

4.3 碎冰半径对冰水阻力及冲量影响

保持航速1.281 m/s、密集度为80%，碎冰厚度0.04 m不变的情况下，碎冰粒子的半径分别为0.05、0.08、0.10、0.12 m，表5为不同碎冰半径下的数值计算结果。

半径由0.05 m增大至0.12 m，冰载荷纵向均值增大约3倍，横向与垂向冰载荷变化较小，但垂向冰载荷呈下降趋势。纵向冰载荷仍为主要占比，垂向冰载荷次之，横向冰载荷占比最小。总冰载荷随碎冰半径的增大而增大，船冰间碰撞更为充分。图12为不同碎冰半径下船体受冰载荷时历曲线。

当碎冰半径由0.05 m增大至0.12 m时，冰载荷均值由10.36 N增加至25.12 N，增加了约2.5倍。碎冰半径的增大使得冰与船舶接触面积增大。冲量量度着动量变化的多少，随着航速的增大，冲量峰值点明显增大。因此极地航行时，考虑到船舶的结构损伤，需将航速限制在合理的范围内。图13为水线面碎冰半径0.05 m和0.12 m碎冰总冲量时历曲线及船舶周围流场分布。

碎冰半径的增大，船舶与碎冰间总冲量增大。且碎冰间及船与碎冰之间能量随碎冰半径的增大而增大。由于碎冰半径的增大使得碎冰质量增加，碎冰之间、船与碎冰间的动能增大。

4.4 碎冰密集度对冰载荷影响

密集度随着温度而变化，冰厚为0.02 m，碎冰半径为0.12 m，航速为1.281 m/s时，密集度分别为20%、40%、60%、80%。表6为不同密集度冰载荷的计算结果。密集度的增大使得碎冰的聚集更为明显。因此船体3个方向上冰载荷增加更显著。表6为不同密集度下冰载荷的数值计算结果。

图14为不同密集度下船舶总冰载荷及冲量时历曲线。

随着密集度的增大，碎冰粒子能量增加，船冰碰撞概率增大，碰撞力也增大。在不同密集度下，船舶水阻力在28.21 N左右，基本保持不变。因此密集度并不是影响船舶水阻力的决定因素，但它对冰载荷大小影响较大。随着密集度的增加，碰撞力在短时间内变化更为剧烈，船冰碰撞产生的冲量峰值点及变化区间范围增大，船舶与碎冰粒子碰撞产生的动能明显增大。

5 结　语

本文采用离散元法，对碎冰水域KCS船舶进行水动力数值预报，研究了不同航速、碎冰大小、密集度对冰载荷及冲量的影响及船冰碰撞时船舶周围的流场分布。将开敞水域与碎冰水域的波面起伏情况进行对比，研究发现：

1）船舶阻力数值模拟值与试验值吻合，构建的数值计算模型正确。船舶航行碎冰区时，由于水面碎冰的缓冲作用，碎冰水域处波面起伏较开敞水域减小。沿船体型宽方向，距船体越远时，波面起伏逐渐减小。

2）船体纵向冰载荷始终与船舶运动方向相反，船体在横向及垂向冰载荷上的受力较复杂。随着航速的增加，动量变化增大，冲量明显呈增大趋势。且碎冰区水阻力大于开敞水域的水阻力。航速较小时，不同航速下横向水阻力时历曲线吻合度较高。

3）随着碎冰半径和厚度的增大，船冰碰撞更充分。冰载荷、冲量及碎冰能量均随碎冰厚度、碎冰半径的增大而增大，冲量峰值点及变化区间范围增大。

4）随着碎冰密集度的增加，船冰碰撞概率增大。碰撞力在短时间内变化更剧烈。船舶冲量、冰载荷及船舶与碎冰粒子碰撞产生的能量均增大。