基于广义接触模型的离散元方法及其对船体冰载荷的分析

孔 帅，季顺迎

(大连理工大学 工业装备结构分析国家重点实验室，辽宁 大连 116023)

海冰广泛存在于地球的南北极及高寒地区，是影响北极航道开发和通航的主要环境要素[1]。船舶在极地航行中会受到来自于海冰和海浪等多方面的动态激励，其中海冰对船舶振动的影响最大[2]；另外船冰之间的碰撞，如冰山撞击，会对船舶结构安全产生极大的安全隐患[3]。

通过数值方法提前确定冰区航行的冰载荷可为船体重点区域安全维护提供技术支持，同时为船舶冰区航行提供预警标准以及相应安全的航行操作指导。在模拟海冰与船体作用时，根据研究需要已开发出众多的仿真技术，如采用LS-Dyna动力学软件平整冰与船体之间的作用过程[4]、根据弹性梁理论推导出的海冰弯曲断裂模式模型[5]、依据现场观测数据建立的二维平面内环形裂纹破坏模型[6]等不同数值方法。

离散元方法由Cundall等[7]提出，最早是应用于岩石、谷物等散体材料的动力特性分析方法。由于浮冰大量存在于海冰初始生成及融化阶段，因此可采用离散元方法对浮冰进行模拟分析，如李紫麟等[8]将莲叶冰离散成圆盘单元，用于分析船体结构冰与莲叶冰之间的相互作用，Sun等[9]模拟了浮冰、波浪和海流共同对直立结构导管架平台的作用。另外，在与船体结构之间的作用中，海冰的破碎过程属于典型的脆性破坏[10]。脆性材料的破碎过程可采用具有黏结-失效效应的离散元模型进行模拟。如模拟海冰与斜坡结构物的多面体单元Hybrid黏结失效准则可表征出渐进性过程[11]，但相比多面体单元，球体单元因其高效的接触搜索判断及较为准确地接触力计算而更便于应用。Potyondy和Cundall[12]基于球形颗粒提出了平行黏结模型，该模型已被PFC3D软件采用并广泛应用于岩石等脆性材料的裂纹生成研究中[13-14]。海冰在结构物作用下的破坏过程也可采用平行黏结模型，如龙雪等[15]将其应用到锥体海洋结构的尺寸参数对平台抗冰性能影响分析中。但现有模拟海冰破碎的平行黏结模型仅考虑了黏结-失效单元整体的力学行为，缺少对球体黏结单元失效渐进性过程的关注。

船舶结构冰区航行过程中的冰阻力特性分析有助于研究船体结构整体运动和航行推进器[16-17]。为更加有效地模拟船舶结构与海冰的作用过程，采用三维广义接触模型对海冰的破碎过程进行模拟，并依据海冰强度试验对离散元中关键细观参数的影响进行了研究。最后将其扩展应用到冰区船舶冰载荷分析中。

1 基于广义接触模型模拟海冰破碎的离散元方法

广义接触模型(3D Generalized contact model，简称GCM-3D)是由GCM-2D发展而来[18-19]，能够很好地应用于岩石材料破碎分析。GCM-3D的黏结单元由局部黏结点组成，根据局部黏结点的失效情况确定整个黏结单元是否失效，进而判断整体破坏是否发生。

1.1 黏结单元的广义接触模型

GCM-3D局部黏结点分布在黏结单元的中间面上，其典型分布形式如图1所示，颗粒单元间的圆点即为局部黏结点。图1(a)为局部黏结点在面(t,n)的分布形式，图1(b)为局部黏结点在面(t,s)的分布形式。

图1 GCM-3D模型的局部黏结点分布

在颗粒的相互作用过程中，考虑单元间相对速度和弹性变形而引起的作用力，颗粒单元之间的接触可采用弹簧-阻尼器-滑块的唯象模型[20]，并且单元间的接触力可解耦为法向分量和切向分量，解耦后法向与切向的接触模型如图2所示，单元间的接触力由各个局部黏结点对之间作用力叠加而成。

图2 黏结单元的广义接触模型

局部黏结点对之间的法向和切向作用力可由线弹性模型计算，即：

(1)

(2)

(3)

(4)

1.2 黏结单元局部黏结点的刚度和强度

GCM-3D模型颗粒单元间刚度由局部黏结点对之间刚度叠加，即：

(5)

(6)

局部黏结点未失效时法向和切向能承受的最大力由式(7)和式(8)定义，局部黏结点失效时的切向力由式(9)定义：

(7)

(8)

Fs=μFn

(9)

综上可知，局部黏结点的接触力本构模型采用如图3所示的本构模型。其中，法向最大接触力由法向黏结强度决定，剪切力服从Mohr-Coulomb准则[21]。当法向接触力或切向接触力达到式(7)或式(8)的数值时，该局部黏结点对便失去黏结效应。当局部接触点发生失效时，局部接触点只考虑接触点之间互相挤压下的法向接触力和切向接触力，其中切向接触力由单元间的摩擦系数确定(式(9))。

图3 局部黏结点对间的本构模型

2 广义接触黏结模型细观参数对海冰力学性质的影响

图4 海冰与船舶结构作用时的破坏模式

GCM-3D模型中影响海冰宏观力学性能的关键细观参数有单元间法向黏结强度、切向黏结强度和接触摩擦系数等。船舶结构在海冰发生持续作用时，海冰的破坏形式复杂，由Lindqvist[22]提出船体结构与海冰相互作用时海冰破坏主要以弯曲破坏为主，局部区域有挤压破坏和剪切破坏。图4为“雪龙号”极地科学考察船艏部区域的海冰破坏图像，分析图像资料可知，弯曲断裂和挤压破坏为海冰的主要破坏模式。因此文中离散元参数影响研究依据海冰弯曲强度试验和单轴压缩强度试验进行。

2.1 海冰单轴压缩和三点弯曲试验模拟

为进行海冰单轴压缩强度和弯曲强度试验的模拟，采用Ji等[23]在渤海海域所做的海冰强度时的模型尺寸。海冰的单轴压缩试验中，采用尺寸为70 mm × 70 mm × 175 mm的长方体冰试样，压头压缩速率为2 mm/s，单元数目为7 800，粒径大小为5.328×10-3m，并在竖直方向施加载荷；海冰弯曲强度的获取采用三点弯曲试验方法，试验中海冰为75 mm × 75 mm × 700 mm的长方体试样，压头压缩速率为3.5 mm/s，单元数目为23 940，粒径大小为6.033×10-3m，并在试样上表面施加竖直方向载荷。离散元模拟所需的参数如表1所示。

表1 离散元模拟的主要计算参数

海冰试样在压头持续的作用下会发生局部破坏最终导致试样整体破坏。图5和图6为采用离散元模拟的不同加载时刻海冰在压头持续作用下颗粒间作用力的分布。其中，单轴压缩试验在未发生破坏时颗粒间作用力均匀分布，最终模型整体因发生剪切错动而破坏，与脆性材料力学的单轴压缩试验现象一致；三点弯曲试验中在未发生破坏时其颗粒间作用力较大区域集中于压头作用位置的上下表面，模型在压头持续作用下最终在压头位置处破坏。

图5 海冰单轴压缩试验的离散元数值模拟

图6 海冰三点弯曲试验的离散元数值模拟

海冰单轴压缩和三点弯曲数值模拟中压头施加在加载面上力的变化曲线分别如图7(a)、(b)所示。将海冰试样发生破坏时的最大压应力用于确定单轴压缩强度和弯曲强度，其强度分别为4.38 MPa和1.40 MPa。

图7 单轴压缩和三点弯曲试验模拟中压力变化曲线

2.2 局部接触黏结点对间的强度参数影响

图和对海冰单轴压缩强度的影响

图和对海冰三点弯曲强度的影响

图 和单轴压缩强度与三点弯曲强度比值的影响

根据Timco等[24]对北极海冰强度的统计分析可知，当年冰的单轴压缩强度可以设定在0.5～5.0 MPa之间，弯曲强度在1.0 MPa左右。一般情况下，海冰单轴压缩强度为弯曲强度的2～6倍，数值模拟中的海冰单轴压缩强度与弯曲强度之间比值要与实际情况相吻合[23]。综合本节分析可知，北极海冰的单轴压缩强度、三点弯曲海冰强度可由局部接触黏结对间的强度参数线性插值获得；同时为保证线性插值的稳定性，模拟中接触摩擦系数确定为μc=0.2。

3 船舶结构的冰阻力特性分析

海冰物理力学性质的复杂性、船舶结构的复杂性、海冰-船体之间作用的差异化导致了船舶结构的冰阻力尚未有准确的计算公式，目前的冰阻力计算公式绝大多数是针对平整冰区稳定航行的半经验公式。这些半经验公式结合了船体破冰过程的理论分析和现场实测数据，并将大量相关的影响参数引入到公式中，可以较为准确地反映出不同参数对冰阻力的影响，对破冰船的设计选型具有一定的指导意义。但经验公式不可对破冰过程中的局部冰载荷分布(如线载荷分布、冰压分布)、冰力时程特征、冰载荷极值分布等特征进行描述。因此可通过经验公式验证数值方法的有效性，进而利用数值方法全方位描述冰载荷特征的性能。

3.1 Riska平整冰区船舶冰阻力半经验公式

Riska公式是基于室内模型试验及波罗的海的大量实船数据分析得出的[25]，并被芬兰-瑞典冰级规范作为冰区船舶螺旋桨推进力推算公式[26]，其公式如下：

F=C1+C2v

(10)

式中：v为船速。该公式考虑了船体形状因素(如船体进水角、船艏长度和船舯长度等)、船体航行速度及航行过程中海冰厚度的影响，C1和C2的取值如公式所示：

(11)

(12)

式中：α为水线处进水角，T为冰区水线，B为型宽，L为水线处全船长度，Lpar为船舯长度，Lbow为艏部长度，hi为冰厚，f1～f4、g1～g3等经验系数由试验测定，其取值分别为0.23 kN/m3、4.58 kN/m3、1.47 kN/m3、0.29 kN/m3、18.9 kN/(m·s-1·m1.5)、0.67 kN/(m·s-1·m1.5)、1.55 kN/(m·s-1·m1.5)。

3.2 Lindqvist平整冰区船舶冰阻力半经验公式

Lindqvist冰力公式也是目前使用较多的破冰船冰阻力计算经验公式[22]，主要考虑了海冰在弯曲断裂模式下的冰载荷。涉及到的参数包括主尺度、船型、冰厚、摩擦以及海冰的弯曲强度。该方法中冰阻力包括破冰阻力、海冰弯曲阻力和浸没阻力三部分，其中，破冰阻力Rc为：

(13)

弯曲阻力Rb为：

(14)

式中：E为海冰弹性模量，υ为泊松比，B为船体型宽，ρw为海水密度，g为重力加速度。

浸没阻力Rs为：

Rs=(ρw-ρi)ghtotBK

(15)

(16)

式中：ρi为海冰密度，htot为冰与雪的厚度和，T为冰区水线，L为船体长度。

考虑速度影响下的总冰阻力Ri：

(17)

3.3 基于广义接触模型的平整冰区船舶冰阻力分析

球体离散元是一种简化的计算模型，无法完全体现海冰各向异性的特征，如柱状晶体生长结构，但可通过对如厚度方向、海冰平面内方向设置有差别的黏结强度模拟海冰强度各向异性的特征。现有的离散元模拟中的黏结强度虽未考虑各向异性特征，但仍可模拟出结构物在海冰作用下的冰载荷特征，包括时程曲线、峰值特征及分布规律[15, 20]。同时，离散元模拟海冰与结构物作用过程中，海冰的粒径越小越能表征出海冰破碎的特征，但会导致计算资源大幅度增加，通常选择在厚度方向布置两层球体单元即可较为精确地描述出冰载荷特征[15]。

为分析平整冰区船体冰阻力特性，建立如图11左侧区域的平整冰区，区域大小为240 m×95 m；图11右侧为厚度方向上球形颗粒的分布形式，通过网格搜索后的相邻颗粒之间初始都采用广义接触模型进行黏结。

图11 由球体颗粒黏结组成的平整冰区示意图

表2为采用离散元模拟时破冰船在平整冰区稳定航行时的参数，其中海冰颗粒之间的法向、切向黏结强度为0.5 MPa。经2.2节中细观参数-宏观参数对应关系可知，海冰弯曲强度和单轴压缩强度分别为0.93 MPa和2.79 MPa；同时，海冰弹性模量为4.41 GPa。以“雪龙号”极地科学考察船为分析对象，船体结构由3 344个三角形单元组成，并以1.54 m/s的恒定船速在平整冰区航行。表2也列出了Riska公式[25]和Lindqvist公式[22]中的计算参数，其海冰与结构间的摩擦系数采用了Lindqvist建议的数值区间[0.1,0.16]，文中取μ=0.15。

表2 离散元模拟中的主要计算参数

图12展示了破冰船在0.8 m厚的平整冰区运行时不同角度的破冰状况。从图12(a)和(b)可以看出海冰会在艏部的持续作用下产生弯曲破坏的断裂模式，并有海冰在船肩位置处因挤压的作用发生挤压破坏并发生堆积。因此船舶重点防护和加固的区域应设置在船艏及船肩位置处。

图12 破冰船在平整冰区航行不同角度的展示

破冰船在破冰过程中的冰力主要由两部分组成：一部分来自船艏部位冲击破冰时冰的破碎和弯曲；另一部分来自于船肩、船舯等部位与海冰之间的摩擦作用，并且摩擦力与有效的接触长度有关，长度越长，摩擦作用越明显。图13为破冰船在冰区航行时的冰阻力变化曲线，可以看出船舶结构的冰载荷因海冰不断的发生弯曲断裂破坏，是一个典型的具有周期性的载荷；另外随着船体不断的驶入冰区，其冰阻力也在逐步提升。当全船进入冰后(x=150 m)，其冰阻力也趋于稳定。选取x=150～230 m的稳定破冰阶段冰载荷画于图13的右上角(其固定边界仅在附近5 m范围内即x=235 m才会有影响，因而固定边界对x=150～230 m的冰载荷没有影响)，均值为1 367.8 kN，其均值可用于确定冰区稳定航行时螺旋桨的推进力，即：

Fp=Fice+Fwater

(18)

式中：Fp为螺旋桨冰区稳定航行时所需推进力，Fice为冰区稳定航行时的冰阻力，Fwater为开阔水域航行时的阻力。

图13 冰厚为0.8 m时船舶冰阻力时程曲线

图14(a)为冰厚在0.8～1.2 m时稳定破冰阶段的冰阻力与经验公式确定的冰载荷对比，图14(b)为三种船速(1.03 m/s、1.54 m/s和2.06 m/s)下1.2 m厚冰的冰阻力与经验公式数值对比。可以看出数值仿真的结果整体上与Lindqvist的结果更加接近，与Riska的结果相差较多。Riska公式数值较小主要是其参考波罗的海的破冰船的测量数据，但波罗的海的冰相对较薄；且该公式对冰载荷与船速、冰厚之间进行了线性简化，使其数值普遍小于Lindqvist冰阻力公式[27]。数值仿真的结果从整体上符合船舶结构冰载荷的规律，为确定冰区航行螺旋桨推进力和冰载荷分布规律提供了可行的计算方法。

图14 数值仿真中冰载荷均值与经验公式的对比

4 结 语

采用广义接触模型建立了船舶结构冰载荷计算的离散元方法，同时也分析了不同冰厚和航速下冰载荷的变化规律。为分析广义接触模型中关键细观参数对冰载荷模拟的影响，对海冰单轴压缩试验和三点弯曲试验进行了数值模拟，确定了强度参数对不同海冰强度的影响，并采用线性插值的方式获得北极海冰的强度参数。该模型中的强度参数可以用于根据现场冰情、船舶航行特征确定相应的离散元计算参数。通过对极地科学考察船“雪龙号”在平整冰区航行时的冰载荷分析可知，该数值模型可以较好地对应冰载荷经验公式。该数值方法可用于冰区船舶螺旋桨功率配比及船舶在冰区航行冰阻力特性研究。为提升该数值方法的适用性，后续将会增加更多冰况、航行方式分析船舶结构的整体冰阻力、局部冰载荷分布特性。