破冰船在冰区航行时，船体周围流场受到冰面的干扰，其流场特性必然发生变化，故在求解船舶所受水动力时需额外考虑冰面干扰。水动力导数求解相对于船舶其他方面的研究较少，主要存在经验公式法、约束船模实验法、数值模拟方法等，其中约束船模试验是目前公认的求解水动力导数较为准确的方法。针对敞水水动力导数研究主要有：吴兴亚等^[1]基于民用打捞型船采用CFD软件STAR-CCM数值模拟敞水中船舶斜航运动，验证了CFD方法求解水动力导数的可行性。Ohmori等^[2]等采用有限体积法分别对Esso Osake船模作斜航运动的黏性流场和水动力进行了计算。Alessandrini等^[3]对系列60船模作定常斜航以及回转运动时的黏性流场和所受侧向力和纵向分布进行了计算。Carrica等^[4]采用RANS方程数值模拟船舶较快速度下周围流场分布并计算了部分水动力，后期又针对回转以及Z型操纵运动进行了相关的数值模拟。Miller^[5]使用了CFDSHIP-IOWA程序对斜航、横荡、艏摇等RANS方程进行了求解从而进行了相关数值模拟。Ismail等^[6]对裸船体模型在0°以及12°漂角进行了斜航数值模拟，得到了船舶水动力导数的预报采用2阶TVD格式更精结论。Cura-Hochbaum^[7]对KVLCC1的邮轮进行了PMM运动数值模拟，预报相应的水动力导数，并将其代入操纵模型中与自航试验对比，验证了RANS方法的可靠性。Simonsen等^[8]对KCS船舶做艏摇运动进行了数值模拟，对所受的力以及力矩的不确定进行了分析。Simonsen等^[9]基于求解器CFDShip-Iowa，数值模拟PMM运动试验，求解部分水动力导数值，数值模拟了KCS船舶的Z型操纵运动以及全回转运动。

Kim等^[10]应用有限元软件LS-DYNA模拟了碎冰航道阻力，并与冰水池实验结果进行了对比，验证了方法的准确性。Vroegrijk^[11]利用STAR-CCM+软件数值模拟了船舶在碎冰航道中的航行时相关性能。Lau等^[12]使用离散元手段针对海冰与船舶作用下破碎模型以及船舶阻力进行了研究。刘晨飞等^[13]运用CFD以及重叠网格技术模拟了斜航、纯横荡运动分别计算了月池封闭以及打开时船舶水动力导数，得到了水动力导数会因为月池存在有所增加。孟庆杰^[14]采用湍流模型中两方程K-W SST模型模拟了浅水工况下船体的吸底效应，不同漂角下水动力性能以及船体周围流场的变化。王建华^[15]采用自主开发求解器naoe-FoAM-SJTU数值模拟了船舶斜航运动，得到了不同漂角下的转艏力矩以及侧向以及相应的水动力导数。针对冰区相关的研究主要有：郭春雨等^[16]使用LS-DYNA软件模拟了船舶碎冰阻力以及船- 冰作用过程，得出数值模拟结果在较低速度条件下在定性上与试验结果比较符合。任奕舟等^[17]采用有限元方法建模对冰区船舶阻力进行了数值模拟，得到的冰阻力较为准确，适用于预报破冰船冰阻力。齐江辉等^[18]使用STAR-CCM+软件数值模拟了碎冰区的航行过程，分析了航速对航行阻力的影响以及碎冰在船体周围的运动状态，得出船舶航行阻力随着航速的增加有所提高，碎冰阻力由于存在波浪排开作用随着航速待提高缓慢增加的结论。涂勋程^[19]使用软件LS-Dyna在碎冰环境和平整环境下物探船的冰阻力预报及相关参数敏感性开展了数值模拟研究并与经验公式对比误差较小，验证了耦合方法的准确性。王超等^[20]运用离散元模型结合欧拉多相流，对船体与碎冰之间的相互作用关系进行探索，计算不同航速、不同碎冰密集度下船体的受力情况，从直观上解释碎冰阻力的变化原因。

碎冰区冰水耦合下船舶水动力导数预报的研究缺失，可查到的文献考虑冰水耦合也仅仅用于碎冰阻力方面未考虑耦合后对流体阻力本身的影响，大部分只针对于船舶直航时碎冰阻力以及水阻力独立研究。本文拟采用STAR-CCM+软件中DEM模块下双向耦合模式对碎冰区中某型船舶进行不同漂角下的斜航运动模拟，求取船舶所受的侧向力以及转艏力矩，进而求解出冰水耦合下的部分水动力导数，为今后考虑耦合下流体增阻后的水动力导数以及操纵性能预报方法研究提供理论基础和技术支撑。

1 基础理论和计算模型

本文主要采用CFD软件STAR-CCM+手段获得，需要满足连续性方程以及动量守恒方程。选用目前广泛使用的标准的k-ε模型处理标准的两方程中湍流问题。依靠STAR-CCM+中DEM模块下的双向耦合模式，进行冰块与水的动量、能量交换达到耦合作用，监测耦合后的船舶所受的水动力。由于双向耦合基础理论内容较多请详见STAR-CCM+软件帮助文档，此处不进行详述。

1.1 DEM基础理论

碎冰工况航行时存在船- 冰- 水耦合作用的过程，本文采用CFD软件STAR-CCM+中DEM模块模拟碎冰粒子，碎冰粒子由若干DEM基本颗粒组成。拉格朗日构架下，通过显式计算进行计算域中所有粒子的轨迹跟踪。CFD-DEM耦合计算包括冰粒子之间的相互作用以及离散相冰与流体相水之间的相互作的计算，离散相DEM冰粒子的运动通常包括平行运动和转动运动两种方式，需满足平动动量守恒以及旋转角动量守恒方程具体形式为：

$ {{m_i}\frac{{{\rm{d}}{v_i}}}{{{\rm{d}}t}} = \sum\limits_j {{F_{ij}}} + {F_g} + {F_{{\rm{fluid }}}}} $

(1)

$ {\frac{{\rm{d}}}{{{\rm{d}}t}}{I_i}{\omega _i} = \sum\limits_j {{T_{ij}}} } $

(2)

式中: m_i、v_i和w_i分别为离散冰项粒子编号i的质量、速度和角速度; I_i为编号i的粒子的转动惯量；R_i为粒子i的半径；F_g=m_ig为粒子i的重力；F_fluid为流体对粒子i的作用力；F_ij为粒子i与粒子j或壁面之间的碰撞力以及其他作用在粒子上的非接触力；T_ij为接触力矩，表示除粒子重力以外的接触力在粒子上产生的扭矩，T_ij=r×F_ij－μ_r|F_ij|ω_i/|ω_i|, 其中r为粒子重心到粒子接触点的矢量。

1.2 数学模型的建立

单一来说冰粒子主要属于固体力学研究范畴，流体力部分属于流体力学研究范畴，软件STAR-CCM+中双向耦合模块可实现二者耦合。DEM参数的设置具有开放性，可进行变参数研究。围绕船体建立数值水池，为模拟无限水域，速度入口距离船首部3.0 L，压力出口距离船尾部3.5 L。水池左右两侧距离船中3.3 L，水池顶部以及底部分别距离船舶2.0、3.5 L，如图 1所示。

1.3 几何模型建立

碎冰区冰块厚度范围在0.3~3.0 m^[19]，小型碎冰跨度在2.0 m左右，大型碎冰尺寸跨度为0.5~10.0 km^[19]。考虑模拟采用的碎冰粒子是用DEM复合颗粒填充形成的碎冰粒子，尺寸过大喷射碎冰粒子耗时多，综合考虑本次碎冰模型采用小型碎冰，长度、宽度、厚度分别为3.0、3.0、1.5 m，按缩尺比40进行建模，碎冰如图 2所示。碎冰的DEM粒子填充数为120。采用STAR-CCM+软件中探针进行喷射碎冰，其中探针的长度为6.0 m, 分辨率为38, 如图 3所示，喷射时颗粒流率为每秒喷射9个。某型冰区船模主尺度以及碎冰物理性质见表 1，其中表 1中碎冰的物理参数指的是依据实尺度冰属性与缩尺的关系换算得到的模型冰参数。船模几何图形如图 4所示。

1.4 网格划分

斜航工况漂角选择0°, 2°, 4°, 6°, 8°网格划分方法为各工况网格尺度不变，需要将船体分别旋转2°, 4°, 6°, 8°将船体与计算域进行布尔减运算，然后划分网格。模型网格划分选用切割体网格、棱柱层网格以及表面重构相结合的方式。总体网格基础尺寸为2.3%L，船体表面网格大小选取为基础尺寸的12.5%，船舶尾部进行网格加密，尺寸为基础尺寸的6.25%，为使船体周围的流场过渡缓和，设置棱柱层数6层，棱柱层延伸率1.2。船体网格图如图 5(a)所示，尾部加密图如图 5(b)所示。使用DEM模拟冰粒子时粒子求解是否收敛很大程度上取决于水线面处网格密度，考虑防止自由液面发散，分别对自由面进行x-y-z向加密，z向采用两层加密，网格尺度分别为基础尺寸的25.0%、12.5%。x-y方向网格尺度均为基础尺寸的80%。以斜航0°为例。计算域整体网格图如图 5(c)所示，z向加密图如图 5(d)所示。

2 计算方法的验证

本次数值计算方法的验证依托哈尔滨工程大学船模拖曳水池中进行的76K冰区加强散货船^[21]的碎冰实验。76 K船模网格设置以及DEM相关设置与本文某型破冰船数值模拟设置一致，通过验证76 K数值模拟的网格以及DEM设置的合理性从而验证本次数值模拟方法的准确性。本次验证选择一种碎冰密集度40%、一个航速点0.613 m/s，进行敞水船模阻力、碎冰阻力以及碎冰姿态验证。进而说明本次计算网格设置的合理性。阻力图如图 6所示，图 7为冰块与船艏相互作用的数值模拟与实验测量的对比图。图 8为船体航行后尾部区域碎冰分布图。静水阻力取图 6(a)后5 s均值为-5.157 N，静水阻力实验值为-5.160 N，可看出数值模拟结果与试验吻合度非常高。碎冰阻力值取图 6(b)25 s后稳定段均值为-4.950 N，实验值为-4.330 N，误差控制在15%以内，满足目前行业内冰阻力预报误差不超20%的要求。图 7中可看出76 K冰区加强散货船以及某型破冰船的数值仿真结果与76 K试验结果在船舶首部均出现了碎冰堆积现象。图 8中可以看出76 K冰区加强散货船以及某型破冰船的数值仿真结果与76 K试验结果在船舶尾部均出现了碎冰向船两侧移动形成的尾迹图。综上所述本次数值模拟具有一定的准确性，可以较为真实地模拟船舶在碎冰区航行时相关的阻力性能以及冰块运动姿态。

3 计算结果及分析 3.1 敞水工况水动力导数数值计算结果

斜航运动选择漂角为0°, 2°, 4°, 6°, 8°，数值模拟采用相对运动模拟船舶斜航，计算来流速度为0.8 m/s，分别计算船舶的侧向力以及转艏力矩如图 9所示，其中船舶的侧向力是指沿着船宽方向的力。用最小二乘法分别拟合无因次后的侧向力F_y、转艏力矩N_z，所得曲线在原点处斜率即为水动力导数。侧向力F_y、转艏力矩N_z的计算结果见表 2。

将船舶在不同漂角所受的侧向力以及转艏力矩依据式(3)进行无因次化后拟合得到的曲线如图 10所示。

$ \left\{ {\begin{array}{*{20}{l}} {N_z^\prime = \frac{{{N_z}}}{{0.5\rho {V^2}{L^3}}}}\\ {F_y^\prime = \frac{{{F_y}}}{{0.5\rho {V^2}{L^2}}}}\\ {{v^\prime } = \frac{v}{V} = - \sin \beta } \end{array}} \right. $

(3)

斜航敞水工况水动力导数计算结果以及统计公式计算结果(采用井上模型)见表 3，可以看出Y′ _v、N′ _v相对偏差为20%、5.2%，Y′ _v误差偏大的原因是由于处于斜航运动模式时变化漂角时船舶力矩的变化较侧向力更为敏感求解结果更为准确，并且统计公式采用的是井上模型，井上模型是根据10艘各种类型的船模(油轮3艘，货船3艘，集装箱船、液化天然气船、滚装船、汽车轮渡各1艘)进行相关试验整理出的水动力导数计算公式，本次计算模型采用的是冰区某型破冰船其型线以及舱室布置、尾部形状都进行了改进与传统船型存在一定的差别，综上考虑文章数值计算方法具有一定的有效性。

3.2 碎冰工况冰水耦合后水动力导数数值计算结果

碎冰工况的斜航运动参数设置与敞水工况一致，漂角为0°, 2°, 4°, 6°, 8°。船体表面设置为壁面，物理模型中选择拉格朗日混合模型中多相耦合模块达到水与冰块流固耦合作用。多相相互作用设置冰- 冰接触、冰- 船接触、冰- 水接触、冰- 空气接触和水- 空气接触。通过冰- 船耦合作用改变冰块运动姿态进而冰块对船体周围流场产生干扰从而达到船冰水相互耦合。阻力模型选用Hertz Mindlin，阻力系数模型使用Schiller-Naumann模型，考虑计算耗时问题，本次模拟为碎冰低密集度40%下水动力相关模拟，经调研某型破冰船在浮冰工况下(浮冰尺寸大于50~60 cm)，连续破冰航速是：连续自破冰下，航速不超过18.52 km/h，0.2~0.3 m的浮冰基本不受影响，考虑本次模拟情况为低密集度小型碎冰并且不存在破冰船引航，选取船舶速度为18.149 6 km/h，依据弗汝德数相似换算得来流速度为0.8 m/s。分别计算冰水耦合后船舶的侧向力以及转艏力矩，如图 11~13所示。

3.2.1 冰水耦合下斜航侧向力

从图 11中可以看出在每一个漂角工况下，船舶的侧向力都呈现了上、下波动的情况，这是由于冰块与船发生碰撞作用时冰块的运动姿态发生变化，部分会翻转以及沿着船表面滑行还有一部分会被波浪排开，这些都对流场造成一定规律性扰动。本次数据处理考虑冰块的干扰具有随机性所以采用干扰的最大值以及最小值分别求解对应的水动力导数，从而形成一个水动力波动区间，更为真实地预报冰水耦合下的船舶水动力导数。其中最大以及最小值求取方法采用将15~33 s时间分为12个时间段，第i个时间段水动力波动的峰值为A_i、B_i，则整个时间段水动力波动峰值为$A = \frac{1}{{12}}\sum\limits_{i = 1}^{12} {{A_i}} 、B = \frac{1}{{12}}\sum\limits_{i = 1}^{12} {{B_i}} $侧向力具体值见表 4。不同漂角所受的侧向力进行无因次化拟合后得到的曲线如图 12所示。碎冰工况冰水耦合水动力导数计算结果见表 5。

3.2.2 冰水耦合下斜航转艏力矩

从图 13中可以看出在每个漂角工况下，船舶转艏力矩同样会由于冰块的运动姿态发生改变(主要包括部分会翻转、沿着船表面滑行以及一部分会被波浪排开)对船体周围流场造成干扰从而导致力矩值上下波动的情况。力矩数据处理方法与侧向力方法一致。转艏力矩具体值见表 6。将浮冰工况下船舶在不同漂角所受的转艏力矩进行无因次化后拟合得到的曲线如图 14所示。

由表 3、表 5和表 7中水动力导数计算结果可以看出碎冰工况冰水耦合下最大侧向力得到的线性水动力导数值F′ _ymax小于敞水计算结果，最小侧向力得到的线性水动力导数值F′ _ymin大于敞水计算结果。最大转艏力矩得到的线性水动力导数值N′ _zmax小于敞水计算结果，最小转艏力矩得到的线性水动力导数值N′ _zmin略大于敞水计算结果，线性水动力导数区间端点的绝对值存在一端数值大于敞水计算结果一端数值小于敞水计算结果的现象，并且非线性水动力导数计算结果呈现正负值。这是由于船舶在碎冰中航行时冰块与船体作用使得冰块运动姿态发生变化，冰块会在船体周围翻转、在船艏的堆积以及由于压力差紧贴船体表面的滑行，使得冰面对船体周围的流场造成了随机性的干扰，非线性水动力导数计算结果呈现正负值可看出对流场的干扰更具有随机性，并且非线性水动力导数值在冰水耦合下不容易捕捉求解，波动灵敏性太大。

4 结论

1) 本文通过与本次数值模拟具有相同网格以及DEM相关参数设置的76 K冰区散货船碎冰区的数值模拟，从静水阻力、碎冰阻力以及碎冰与船舶艏尾运动姿态与实验进行了对比验证，并将敞水下船舶进行斜航模拟得到的水动力导数值与井上模型统计公式数值对比，数值较为接近，进一步说明了数值方法的可靠性。

2) 船舶的侧向力以及转艏力矩在冰水耦合的情况下都呈现了上、下波动的情况，这是由于冰块与船发生碰撞作用时冰块的运动姿态发生变化，部分会翻转以及沿着船表面滑行还有一部分会被波浪排开，这些都对流场造成一定规律性扰动。

3) 水动力导数计算结果可以看出最大侧向力以及最大转艏力矩得到的线性水动力导数值F′ _ymax、N′ _zmax小于敞水计算结果，最小侧向力以及最小转艏力矩得到的线性水动力导数值F′ _ymin、N′ _zmin大于敞水计算结果，这是由于船舶在碎冰中航行时冰块与船体作用使得冰块运动姿态发生变化，冰面对船体周围的流场造成了随机性的干扰。随机性取其区间值能体现冰块对流场干扰的随机性。

4) 本文使用CFD-DEM方法虽然一定程度上解决了冰水耦合的问题，但还有许多设置需要改进，同时受到时间和计算能力等的限制，本文只计算了一种密集度下冰水耦合的情况，后期将针对于不同密集度以及碎冰的不同尺寸部分进行研究。