冰区航行中船体结构冰压力分布特性的离散元分析

刘 璐，曹 晶，张志刚，季顺迎

（1.大连理工大学工业装备结构分析国家重点实验室，辽宁大连116024；2.中国船级社上海规范研究所，上海200135）

0 引 言

随着全球气候变暖，人类在极区的自然资源勘探考察和开发逐渐增多。其中，北极航道的开辟大大减少了海运航程，也避免了航经不稳定地区带来的安全性因素。然而，船舶在冰区的运行面临破冰过程中的结构安全问题，急需研究冰载荷在船体结构表面的分布情况，从而识别船体结构的危险区域并分析结构的冰激疲劳问题，从而确定冰载荷监测的最优位置[1]。

自上世纪80 年代起，国外在极区开展了大量的船体结构冰压力现场测量工作，积累了大量一年冰和多年冰与船体结构相互作用的相关数据[2]。为研究局部冰压力的分布特性并从海冰的物理力学性能角度分析冰压力成因，国外也开展了大量的冰压力现场试验[3]。然而，现场的试验测量经济成本较高，且现场试验结果较难获得重复性验证，难以获得更加广泛的规律。因此，借助于海冰物理力学性质的相关研究成果实现对冰压力的数值预测，可为深入探究冰压力的形成机理和冰区船舶的工程设计提供必要的手段[4]。对冰-船相互作用的物理过程进行系统的总结和归纳，形成海冰断裂过程的物理模型并确定模型参数，同时结合大量的实测数据结果，通过简化的海冰断裂模型发展半解析半数值方法预测船体结构冰载荷和冰压力分布情况[5]。而采用有限元、黏聚单元、离散元等方法分析海冰与结构作用过程则更具有工程实用性，可深入分析不同海冰工况下海冰与结构相互作用过程中的断裂机理[6]。

近些年来，国内逐步开展了破冰船冰载荷的现场监测和数值计算工作[7]。在海洋工程中使用较多的压力盒测试法可以直观地测得结构冰载荷，但其成本较高，不适合船-冰接触面积较大的船体结构。根据国外的相关经验，采用时域法中的影响系数矩阵法可通过船体结构上的应变测试结果反演冰载荷[8]。在寒区海洋工程领域，也开展了相应的海冰模型试验工作[9]。在破冰船冰载荷的数值预测领域，国内开展了较多的工作[10]，在船-冰作用模型、冰压力预测等方面发展迅速。值得注意的是，在数值模拟中大多采用商业软件如ABAQUS、LS-DYNA，或将高度经验化的数学模型引入数值方法中，不能较好地满足数值预测的广泛适用性和对复杂条件的兼容性要求。

离散单元法最早用于离散介质的动力过程模拟，也可用于碎冰的堆积成型过程模拟[11]。随着其粘结模型的不断发展，该方法也可用于连续介质的破坏过程，在岩土工程、地质体演变等领域应用广泛[12]。离散元法在海冰领域的应用较早，在冰脊成型、海冰与结构相互作用等方面发展较快[13]。采用扩展圆盘单元可以模拟船体在碎冰区的航行过程，可对碎冰区船舶结构冰载荷的影响因素进行定量化的分析[14]。通过引入基于GPU（Graphics Processing Unit）的并行计算可实现大规模离散元计算，提高对海冰断裂过程的分析精度和效率[15]。在海冰离散元中，海冰单元的断裂完全取决于外力的传递过程和材料自身参数，体现出了更好的自适应性和健壮性。不同于前述数值模型专注于船体结构的受力分析，海冰离散元方法同时可分析海冰的动力和断裂过程，包括堆积、阻塞等情况。因此，离散元方法是一种能够满足广泛需求的船-冰相互作用分析及冰载荷预测的方法。

本文通过带有粘结-破碎模型的海冰离散元方法模拟破冰船与平整冰作用过程，同时采用基于GPU 的大规模并行技术大幅提高计算效率和规模。针对船体结构在平整冰区直行破冰过程，分析了船体结构上的总体冰载荷及局部冰压力，并通过规范验证了冰压力模拟的可靠性。

1 破冰船冰区航行的海冰离散元方法

在平整冰与船体结构的相互作用过程中，将海冰离散为均匀大小和质量的球体单元。海冰单元的受力考虑单元之间粘结力、接触力以及海水对单元的浮力和拖曳力[16]，这里主要阐述粘结-破碎模型的实施方式。同时，通过对船体结构表面节点力的统计合理建立船体结构表面冰压力的评估方法。

1.1 海冰单元的粘结-破碎模型

采用平行粘结模型对离散球体单元进行冻结[17]。平行粘结模型假定相邻的两个球体单元之间通过圆形截面梁进行粘结，该模型可以同时传递粘结力和力矩，如图1所示。

图1 球体单元之间的平行粘结模型Fig.1 Parallel bond model between spherical elements

粘结力考虑因法向和切向弹性变形而引起的弹性力模型，这里采用与弹性接触力模型统一的力模型计算，根据梁模型计算法向和切向的力矩。将力和力矩分为拉力、剪力、扭矩和弯矩，即

式中，F 和M 分别代表力和力矩；下标b 代表粘结，上标n 和s 分别代表法向和切向分量。在粘结模型中考虑梁的拉伸、弯曲和扭转作用，两个球形单元之间的最大法向应力和剪切应力为

式中，A、J 和I分别为梁的截面积、极惯性矩和惯性矩，且有A = πR2，J = πR4/2，I = πR4/4，R 为球体单元的半径。假设拉力为正，压力为负，采用Mohr-Coulomb理论判断粘结的断裂失效，

式中：τb即为切向粘结强度；C 为材料粘聚力强度；σ 是法向应力；μ表示内摩擦系数，μ = tanθ，θ 表示曲线倾角，如图2 所示。粘结的失效准则可以表示为

图2 海冰离散元中的Mohr-Coulomb失效准则Fig.2 Mohr-Coulomb criterion in DEM for sea ice

式中，σb为法向拉伸强度。图2中实线外侧即为断裂区域，即上式表示的失效准则。

在海冰离散元的实际应用中，如何根据已知海冰的强度等力学性质合理选择离散元参数是数值模拟准确性的关键。通过单轴压缩和三点弯曲试验的离散元模拟，研究颗粒单元尺寸、粘结强度和内摩擦系数对离散元模拟结果的影响，进而可确定海冰宏观强度与以上三者之间关系[18]。海冰的弯曲强度与颗粒单元的粘结强度紧密相关，可通过对海冰的三点弯曲试验确定海冰宏观强度σf与离散元粘结失效参数σb之间的关系，

式中，h 为海冰厚度，D 为颗粒单元直径。另外，一般在具体的模拟中，令μ = 0.2 且σb= C[15]。以上关系式与渤海现场试验中测得的海冰弯曲强度对比验证结果良好[19]，可用于渤海海冰力学性质的模拟分析，是海冰离散元参数选取的重要依据。

1.2 船体结构的冰压力计算

将船体的外壳结构划分为若干三角形单元。通过海冰单元对船体结构的碰撞力插值到三角形单元的三个节点上，并考虑三角形单元的面积，即可得到每个节点上的冰压力。

由于离散元方法的积分步长非常小，通常在10-9～10-5s 范围内，直接统计在某个时间步上的瞬时冰压力显然是不合理的。这里通过统计平均的方式，将一段时间内的冰压力平均值作为这段时间的冰压力，那么第i个统计时间步某一节点上的冰压力可表示为

式中，pj是第j 个时间步的冰压力，r 和s 分别代表该段时间的起始和结束时间步数。式中表示的冰压力即为瞬时冰压力，它体现了船体结构在较短时间内的冰压力情况。但是船体结构的瞬时冰压力不能很好地体现航行过程中的统计特性。为了克服瞬时冰压力在表述上的缺点，更好地体现船体结构上冰压力较大且作用频率较高的区域，即高压和易疲劳区域，这里引入累计最大冰压力和累计平均冰压力。第N个统计时间步某一节点上的累计最大冰压力和累计平均冰压力可表示为

从船体结构上的所有节点来看，累计最大冰压力体现了船体结构上的高压危险位置，该区域可能会产生结构强度破坏；累计平均冰压力体现了船体结构上冰压力较大且作用频率较高的位置，会产生较为持续且强度较高的冰压力作用，容易产生材料的疲劳强度破坏。

2 破冰船在平整冰区的离散元模拟

采用海冰离散元模拟船体结构在平整冰区破冰时与海冰的相互作用过程，分析船舶破冰过程中冰压力在船体结构表面的分布形式。船体结构模型采用某极地科学考察船，计算中将船体结构视为刚体，船体长、宽和高分别为167 m、22.6 m 和13.5 m，吃水深度为9 m。图3 所示为该科学考察船船体表面划分为三角形的结构形式，模拟中三角单元的平均尺寸为1.5 m。

计算中主要的计算参数列于表1 中。其中，海冰的弯曲强度是海冰物理力学参数的重要表征，直接影响到海冰的破坏模式，进而影响海冰与结构物相互作用过程中的冰载荷[19]。在海冰与结构相互作用的离散元模拟中，需考虑将海冰的宏观强度转化为离散元中球体单元间的微观强度，确保模拟中海冰的宏观强度与实际保持一致，该过程中还需综合考虑单元的尺寸效应[20]。海冰受到海流的拖曳力模型参考Sun和Shen（2012）[16]的相关工作。

图3 科学考察船的船体结构三角单元划分Fig.3 Triangular meshes of ship hull for the scientific research ship

表1 船体结构破冰过程模拟的离散元参数Tab.1 DEM parameters in the simulation of icebreaking process of ship hull

2.1 船体结构破冰过程的离散元分析

连续的直行破冰是破冰船冰区航行中的主要破冰操作模式，这里采用离散元方法分析船体结构直行破冰过程中的总体冰载荷和局部冰压力，主要计算参数列于表2中。在模拟中，平整冰区计算域的长×宽为600 m×200 m。计算域边界和海冰单元之间采用弹簧约束，弹簧刚度与单元之间的接触刚度相同。这里参考ISO19906，取极地海冰的弯曲强度σf=0.7 MPa[21]。当海冰厚度为1.5 m 时，两层球体单元的总单元数约为41 万，且球体单元粒径为0.826 m。海冰弯曲强度为0.7 MPa 时，根据式（7）算得单元粘结强度为1.315 MPa。采用基于GPU 的并行技术，破冰船的直行破冰过程模拟如图4 所示。从图中可以看出，在船体结构的作用下海冰发生断裂，在冰面上形成边缘不规则的断裂从而形成开阔水道。

表2 船体结构直行破冰模拟中的相关计算参数Tab.2 Main parameters in the simulation of ship icebreaking in straight line

图4 船体结构破冰过程的离散元模拟Fig.4 DEM simulation of the icebreaking process of ship hull

图5 船体结构破冰过程中的海冰破坏模式Fig.5 Failure mode of sea ice during icebreaking of ship hull

图5 是船体结构破冰过程中的海冰破坏模式对比。从图中可看出，离散元模拟中船体结构两侧出现了“锯齿”状连续的断裂，这种破坏模式是典型的弯曲破坏。在实际船舶走航的海冰断裂监测中，如雪龙号和芬兰某型冰区船，也出现了类似的破坏模式。但是实船监测中海冰还呈现“环状”的破坏，并不是规整的“锯齿”状。在离散元模拟中船体结构是刚体，只有x方向的刚性平动；而在实船破冰过程中，由于船体两侧非同时的海冰结构和破坏，船体还会产生横摇，从而对船体两侧的海冰造成挤压作用。因此，在实船破冰过程中海冰还会发生剪切破坏，其“环状”的裂纹是剪切和弯曲破坏等复杂作用下的结果[22]。

船体结构破冰过程中的总体冰载荷如图6 所示。从图中可以看出，x 和z 方向的冰载荷具有类似的特点，即在较高的水平线上下波动，这种波动体现了海冰的断裂过程，同时说明了船体结构运动过程中海冰的持续作用。y 方向的冰载荷在0 附近上下波动，且会出现较大的峰值载荷，说明了船舶在破冰过程中海冰会从两侧对结构产生挤压，而这种挤压作用会对结构安全性造成严重隐患。因此船舶两侧的冰载荷影响还需要对局部的冰压力进行深入的分析。

图6 船体结构上的总体冰载荷时程Fig.6 Time history of ice loads on ship hull

船体结构破冰过程中的局部冰压力如图7 所示。从图中可以看出，船舶破冰过程中最大的冰压力主要集中在船肩处，该部位应是船体结构安全的重点监测部位。由于破冰部位主要为船艏柱，该部位会与海冰发生较多的接触，所以累计平均冰压力集中于此，也说明船艏柱是冰激结构疲劳分析重点关注的部位。

图7 船体结构上的冰压力分布Fig.7 Distribution of ice pressure on ship hull

2.2 船体结构冰压力的对比验证

国际船级社组织（International Association of Classification Societies,IACS）对特定工况下船体结构的冰压力给出了规范计算方法[23]。该计算针对船体结构与大块浮冰的自由碰撞过程，通过能量守恒原理并结合一定的经验处理，可总结该碰撞过程的理论模型。在理论模型中根据碰撞点处的船型参数给出了等效的压力板大小，根据该压力板上的冰载荷即可求出碰撞点上的冰压力[24]。通过离散元计算了船体结构与大块浮冰的碰撞过程，分析了碰撞点上的冰压力并与IACS 的规范进行对比，从而验证本文离散元方法的合理性。采用IACS 冰级PC5的相应参数作为离散元计算参数，主要参数列于表3中，其他相关参数参照表1。其中，船速为2.25 m/s，海冰厚度为3.0 m，海冰弯曲强度为1.0 MPa[25]。根据式(7)可由海冰厚度和弯曲强度确定粘结失效强度的具体数值。船体结构与大块浮冰自由碰撞的离散元模型如图8所示。设置船体结构与海冰的碰撞点为距船艏柱5 m、15 m、25 m和35 m位置，并分别根据规范推荐方法计算对应碰撞点上的冰压力。

表3 船体结构冰压力验证分析中的相关计算参数Tab.3 Main parameters in the validation of ice pressure on ship hull

图8 冰压力规范验证的离散元模型Fig.8 Sketch of DEM simulation for the validation with IACS standard

船体结构与大块浮冰碰撞的离散元模拟过程如图9 所示。从图中可以看出，船艏部位与海冰发生碰撞造成海冰的挤压破碎。随着船体结构的运动，海冰不断发生破碎且破碎后的海冰呈粉末状浮在水面上。图10 为四个算例中碰撞点处的冰压力与规范对比情况。其中，x 为碰撞点距船艏柱的距离，alpha和beta是碰撞点处的船型参数，w和b为根据规范计算的压力板宽度和高度，P为离散元模拟的冰压力，Pstd为IACS规范值。

图9 船体结构与大块浮冰碰撞的离散元模拟Fig.9 DEM simulation of collision between ship hull and ice floe

图10 碰撞点处船体结构的冰压力Fig.10 Ice pressure on ship hull on the collision position

离散元模拟结果与规范进行对比的情况列于表4 中。离散元的计算结果与规范值相比存在上下误差，压力值保持在相同量级范围内，误差范围在6.7%～18.1%之间。因此，本文的离散元方法对船体结构的冰压力可进行准确可靠的分析，具有良好的工程适用性。

表4 离散元模拟结果与IACS标准对比Tab.4 Comparison between DEM simulation and IACS standard

3 结 语

本文基于平行粘结模型建立了基于球体单元的海冰离散元方法，并通过Mohr-Coulomb 准则实现了单元间的断裂，从而可对海冰与结构相互作用的破碎过程进行模拟；采用海冰离散元方法模拟了科学考察船的破冰过程，研究了船体结构上的总体冰载荷时程；同时，结合船舶破冰的特点，对累计最大冰压力和累计平均冰压力进行了有效的分析；为验证本文计算结果的可靠性，通过IACS 规范对本文离散元方法计算的船体结构冰压力进行了校核。结果表明，本文离散元方法的冰压力计算误差范围在6.7%～18.1%之间，充分说明了该方法的合理性。