破冰船引航下极地船舶结构冰荷载的离散元分析

刘 璐，胡 冰 ，季顺迎

(1. 大连理工大学 工业装备结构分析国家重点实验室，辽宁 大连 116024；2. 上海交通大学 海洋工程国家重点实验室，上海 200240；3. 中国远洋海运集团有限公司，上海 200027)

北极航道大大减少了我国与欧洲的海运航程，同时避免了航行经过不稳定地区带来的安全性问题。采用极地船舶独立进行冰区航行会大大提高航运成本，而采用破冰船引航的方式则可有效降低成本[1]。船舶在冰区的运行面临破冰过程中的结构安全问题，对船舶设计和冰区航行保障的要求较高，需采用有效方法对冰荷载进行分析并总结不同航行条件下的冰荷载特点[2]。

目前，冰荷载的研究方法主要包括现场监测、模型试验、经验公式和数值模拟。国外开展船舶结构冰荷载的现场测量工作较早，积累了大量复杂冰况的冰荷载数据[3]。对船冰相互作用的物理过程进行系统总结和归纳并结合大量的实测数据结果，可建立经验化的直接计算公式[4]。近些年来，国内逐步开展了船舶结构冰荷载分析工作[5-6]，在冰荷载现场原型测试、数值模拟等领域开展了大量工作[7-8]。数值模拟方法具有较好的参数可控性、结果可重复性，经济成本也较低，还可对海冰破碎、冰块碰撞重叠等过程进行细化分析，其中离散元方法使用较为广泛[9]。离散元方法可对船冰相互作用过程中的海冰断裂机理进行深入分析，通过引入流固耦合模型可模拟船冰作用过程中流体对海冰和船舶结构的影响[10]。

基于闵可夫斯基和(Minkowski sum)的扩展多面体单元方法可准确描述海冰的非规则形态，对海冰与结构的碰撞过程中的冰载荷进行有效分析[11]。在海冰离散元方法中，海冰单元的运动和断裂完全决定于单元间的作用力传递和材料的自身参数，体现出更好的自适应性和健壮性。本文采用基于闵可夫斯基和原理构造扩展多面体单元，建立扩展多面体离散元方法。采用该方法计算海冰与船舶结构的相互作用过程，分析船舶结构的冰载荷并与经典的Lindqvist公式和Riska公式对比验证结果的可靠性。采用离散元方法模拟破冰船引航条件下的船冰作用过程，分析不同航速和船宽比条件下的冰载荷特点。

1 海冰的扩展多面体离散元方法

根据闵可夫斯基和原理构造扩展多面体单元，采用Voronoi切割算法构造扩展多面体形态的海冰单元。同时，发展了扩展多面体单元间的粘结失效模型，并采用非线性的接触力模型计算粘结失效后单元间接触力，从而建立海冰的扩展多面体离散元方法。

1.1 海冰的扩展多面体单元构造

对于空间中的两个几何体A和B，二者的闵可夫斯基和可以表示为：

采用2D-Voronoi切割算法划分平面单元，在此基础上生成具有一定厚度的海冰单元。Voronoi单元的几何形态可通过规则度参数进行控制[12]。如图2所示，首先将2D-Voronoi单元按照扩展半径进行“腐蚀”形成新的多边形，通过给定厚度生成具有一定厚度的三维多面体单元，最后通过闵可夫斯基和运算生成扩展多面体单元。

图 1 基于闵可夫斯基和的扩展多面体Fig. 1 Dilated polyhedron based on Minkowski sum

图 2 2D-Voronoi切割算法生成海冰单元Fig. 2 Sea ice element based on 2D-Voronoi tessellation algorithm

1.2 扩展多面体的粘结失效模型

考虑粘结作用在单元间的分布，本文采用在单元交界面上设置粘结点的方式粘结相邻且共面的单元。如图3所示，在粘结面上设置粘结点，粘结点之间的法向应变εn可写作：

式中：d为两粘结节点间的距离向量；n为交界面的法向，这里取为两单元粘结面的外法向之差，即n =unit(ni- nj)；Cij为两个粘结单元之间的特征长度，与刚体有限元中的定义相同[13]，即Cij= hi+ hj。考虑切向方向t = unit (d - (d·n)n)，则切向应变可写作：

图 3 相邻单元在交界面上的粘结作用Fig. 3 Bonding on interface between adjacent elements

根据三维条件下的弹性矩阵，交界面上两单元间的弹性应力σe可写作：

式中：E为材料弹性模量；ν为材料泊松比；ε为考虑粘结节点之间的法向和切向应变，即。

式中：A为交界面面积；n为一个粘结面上的粘结节点个数。将粘结力转换到单元质心上，结合接触力等其他受力，即可计算每个单元质心受到的合力。

针对每个粘结点受到的作用力，可采用不同的失效模型控制粘结失效。这里采用混合的断裂能和判断粘结失效，即同时考虑法向拉伸和剪切计算临界强度σn和σs，在粘结力超过临界应力强度后考虑刚度的软化过程。结合拉伸和剪切破坏并采用混合的断裂能模型计算断裂能，从而根据断裂能计算粘结点间的临界变形[14]。若粘结点间的变形超过该临界变形，则该粘结失效。本文发展的粘结失效准则更全面地考虑了脆性材料的破碎特性，可对海冰的破坏过程进行更精确模拟。

1.3 扩展多面体的接触模型

扩展多面体间的接触可采用基于近似包络函数的接触搜索算法计算接触点和接触重叠量[11]。这里采用非线性的接触力模型简化计算扩展多面体之间的接触力[15]。其中，法向接触力Fn包括法向弹性力和法向黏滞力，即：

式中：kn为法向接触刚度，且；E*和R*分别为等效弹性模量等效单元半径[16]；为两单元的法向接触重叠量；为法向相对速度；cn为法向阻尼系数，且有，mAB为等效质量[16]，ζn为无量纲阻尼系数，e为回弹系数[17]。

切向接触力与切向重叠量和摩擦力有关，一般切向力满足库伦摩擦定律。这里采用简化的方式避免与摩擦力的频繁比较和判断。其中，切向弹性力写作：

式中：ct为切向阻尼系数，这里取[18]；为切向相对速度。该切向力计算方法可避免在每个时间步中与摩擦力进行对比，提高了计算效率。

2 破冰船辅助航行中冰荷载的离散元模拟及验证

采用扩展多面体离散元方法模拟船舶结构在平整冰区的航行过程，并通过船舶结构冰荷载的Lindqvist公式和Riska公式验证模拟冰荷载的准确性。通过模拟破冰船辅助航行时货船的冰荷载特点，分析不同参数条件下冰荷载的变化趋势。

2.1 离散元模拟参数

图4(a)为“雪龙号”全比例船型，船体结构的长度为167.0 m，型宽22.6 m，吃水9.0 m。将船体结构表面划分为三角片元用于海冰单元与船体结构之间的接触搜索，模拟中船体结构视为不变形的刚体结构。平整冰区的3个方向采用固定边界约束，船舶以速度vs在平整冰区航行。

图 4 雪龙号破冰船与平整冰作用的离散元模拟Fig. 4 Simulation sketch of the interaction between icebreaker Xuelong and level ice

扩展多面体离散元方法可合理分析海冰与结构相互作用中的弯曲破坏模式，在与ISO199006的冰力规范公式对比中也体现出了良好的准确度。模拟中主要的计算参数列于表1[19]，考虑海冰受到水的拖曳作用，拖曳力和拖曳力矩可写作[20]：

表 1 船舶结构破冰过程模拟的主要计算参数Tab. 1 Main parameters in the simulation of icebreaking process of ship hull

2.2 船舶结构冰荷载模拟结果及验证

采用扩展多面体离散元方法模拟单个船体在冰区航行过程，将稳定段冰载荷时程的均值作为冰阻力，并通过经验公式验证分析本文数值方法的准确性。在模拟中，平整冰区的尺寸L×W为600 m × 150 m。

图5是船舶结构上3个方向的冰荷载时程曲线。可见，船体刚进入冰区时x和z方向上的冰荷载经过一段时间的上升，之后冰荷载趋于稳定且在一定水平上持续波动。由于船体结构破冰过程中海冰会经历由船艏至船舯型线不断收缩变小的过程，海冰会在船肩处剧烈地由两侧挤压船体，因此，y方向上船体结构也会受到较大的冰荷载作用。这里主要考虑船体是否具有向前的破冰能力，即重点考虑x方向上的冰荷载。在后面的分析中，冰荷载和冰阻力均指x方向的结构受力。

图 5 3个方向的船体冰荷载时程曲线Fig. 5 Time history of ice loads on ship hull in three directions

采用稳定阶段的冰载荷均值作为冰阻力，将冰阻力与船舶结构冰阻力的Lindqvist公式和Riska公式计算结果对比。Lindqvist公式将船舶结构冰阻力分为船艏处海冰挤压冰力、海冰弯曲破坏冰力和破碎海冰的浸没阻力三部分，在冰区船舶设计中使用较为广泛[4]。Lindqvist中每部分冰力均采用理想化的模型总结计算公式，因此其计算的冰阻力较为保守。Riska公式是在Lindqvist基础上总结了波罗的海实船测量数据而得出的经验公式，在冰区船舶的冰荷载计算中具有代表性[21]。图5(a)中的虚线为对应工况下Lindqvist公式计算的冰阻力。图6是采用扩展多面体离散元方法计算的冰阻力随冰厚的变化，及其与Lindqvist公式和Riska公式对比结果。从图6可见，三者计算的冰阻力均随冰厚增大而增大，且具有明显的线性增加趋势；在具体数值上，离散元结果基本介于Lindqvist公式和Riska公式之间。因此，通过Lindqvist公式和Riska公式可以充分说明本文扩展多面体离散元方法计算船舶结构冰荷载具有良好的准确性和可靠性。

图 6 船体冰荷载的离散元计算结果与经验公式对比Fig. 6 Comparison between DEM result and empirical formulas of ice loads on ship hull

图 7 货船结构Fig. 7 Hull structure of cargo ship

2.3 破冰船引航下船舶结构冰荷载的离散元分析

采用扩展多面体离散元方法模拟破冰船辅助航行中破冰船和商船的冰荷载，即破冰船在前、货船紧随其后。平整冰区的尺寸为1 000 m × 120 m。货轮长度为200 m，型宽27.8 m，吃水12.3 m（图7）。

图8(a)是无引航条件下破冰船和货船的冰载荷时程曲线，虚线是稳定阶段冰载荷的均值。从图8(a)可以看出，货船的冰荷载则与破冰船差别较大，进入冰区后冰荷载会出现类似脉冲形态的波动，但是其整体水平较破冰船要小。由于货船较破冰船要长，其冰荷载依然比破冰船小很多。图8(b)是无引航条件下货船的冰荷载时程曲线，虚线是稳定阶段的冰阻力。无引航条件下货船冰阻力为5.9 MN，比引航条件下的货船冰阻力要大。但由于货船不具备破冰能力，其艏柱倾角较小且水线处存在较为尖锐的艏部结构，容易使海冰发生挤压破碎，这与破冰船艏部大倾角造成海冰弯曲破坏不同。因此，无引航条件下货船冰荷载存在近100 MN的峰值载荷，对货船的结构安全构成严重威胁。因此，在有破冰船引航条件下，货船的冰阻力则可明显降低，且不会存在较大的峰值荷载。

图 8 有无破冰船引航条件下船舶结构冰载荷时程曲线Fig. 8 Time history of ice loads on ship hulls with and without icebreaker pilotage

采用扩展多面体离散元方法分析不同船速和船宽比条件下两船的冰阻力特点。船宽比是货船与破冰船宽度的比，即w2/w1。图9(a)是不同船速条件下两船的冰阻力对比，模拟中两船的宽度比为1.0。从图9(a)可以看出，若采用专业破冰船在破冰能力较弱的船舶之前进入冰区进行破冰作业，即使货船的结构尺寸比破冰船大，其受到的冰荷载阻力也会小很多，能够大大缓解货船在冰区航行的阻力。图9(b)是不同的船宽比条件下两船的冰阻力，模拟中船速为1.5 m/s。从图9(b)可以看出，破冰船的破冰阻力并没有受到货船的影响，基本保持在相同的水平。在货船宽度明显小于破冰船时，即w2/w1= 0.5时，货船冰阻力比破冰船小很多。随着货船宽度增大，特别是在w2/w1> 1时，货船的冰阻力明显增大。

图 9 不同工况下破冰船和货船的冰阻力Fig. 9 Ice resistance on icebreaker and cargo ship under different conditions

3 结 语

介绍了基于闵可夫斯基和的扩展多面体单元方法，并将该方法应用于船舶结构与海冰的相互作用模拟，分析船舶结构上的冰荷载。通过离散元计算的船舶冰荷载与Lindqvist公式和Riska公式对比验证了离散元方法，得出本文的计算结果在数值和趋势上具有良好的准确性。采用本文方法模拟破冰船引航条件下货船的冰荷载，分析了不同船宽和船速条件下破冰船和货船的冰阻力特点。计算结果表明，破冰船引航可以有效降低货船的冰荷载，但在不同船宽比情况下，货船冰荷载有明显的差异。本文方法在海冰的离散元模拟中引入了非规则的计算单元，并在单元间的粘结失效中考虑了混合断裂能和刚度软化过程，从而使海冰的破坏更符合真实的物理过程。