基于SPH 方法的黄河破冰船冰阻力数值模拟分析

郑兴，田治宗，谢志刚，张宁波

1 哈尔滨工程大学 船舶工程学院，黑龙江 哈尔滨 150001

2 黄河水利科学研究院，河南 郑州 450003

0 引 言

黄河由于其所在地理位置纬度较高，冬季严寒，因而在冬季会出现结冰现象，特别是在黄河最北端的内蒙古流域，在封河初期以及开江解冻期，浮冰在水流的作用下易在河流狭窄弯曲段出现堆积，形成冰排、冰坝，从而堵塞河道，造成凌汛灾情，给沿岸人民的生活和财产造成巨大损失[1]。为了应对冰凌灾害，需要采取破冰措施进行疏通，目前的破冰措施主要有爆破破冰[2]、破冰船破冰[3]等方法，其中破冰船因具有机动性强、成本低等优点，能在部分河段的破冰疏通方面发挥重要作用。而采用有效的数值分析方法对黄河破冰船的冰阻力进行估算，对于破冰船的设计和操纵具有重要作用。

近年来，国内外一些学者采用数值模拟方法对冰与结构物的相互作用问题进行了分析，现有的用于模拟冰破坏的数值方法主要有有限元方法（FEM）、离散元方法（DEM）、近场动力学（PD）方法以及光滑粒子流体动力学（SPH）方法等。 Kolari等[4]提出了一种基于模型更新技术的冰连续破坏有限元模拟方法，该方法是采用各向异性连续损伤力学（CDM）模型来预测裂纹的扩展方向，同时也将模型更新技术运用到了有限元模型中裂纹的扩展发展中。Carney 等[5]建立了冰的现象学破坏模型，并采用显式动力分析程序LS-DYNA 中的有限元程序对这种冰的塑性敏感性破坏模型予以了求解。Pernas-Sánchez 等[6]建立了基于Drucker-Prager 屈服准则的弹塑性材料模型，得到了冰在高应力率作用下的破碎行为，同时通过采用有限元程序LS-DYNA 中的求解模型，并将拉格朗日、任意拉格朗日−欧拉（ALE）方法和SPH 这3 种方法集成应用到模型研究中，得到了冲击力对细长冰柱的影响。孔帅等[7]基于高性能DEM 方法，对极区浮式平台的冰载荷进行了数值分析。狄少丞[8]建立了用于模拟海冰与两单元接触的离散元模型，研究了冰的单轴压缩和三点弯曲破坏过程，并与试验结果进行了对比，显示吻合较好。薛彦卓等[9]应用近场动力学理论建立了冰材料的弹脆性破坏模型，并针对冰梁的三点弯曲破坏进行了模拟研究，其数值模拟结果与试验数据的对比显示其一致性较好。

针对破冰船的破冰阻力预测，可以采用SPH方法进行计算分析。SPH 方法属于无网格拉格朗日粒子算法，与传统的网格方法相比，无网格粒子算法的求解不依赖于网格，而是通过一系列的节点（粒子）来求解整个系统，且这些粒子点之间不需要网格进行连接，可以自由运动，特别适合分析大变形和不连续的自由表面问题。另外，SPH 方法在模拟固体力学方面因为粒子的拉格朗日特性，粒子在发生屈服后会自然生成裂纹，因此SPH 方法也非常适合模拟固体大变形以及破坏问题。随着SPH 方法在流体力学和固体力学领域的快速发展，该方法也被用于模拟冰动力学、冰破坏以及冰与结构物的相互作用方面。例如，Das[10]在LS-DYNA 软件中使用SPH 模型对冰梁的四点弯曲破坏进行了模拟，并采用冯米塞斯屈服准则针对冰粒子的破坏失效予以了判别，结果显示采用该方法一旦冰粒子失效，偏应力分量将直接降为0。Zhang 等[11]使用一种改进的SPH方法，并结合D-P 屈服准则和黏聚力软化模型的弹塑性本构模型，对冰的弯曲和压缩破坏特性进行了模拟。Gutfraind 等[12]和Oger 等[13]在SPH 方法中应用基于Mohr-Coulomb 屈服准则的流变学，对在水面上漂浮和在风力作用下移动的碎冰场进行了模拟。Shen 等[14]提出了一个二维河面中冰的漂移和阻塞数值模型，并采用SPH 方法针对河面冰的漂移运动进行了求解与模拟。Ji 等[15]提出了一种新的海冰动力学黏弹塑性本构模型，其中SPH 方法被用于模拟矩形海盆中冰的运动。此外，Ji 等[16]还发展了一种混合拉格朗日−欧拉（HLE）海冰动力学方法，其中海冰覆盖层用具有自身厚度和密集度的SPH 模型表示。Pan 等[17]提出了一个新的SPH 非牛顿模型，用于研究冰盖和冰架的耦合动力学问题。乔岳[18]应用SPH 方法建立了船−冰−水耦合作用数值模型，并对冰阻力预报方法予以了分析，不过其主要关注的还是冰阻力的预报方法，对于冰−水耦合模型也只是采用了简单的排斥力模型，未考虑碎冰区以及有波浪作用下船−冰−波的相互作用情况。随后，Zhang 等[19-20]应用SPH 方法建立了船−冰−波耦合数值模型，模拟研究了冰的破坏过程和波浪−冰的相互作用过程，并预测了船体所受的冰载荷。以上一系列研究表明，使用SPH 方法研究破冰阻力可行且有效。

根据SPH 方法的无网格特性以及其在流体问题、固体材料大变形及破坏问题，以及流−固耦合问题上的适用性与优势，本文将基于SPH 方法建立船−冰−水耦合数值模型，并对黄河破冰船的冰阻力进行预报。该方法的构建包括冰的弹塑性本构方程、船−冰−水数值模型、破碎冰环境以及与黄河河道中流场的作用。通过与数值计算结果的对比验证，对黄河破冰船冰阻力相关因素的影响规律进行分析，以为破冰船的设计和操纵提供指导。

1 数学模型的建立

1.1 SPH 的基本方程及相关数值方法

流体动力学的基本控制方程包括质量守恒方程和动量守恒方程，具体如下：

式中：D/Dt为物质导数；α，β 为x，y，z方向上的笛卡尔分量；ρ 为粒子密度；ν 为粒子速度；g为重力加速度；σαβ为应力张量，xβ为位置坐标在β 方向上分量；t为时间。本文中的连续性方程引入了Antuono 等[21]提出的一种人工扩散项，主要用于消除流体模拟中压力场不合理的高频振荡。为了提高数值计算的稳定性，减小不稳定振荡，引入了Monagham 型人工黏度耗散项。固体相的质量守恒方程粒子可近似表达式为

式中： ρi为具有速度分量vi的粒子i的密度； ρj，mj分别为具有速度分量vj的粒子j的密度和质量；Wij为SPH 方法中的光滑核函数。冰粒子动量方程的离散形式为

1.2 冰的破坏模型

为了模拟冰的破坏过程，本文将弹塑性本构模型应用到了SPH 方法中。结合Drucker-Prager屈服准则，得到具有非关联流动规律的冰模型应力应变方程的应力应变关系为：

在破坏过程中，还需要对应力张量进行塑性修正，需要采用黏聚力软化方法降低塑性流动阶段冰的黏聚力，最终实现对冰破坏过程的模拟。由于冰的破坏模型不是本文的研究重点，故在此不予展开，具体的过程可参见文献[20]。

1.3 船−冰−水耦合模型

在冰−水的SPH 耦合计算模型中，将破冰船当作流体相和冰相的固壁边界来处理，采用虚粒子边界方法。该方法是先对破冰船边界建模，然后划分网格，随后在这些网格节点上布置SPH 粒子以实现对固壁边界的布置。这些船体边界粒子作为边界虚粒子，将分别参与到冰相和流体相的求解中。在求解冰相控制方程时，冰粒子计算域中的船体固壁边界粒子可以作为虚粒子来施加边界条件。这些船体固壁边界粒子将参与冰相的连续性和动量方程求解，这样，船体边界附近的冰粒子将满足应力和速度的连续性。在求解流体控制方程时，流体粒子计算域中的船体固壁粒子也将充当流体相的边界虚粒子。这些固壁粒子同样将被计入流体相的连续性和动量守恒方程求解中，用来满足流体相的压力和速度边界条件。

由于船体边界粒子的速度与船体速度相同，当这些船体边界粒子作为虚粒子来近似船舶边界与冰之间的耦合作用时，需要从其邻域内的冰粒子中插值计算出这些边界虚粒子的应力，用以准确近似边界附近冰粒子的应力梯度。然后，再从流体邻域内的体粒子中插值计算得到这些固壁边界虚粒子的压力，用以近似船体边界附近流体粒子的压力梯度。该计算过程可参见图1。图中：i，f分别为第i号和第f号粒子；kh为粒子支持域的影响范围，其中k= 2.0，取1 倍粒子尺寸，h为光滑长度。

图1 船−冰−水耦合模型示意图Fig. 1 Schematic diagram of ship-ice-water coupling model

当冰层发生破坏时，破坏位置处冰粒子的密度值变化较大，为了保证密度的连续性，在船−冰耦合求解中，船体边界粒子的密度根据其领域内冰粒子的密度插值计算得到。同理，在船−水耦合求解中，船体边界粒子的密度也可以由其邻域内流体粒子的密度插值得到。此外，还可通过忽略船体边界虚粒子与周围邻域内冰粒子和流体粒子的黏性作用，来施加可滑移边界条件。通过上述方法，可实现流体相、冰相与固壁边界的耦合。这种边界处理方法的主要优点是计算过程简单，不需要计算船体边界和耦合界面复杂的几何信息。船体表面的边界粒子应足够密，以防止冰粒子和流体粒子穿透船体壁面边界[19]。需要指出的是，由于直接采用SPH 方法进行大尺度数值模拟计算效率较低，为了提高计算效率，在上述耦合模型中只考虑施加了流体对冰的作用力，忽略了冰对水的作用。此外，数值模型中也未考虑船体的运动响应或是力学响应。

另外需要指出的是，在数值模型中均是采用虚拟粒子法来处理计算域的边界[20]。例如，在后文的数值模拟中，对于平整冰的模拟计算，需沿着平整冰面布置相应的边界虚粒子；而在黄河河道碎冰条件的数值模拟中，则是沿着河道边界布置对应的边界虚粒子。

2 平整冰情况下的数值模拟与验证

本节将采用SPH 耦合计算模型，针对使用区域黄河上冰面的特点，先给出平整冰情况下的冰阻力计算结果，用以为实际河道中破冰船冰阻力的预测提供参考。

2.1 船模模型建立

所研究的轻型黄河破冰船的总布置图如图2（a）所示，其三维模型如图2（b）所示，主要参数如表1 所示。

图2 28 m 黄河破冰船模型示意图Fig. 2 Schematic diagram of a 28 m ice breaker model operating on the Yellow River

表1 28 m 黄河破冰船主要参数Table 1 Main parameters of a 28 m ice breaker operating onthe Yellow River

为了采用SPH 方法进行模拟，首先在船体表面建立均匀分布的固壁粒子来表示船体边界的作用，然后通过求解冰粒子对船体固壁边界粒子的作用力，得到船−冰耦合过程中船体所受的冰载荷。当船体边界粒子作为冰粒子的边界虚粒子时，可以利用相邻冰粒子的应力张量的体积积分来估算冰粒子对船体边界粒子的作用力。冰粒子对船体边界粒子的作用力可以表示为

选取2015年5月～2017年5月我院ICU接受机械通气治疗的新生儿140例作为研究对象，其中，男79例，女61例，平均胎龄（35.36±2.77）d，平均出生体重（2014.05±160.77）g，呼吸窘迫综合征57例，病理性黄疸49例，高胆红素血症21例，缺氧缺血性脑病13例。将其随机分为研究组与对照组，各70例，两组一般资料比较，差异无统计学意义（P＞0.05）。

2.2 平整冰面冰阻力的模拟与验证

在给出SPH 方法的计算结果之前，将引入Lindqvist 经验公式对所计算的冰阻力进行验证。基于冰阻力与速度呈线性关系的假设，Lindqvist公式将冰阻力分为3 个部分进行计算，包括挤压阻力、浸没阻力和弯曲破坏阻力[23]，相关公式可参见文献[23]。Lindqvist 计算公式的主要参数如表2 所示。

表2 Lindqvist 计算公式的主要参数Table 2 Main parameters of Lindqvist formula

表3 V = 4 kn，hi = 0.3 m 时不同弯曲强度情况下的破冰阻力预报精度比较Table 3 Comparison of prediction accuracy of ice breaking resistance with different bending strength when V = 4 kn and hi = 0.3 m

图3 V = 4 kn，hi = 0.3 m 时不同弯曲强度下破冰总阻力的时间历程比较Fig. 3 Time histories comparison of ice breaking resistance with different bending strengths when V = 4 kn and hi = 0.3 m

为了研究不同冰厚情况下的破冰阻力，图4给出SPH 方法计算的不同冰厚和不同弯曲强度下冰厚对冰阻力的影响。由图4 可以看出，冰阻力的变化与弯曲强度大致成线性增长，且冰阻力随冰厚的变化也较为明显，当冰厚增大1 倍，其他条件不变时，冰阻力将增大2.5 倍左右，这充分说明冰厚是影响船舶破冰阻力的最重要因素。图5（图中， ε¯P为应变率幅值）给出了不同冰厚情况下当t= 20 s 时平整冰情况下的破冰形状。总体来看，破冰船在破冰过程中在船艏前会形成周向裂纹，随后，冰层会遵循相同的破裂模式，这与图3中冰阻力呈周期性变化的特点是相对应的。此外，在船体的左舷和右舷也发生了一些局部冰被挤压破碎的现象，这与实际破冰船在平整冰中破冰过程中的类似现象相符[24]。另外，在船后形成的水道中也存在一些碎浮冰。结果表明，破冰过程并不总是两侧对称的。由图5 可以发现，当冰厚较小时，船艏附近冰的裂纹扩展不明显，破碎后的块状浮冰边界不清晰；但随着冰厚的增加，船艏附近冰的破碎程度增加，在出现断裂后有清晰的块状浮冰。这说明采用SPH 方法不仅可以预报冰阻力，还可以预报破碎后的冰裂纹情况。

图4 不同弯曲强度下冰厚对冰阻力的影响Fig. 4 Influence of ice thickness on ice breaking resistance at different bending strength

图5 不同冰厚情况下t = 20.0 s 时的破冰形状Fig. 5 Ice breaking shapes at different ice thickness when t = 20 s

3 河道中冰阻力模拟与分析

为了研究实际河道中破冰船的冰阻力，本节将对河道中碎冰区的船舶冰阻力进行分析。首先，根据黄河某段河道的卫星图片（图6（a））绘制出河道的边缘模型（图6（b））；随后，为了建立三维河道的剖面图形，先采用Solidworks 软件建立河流出口与入口这两处河道的断面，然后通过河道两侧地形进行诱导拉伸，得到每个不同位置处的剖面形状。在模拟河道中冰的堆积与运动情况时，需要有关河道流场流速方面的分布信息。由于直接采用SPH 方法对大尺度的河道流速进行求解有困难，本文采用OpenFOAM 中的icoFoam求解器予以求解。用该求解器求解非稳态不可压缩N-S 方程，离散格式采用有限体积法。河道网格采用blockMesh 工具生成六面体结构网格，使用snappyHexMesh 工具读取河道表面stl 文件后画出河道轮廓并细化网格。随后，在河道进口处给定流速U=0.5，0.7，1.0 m/s，即可获得给定流速对应的河道流场分布。

图6 河道模型的建立Fig. 6 Setting up of the river channel model

建立好河道边界后，需在河道区域内建立浮冰的SPH 粒子模型。该碎冰模型可采用Voronoi图来建立，其建模原理如下：首先，根据碎冰航道的面积以及碎冰密集度和浮冰的平均尺寸，计算出该碎冰区域内浮冰块的数量，其中碎冰密集度根据实际的碎冰覆盖面积与相应冰区总面积之比得到，浮冰平均尺寸为碎冰覆盖总面积与碎冰数量之比；然后，在所需计算的河道模型区域内生成相同数量的随机样点，根据随机点生成Voronoi多边形，以保证多边形内任何位置点到该多边形随机样点的距离最近；最后，根据选定的碎冰密集度对Voronoi 多边形进行收缩，以确定SPH 粒子的填充区域。根据该方法生成的河道碎冰区计算模型如图7 所示。

图7 碎冰建模示意图Fig. 7 Setting up of broken ice model

图8 所示为碎冰密集度C= 70%时破冰船在河道中作业时的模拟结果。碎冰计算模型中，浮冰的平均尺寸为12.1 m2。该船从河道入口出发，行驶轨迹与河道主航道重合，航速保持V= 4 kn不变。由图可见，采用SPH 方法能够对碎冰在河道中的漂移进行有效模拟，而且在破冰船的作用下，船行驶后尾部有明显的开阔水道。为了对不同碎冰密集度下的破冰阻力进行分析，图9 给出了航速V= 4 kn、冰厚hi= 0.3 m 时不同碎冰密集度下破冰船冰阻力的时间历程比较。由图9 可以看出，当碎冰密集度较小（C= 50%）时，浮冰对船体破冰阻力的作用比较小，当碎冰密集度增大到C= 90%时，浮冰对船体破冰阻力的作用明显增大。另由图9 所示的时历曲线还可以发现，当碎冰密集度较小时，冰阻力的作用时间比较短，且每次作用时冰阻力的最大幅值也较小。表4 给出了不同碎冰密集度下的冰阻力平均值。可见随着碎冰密集度的增加，冰阻力的作用时间逐渐变得连续，冰阻力的平均幅值也逐渐增大。碎冰密集度对破冰船破冰阻力作用时间的分布有着明显的影响。

图8 碎冰密集度为70%时破冰船在河道中作业的模拟结果Fig. 8 Numerical results of icebreaker working in the river when the concentration of broken ice C = 70%

图9 V = 4 kn，hi = 0.3 m 时不同碎冰密集度下破冰船总阻力的时间历程比较Fig. 9 Time histories comparison of breaking ice resistance with different concentrations of broken ice when V = 4 kn and hi = 0.3 m

表4 不同碎冰密集度情况下SPH 计算的冰阻力平均值Table 4 Average ice resistance by SPH with difference concentrations of broken ice

为了研究碎冰情况下不同冰厚对破冰船冰阻力的影响。图10 给出V= 4 kn、碎冰冰密集度C= 70%时，采用SPH 方法计算得到的不同冰厚下破冰船冰阻力的时间历程比较。图中数值结果对应的模拟时间历程为120 s。由图可见，冰厚对冰阻力的曲线形状变化影响不明显，但对冰阻力的幅值变化影响显著。当hi= 0.3 m 时，冰阻力的最大幅值为2.73×105N；当hi= 0.4 m 时，冰阻力的最大幅值为3.11×105N；当hi= 0.5 m 时，冰阻力的最大幅值为4.93×105N，当hi= 0.6 m 时，冰阻力的最大幅值增大为5.81×105N。由图可知最大冰阻力的增幅是冰层厚度增幅的2.7 倍，说明冰厚是影响破冰阻力的主要因素。表5 示出了不同冰厚对冰阻力的影响。当hi= 0.3 m 时，冰阻力的平均值为31.88 kN；当hi= 0.4 m 时，冰阻力的平均值为33.53 kN；当hi= 0.5 m 时，冰阻力的平均值为48.48 kN；当hi= 0.6 m 时，冰阻力的平均值为61.03 kN。可见冰阻力的平均值相对于最大幅值降低了1 个数量级，此时，平均冰阻力的增幅是冰层厚度增幅的1.91 倍。

表5 不同冰厚情况下SPH 计算冰阻力平均值Table 5 Average ice resistance by SPH with difference ice thicknesses

图10 V = 4 kn，C = 70%时不同冰厚下破冰船总阻力的时间历程比较Fig. 10 Time histories comparison of breaking ice resistance with different ice thicknesses when V =4 kn and C = 70%

为了验证SPH 结果的可靠性，采用文献[25]中有关碎冰中的船舶冰阻力进行了验证，计算公式如下：

表5 关于不同冰厚条件下SPH 计算结果与经验公式的比较说明，采用SPH 方法模拟河道中碎冰的冰阻力具有一定的合理性。

4 结 语

本文首先对SPH 方法进行了拓展，通过结合冰的弹塑性本构方程和河道的建模方式，完成了船−冰−水耦合数值模型的建立；然后，采用该数值模型对平整冰面中黄河破冰船的冰阻力进行预报，并与Lindqvist 经验公式进行比较验证了其有效性；接着，将该数值模型拓展应用到了黄河河道中破冰船的冰阻力预报中；最后，通过对不同情况下冰阻力计算结果的比较，分析了黄河破冰船冰阻力的规律和特性。主要得到如下结论：

1） 采用SPH 方法能够较好地预报黄河破冰船在平整冰面情况下的破冰阻力，通过与Lindqvist经验公式的对比，显示采用SPH 方法得到的冰阻力误差结果小于17.6%。

2） 可以采用Voronoi 图的方法建立破碎冰面的SPH 粒子模型，该模型可以通过对碎冰密集度的调整来实现对不同碎冰冰面初始时刻的划分，能为河道中碎冰的运动及其对船舶的碰撞作用提供研究基础。

3） 借助icoFoam 求解器求解的河道流场，采用SPH 方法能够实现对真实河道中碎冰运动和冰阻力的模拟计算。通过对不同碎冰密集度情况下冰阻力的对比研究，发现在该河道中碎冰条件下冰阻力的作用方式与平整冰面不同，河道中破冰阻力的作用时间历程变化非常剧烈，且作用时间较短。此外，冰厚也是影响冰阻力最重要的一个因素。

以上规律和特性可为黄河破冰船的实际操作提供参考。