冰区船艏部连续破冰模式下破冰阻力数值计算研究

2022-07-25刚旭皓田于逵季少鹏余朝歌

船舶力学 2022年7期

刚旭皓，田于逵，季少鹏，寇莹，国威，余朝歌

（中国船舶科学研究中心，江苏无锡 214082）

0 引言

冰区航行船舶是北极探索和开发过程中不可或缺的基本装备。冰区船舶冰阻力模拟既能为冰区航行船舶的初步设计提供参考，又能为船-冰相互作用研究提供依据，避免船舶在航行过程中因主机功率不足而随流冰漂移或倾覆，因此准确模拟船舶在冰区航行时的破冰过程并预报船舶冰阻力具有重要的理论和现实意义。

随着计算技术的不断进步，数值计算方法逐渐成为冰区船舶冰阻力性能研究的有效工具。Wang[1]将破冰过程理想化为连续的接触-挤压-弯曲破坏的循环，提出了时域内数值流程的解决方法；Sawamur[2]提出了一种用于循环计算冰阻力的连续接触流程，模拟了碎冰受到外力作用后的动态响应；Su[3]将破冰的过程分为破冰和碎冰运动两个部分，建立了冰载荷与船体运动的三自由度耦合方程，并将计算结果与破冰船AHTS/IB 试验数据对比，显示了很好的一致性。海冰的离散元模型从20世纪80年代发展起来并不断完善，并在涉冰研究领域取得了一定的研究成果。Hopkins[4]等将冰块离散为随机尺寸的圆盘，模拟了海冰在斜坡式结构前的堆积过程和压力冰脊的形成；Lau[5]在此基础上，对平整冰与锥体结构相互作用中的非线性大变形和断裂过程进行模拟，编写了三维离散元程序DEMICE 3D；季顺迎和狄少丞等[6-7]将基于GPU的高性能计算方法引入离散元数值模拟中，显著提高了DEM方法计算规模和计算效率,并将其应用于破冰船的研究。除了上述数值计算方法之外，理论分析方法和边界元方法也在船-冰相互作用领域有一定的应用，但由于船-冰相互作用过程的复杂性和随机性较强，还需进一步完善。国内冰阻力预报数值研究已经取得了较大的进步，但由于模型试验和实船验证不充分，预报的准确性和适用度相对不足。只有将冰阻力预报与海冰物理力学特性及试验验证紧密结合，才能建立起有效的冰阻力预报方法，为冰区船航行性能与船型设计优化等工程应用提供扎实、有力的技术手段。

本文基于中国船舶科学研究中心小型冰水池（CSSRC SIMB）中开展的连续破冰模式下冰区船艏破冰阻力模型试验[8]，根据离散元模型的基本原理，结合元胞数组方法，利用冰力学试验数据精细化参数，建立离散元数值计算模型Bri-DEM，对船艏模型的破冰阻力及冰层破坏模式进行数值计算，分析船艏倾角和冰厚变化对破冰阻力的影响，获得破冰阻力相应的变化规律。

1 基于离散元方法的船艏破冰阻力数值计算

1.1 基本原理

在离散元方法中，海冰单元的力学行为是牛顿第二定律与接触处力-位移定律交替应用的结果：牛顿第二定律用于确定各质点在接触力、阻尼力和物体力作用下的平动和旋转，而力-位移定律则用于更新各接触点相对运动产生的接触力以及单元的位置信息。为模拟海冰的连续特性，相邻单元之间采用平行粘结模型粘连。平行粘结模型的基本原理是设想一组具有恒定法向和切向刚度的线性弹簧，以两单元的接触点为中心，均匀分布在粘结圆盘上，完成单元间作用力和力矩的传递。单元之间的粘结失效通过梁理论进行求解：

式中，Fn、Fs分别为作用在粘结圆盘上的法向力和切向力，Mn、Ms分别为作用在圆盘上的法向力矩和切向力矩，σmax、τmax分别为作用在粘结圆盘上的最大的拉应力和剪应力，I、J分别为粘结圆盘的惯性矩和极惯性矩，A为粘结圆盘面积。当圆盘上的最大拉伸应力σmax超过法向粘结强度，或最大剪切应力τmax超过切向粘结强度时，颗粒间的粘结键将会被破坏，由此模拟内部裂纹产生的情况[9-10]。

在单元粘结处，单元受力和位移的关系可以通过以下参数来表示：法向刚度Kn和切向刚度Ks、摩擦系数μ、法向阻尼和切向阻尼，以及粒子重叠时在重叠区域中心形成的接触区。

单元的法向接触刚度和切向接触刚度对单元间的受力具有重要的影响，单元的接触刚度与接触处的有效弹性模量相关。在冰单元与结构物相互作用时，单元与结构之间的有效弹性模量Ee可表示为

式中，Es为结构物的弹性模量，Ei为冰的弹性模量，v为粘性系数。冰单元的法向接触刚度Kn可表示为

式中，hi、ρi表示冰厚和冰密度，l为冰特征长度，g为重力加速度。根据Ji[11]的研究，Ks与Kn之间存在线性关系，切向刚度Ks可表示为

式中，r根据弹性模量和切变模量的关系确定，取r=1/2（1+v）。

对接触处进行力学分析，接触力Fi可分解为切向力和法向力：

式中，Fn、Fs分别表示接触力的法向和切向分量，ni、ti为单位向量。根据Hertz接触理论，法向力Fn可以表达为

式中，xn表示单元之间的重叠量，vn为法向相对速度，Cn为单元间的法向阻尼系数，可表示为

式中，ξn为单元间的阻尼比，M为单元的等效质量，ei为回弹系数。ei可参考Ramírez[12]提出的不同碰撞速度下的回弹系数计算，如图1所示。

图1 不同碰撞速度下的回弹系数[11]Fig.1 Rebound coefficients at different speeds

单元间的切向力Fs一般以增量的形式进行计算。剪应力的增量为每个时间步内切向方向的变化量与切向刚度的乘积：

式中，Δxs和vs分别为两个单元在接触点处的切向位移增量和相对切向速度，Cs为切向阻尼系数。假设单元间的接触力满足摩尔-库伦定律，则有

1.2 船艏及冰层离散化

1.2.1 船艏离散化

设计破冰型冰区船艏部模型，利用三维建模软件UG进行数值化建模，利用网格划分工具Hypermesh 将模型合理地划分为一系列三角离散单元，图2 是船艏模型离散化示意图，表1 为船艏模型基本参数。提取网格节点信息导入数值计算模型，通过判断冰单元与船艏模型网格单元之间的接触来判断冰单元受到的接触力和力矩，以此对船艏受到的作用力进行求解，进而计算破冰阻力。

表1 船艏模型基本参数Tab.1 Basic parameters of bow model

图2 船艏模型离散化示意图Fig.2 Schematic diagram of bow model discretization

1.2.2 冰层离散化

利用矩形单元采用最密排列的方式对冰层进行构造。对冰层进行离散化后，将单元编号为k1、k2，…，kn，采用元胞数组储存颗粒的位置和运动信息，随着时间步更新，利用元胞数组特性对粒子的位置和运动信息进行实时更新，可实时对粒子间的相互作用力进行计算，简化了颗粒搜索过程，同时也利于对冰单元kn与船艏相互作用力以及单元与单元之间相互作用力的计算。

冰单元kn与船艏模型相互作用时，作用力主要有单元与船艏之间的挤压力以及单元在船艏表面滑移产生的阻尼力以及摩擦力。对冰单元之间力的传递方式进行设置，冰单元与船艏接触的位置不同，力在冰单元之间的传递也会有所不同，图3表示不同接触位置接触力的传递方式。当船艏与冰层左侧的冰单元接触时（k1、k2），冰单元之间的接触力自左侧向右传递和向后传递；当船艏与冰层右侧的冰单元接触时（k4、k5），接触力自右侧向左传递和向后传递；当船艏与中部的冰单元接触时（k3），冰单元之间的接触力自中间向两侧传递及向后传递。在完成第一步传递之后，接触力逐渐拓展，直到覆盖到整个冰层，计算每一个冰单元受到的力。通过对接触力的不同传递方式进行模拟，以此来计算冰单元之间力的作用以及冰层的裂纹生成和断裂。

图3 不同接触位置接触力的传递方式Fig.3 Transmission mode of contact force at different contact positions

1.2.3 计算参数输入

考虑到冰的力学特性对冰单元之间的粘结作用有较大的影响，为了确定计算参数，对模型冰的力学特性在不同加载速率下的变化进行了研究。研究发现：冰的弯曲强度随加载速率的增加基本没有发生变化；而冰的压缩强度，在试验速度范围内呈现出随着加载速率增加而逐渐减小。图4为弯曲强度与压缩强度随加载速率的变化。为避免偶然性结果，对同一工况进行了重复性实验，图中使用不同形状标记表示。根据不同速度下冰强度值确定单元间的粘结强度，将冰的物理力学参数和计算得到的粘结强度输入计算模型，具体的计算参数如表2所列。

表2 船艏模型破冰阻力数值计算输入参数Tab.2 Input parameters for ice-breaking resistance calculation of bow model

图4 弯曲强度与压缩强度随加载速率的变化Fig.4 Variation of bending strength and compression strength with loading rate

本次数值计算冰厚取40 mm，船艏模型与冰相互作用的速度设置为0.01 m/s、0.05 m/s、0.10 m/s、0.15 m/s 和0.20 m/s，冰层的边界条件设置为三边约束、一边自由的状态，自由边为船艏模型进入冰层的方向，与冰水池模型试验的情况一致。在数值计算模型中，将船艏模型的网格模型导入计算程序中，设置船艏与冰层的相对位置，如图5所示。

图5 三维船艏与冰层的设置形式Fig.5 3D setting of the bow and ice sheet

1.3 计算结果分析

利用Bri-DEM 程序对船艏与冰相互作用过程进行模拟，获得了船模破冰阻力，同时对破冰过程及碎冰运动进行了模拟，图6为不同时刻冰单元的运动情况。

图6 不同时刻冰单元的运动情况Fig.6 Movement of ice units at different times

破冰过程中，船体艏柱线首先与冰层接触，对冰单元产生挤压作用，引起冰单元旋转，同时冰单元之间的相互作用力逐渐向后、向两侧延展，引起冰单元运动，直至冰单元之间的作用力超过粘结强度，引起冰单元断裂。船艏与冰层之间的接触随着船艏前进，由单点接触转变为多点接触，形成接触-挤压-破坏的周期性过程。由于未考虑水动力的影响，冰单元在冰层上断裂之后，会直接向下运动，且不会与冰层再次粘结。

船艏在冰层中破冰航行时产生接触-挤压-弯曲破坏的循环过程，从船艏受力的时历曲线中可看出，曲线具有明显的脉动特性和周期性，图7 为速度0.15 m/s 时船艏受力时历曲线。冰单元的受力与船艏与冰的重叠距离有关，在单一时间步内，随着船艏与冰单元之间的重叠面积增加，引起冰单元受力增加，船艏受力增加，但计算得到的重叠距离具有随机性，因此不同时间步的力值出现较大的波动。

图7 速度0.15 m/s时船艏受力时历曲线Fig.7 Time-history curve of bow force at a speed of 0.15 m/s

船舶破冰速度的大小决定了船-冰相互作用的速度和频率，是影响船体破冰阻力的一个重要因素，分别对不同速度的计算结果绘制力-速度曲线，如图8所示。随着破冰速度的增加，破冰阻力呈增大趋势。为验证计算模型准确性，将采用模型试验数据和经验公式对计算结果进行校验。

图8 不同速度下船艏破冰阻力Fig.8 Ice-breaking resistance of ship bow at different speeds

2 船艏破冰阻力数值计算结果验证

在冰阻力的数值预报中，能否准确地预报冰阻力值以及破冰模式，是评价数值计算模型优劣的重要指标。将Bri-DEM 计算结果与模型试验数据进行比较，可验证数值计算模型的可靠性。在此基础上，使用验证后的模型进行拓展计算，进一步对船艏冰阻力进行研究。

2.1 破冰效果比较分析

计算模型可将冰单元的运动及位置信息输出，获得冰层断裂和碎冰旋转浸没的过程。当船艏与冰层接触时，船体艏柱与冰层最先产生接触，冰层产生挤压破坏，随着船艏继续破冰前进，垂向力逐渐增加，迫使冰层向下弯曲破坏，碎冰从冰层上脱落；同时会对其相邻的冰单元产生影响，使其发生不同程度的变形，这种现象与试验过程中船艏与冰层接触产生挤压破坏并产生裂纹、最终发生弯曲破坏的过程较为一致，见图9。

图9 弯曲破坏阶段对比Fig.9 Comparison of flexural failure stages

随着船艏继续破冰前行，冰单元与船艏的接触逐渐由单点接触逐渐转换为多点接触，同时冰层内部单元之间的粘结作用会对航道两侧的冰层产生影响，如果航道两侧的冰层受力未达到断裂强度值，船艏破冰后，将恢复到原来的状态。与试验照片进行比对，数值模拟能够反映冰层的破坏与船艏对破冰航道外冰层的影响，如图10所示。

图10 冰层破坏模式对比Fig.10 Comparison of ice destruction patterns

由于在数值计算过程中未考虑水动力以及水的浮力的作用，冰单元断裂之后会直接脱落，而在模型试验过程中，碎冰会逐渐覆盖船底，并随船向后运动，如图11所示。

图11 冰层破冰效果对比Fig.11 Comparison of ice breaking effects

利用Bri-DEM 程序对冰层的破坏及碎冰的运动进行模拟，能够较好地反映出冰层的破坏过程，并可模拟出冰层首先破坏的位置及其对周围冰层的影响。

2.2 破冰阻力对比分析

在船艏与冰相互作用的过程中，随着破冰速度的增加，船艏受到的破冰阻力也在不断地变化。统计比较Bri-DEM 计算结果与模型试验结果，如图12 所示。在与模型试验比较的同时，利用Jeong 经验公式估算了相同工况下的破冰阻力。Jeong[13]公式是在Spencer公式的基础上提出的一种冰阻力经验估算方法，计算公式表达如下：

图12 数值计算结果与试验数据对比Fig.12 Comparison of numerical calculation results and test data

式中：Ri(v)、Ro(v)分别为冰阻力和敞水区域阻力；Cb、Ch和Cr分别为浮力系数、清冰力系数和破冰力系数；Fh和Sn分别为Froude 数和强度因子；p和q分别为Froude 数和强度因子的幂指数。可以看到，Bri-DEM 能够模拟出破冰阻力的变化趋势，即随着速度的增加，破冰阻力逐渐增大，Bri-DEM 计算结果与模型试验结果、Jeong 经验公式估算结果大体相近，且Bri-DEM 计算结果与模型试验结果二者较为一致。

3 多方案破冰阻力数值预报

冰区船在冰层中航行的破冰阻力受到诸多因素的影响，如冰的物理力学特性、船体形状、船舶航态以及航速等，选取船艏倾角和冰厚作为典型影响因素，研究其对破冰阻力的影响。

3.1 船艏倾角的影响

船艏倾角大小对破冰模式以及破冰时主要受力类型有重要的影响。当船艏倾角较小时，破冰的主要类型为弯曲破坏，而船艏倾角较大时，破冰模式为挤压破坏和弯曲破坏的共同作用。为了研究船艏倾角改变对冰阻力的影响，计算了不同船艏倾角下的破冰阻力，如图13 所示。随着船艏倾角的增加，船艏的破冰方式中挤压破坏的比例逐渐升高，弯曲破坏的比例下降，船艏受到的破冰阻力逐渐增加。在现代破冰船舶设计中，船艏的外倾角增加和水线角与船艏倾角减小可以提高船舶的破冰性能，本文的模拟结果反映出这种随船艏倾角减小破冰阻力相应减小的效果。

图13 不同船艏倾角下破冰阻力比较Fig.13 Comparison of ice-breaking resistance at different bow angles

3.2 冰厚的影响

冰厚是海冰最重要的物理特性之一，不同厚度的冰层具有不同的力学性质，其抵抗破坏的能力也会不同。海冰条件对船舶冰阻力的影响，也以冰厚最为显著。为了分析冰厚对冰阻力的影响，计算了不同冰厚下的破冰阻力，图14 为不同冰厚下船艏破冰阻力比较。随着冰厚的增加，船艏的破冰阻力增加，并且速度越高，这种趋势越明显。因此冰区船在不同冰区航行时，应根据所处冰区的冰厚，在考虑效率和推进功耗的基础上，合理地选择最佳航行速度。