基于高性能离散元方法的极区浮式平台冰载荷数值分析

孔帅，季顺迎，季少鹏，王迎晖，刚旭皓

1 中国船舶科学研究中心, 江苏 无锡 214082 2 大连理工大学 工业装备结构分析国家重点实验室, 辽宁 大连 116024

0 引 言

近 40 年来，受全球气候急剧变暖的影响，极地海冰覆盖面积持续减少，引起了各国开发利用北极自然资源的热情[1-2]。目前，我国在固定式海油钻井平台的运营方面已具有较为完备的设计建造基础[3]。移动式结构中的浮式结构适用于水深较深的海域，更加便于在极区的移动作业[4-6]，但我国在冰区浮式结构运营方面尚处起步阶段。因此，建立一套科学、有效的冰载荷预报技术是冰区浮式结构优化设计和冰激疲劳损伤特性分析的关键。

近年来，基于单元法的冰载荷预报分析技术已受到各国学者的关注，例如有限单元法[7]、近场动力学[8]和离散元法（discrete element method，DEM）[5,9-11]。离散元法由Cundall 等[12]提出，适于模拟散体材料在准静态或动态载荷激励下的运动、变形及破坏过程。寒区海域中的海冰呈现出很强的离散分布特性，离散元模型已被成功应用于碎冰动力过程[5,9-10]及海冰流变学的研究[11]。

为进一步精细化模拟浮式结构−冰相互作用的过程，离散元法需要在海冰的破碎分析、海冰拓扑形态表征及计算效率等方面有所突破。在海冰与浮式结构的相互作用过程中，海冰的破碎属于典型的脆性破坏过程[13]。Potyondy 等[14-22]基于球形单元提出了平行黏结模型，该模型通过在接触单元之间建立梁单元模型来传递单元之间的力和力矩，并根据梁单元上的最大应力是否达到预设强度来判断材料是否破碎。平行黏结模型的工程应用性较强，适用于海冰等材料的脆性破坏分析[17-22]。Voronoi 图是几何学中重要的图形概念，该图可以有效利用限定空间内特征点的分布形式来表征自然状态下物质的离散特性。该方法已被朱红日等[22]应用于极区浮碎冰的模拟生成中，具有生成效率高、复现能力强的优点。在离散元数值模型中，接触判断及单元搜索占据了较多的计算资源，若采用串行处理技术的传统离散元数值分析方法，浮式结构冰载荷的分析效率将会极低。龙雪等[17,20-22]运用并行处理化技术，建立了高性能离散元算法并用于分析导管架平台及船体结构的冰载荷。

为了分析浮式结构冰区作业时的冰载荷，本文将首先基于高性能离散元分析模型及采用平行黏结模型对海冰的破碎过程进行模拟，然后，再模拟构造冰情，分析平台与海冰的作用过程，最后，利用实测数据对离散元数值结果进行验证。

1 模拟海冰动力过程的离散元法

基于单元模拟技术的离散元法可以较为便捷地从细观角度描述海冰单元间的动力行为，且可有效构造复杂的冰型，例如碎冰[22]和冰脊[23]等。对于海冰的破碎过程，可通过引入平行黏结模型进行模拟[14,17-18]。图形处理器(graphics processing unit，GPU)拥有数量可观的核且每个核又包含多个线程，近年来，英伟达公司的CUDA-C 产品利用GPU集群处理技术为计算模型的并行化处理提供了新的研究思路。离散元算法中，接触判断和单元搜索计算资源占据了绝大部分的计算资源，为此，可通过结合单元邻居链表算法及背景网格搜索法进行并行化加速计算处理，计算效率是串行程序的数十倍[24]，在复杂作业环境模拟、多介质场耦合分析及超大计算域分析等方面具有较强的计算优势。

1.1 海冰单元接触模型

在离散元法中，通常把散体看作具有一定形状和质量颗粒单元的集合，单元之间的运动方程相对独立，但可以发生接触并产生相互作用力。在颗粒的相互作用过程中，考虑因单元之间相对速度和弹性变形而引起的作用力。颗粒单元之间的碰撞过程可采用弹簧−阻尼器−滑块的唯象模型进行模拟[19]，并且单元之间的接触力可解耦为法向分量和切向分量，解耦后的法向与切向接触模型如图1 所示。因此，第i个接触对中单元之间因接触而产生的力Fi可 由法向接触力Fni和切向接触力Fsi叠加计算，即

图1 单元间法向与切向接触模型Fig. 1 Normal and tangential contact model between elements

1.2 模拟海冰破碎的平行黏结模型

海冰的破碎过程直接影响船体结构冰载荷的频率和数值，海冰在船体的持续挤压过程中将变为尺寸更小的冰块。通过分析观测得到的北极海冰资料，可知其破碎冰多呈多边形的几何形状。图2所示为采用Voronoi 算法中的种子点分布特性控制海冰几何规则度而生成的碎冰域，该算法可有效控制碎冰场中的海冰面积及其分布特征。考虑破碎冰的再破碎，本文采用平行黏结模型模拟冰层的冻结作用。海冰内部不仅存在力的传递，还存在力矩的传递，因此，对海冰进行离散元分析时常采用平行黏结单元[17-22]。平行黏结是将2 个球体黏结在一起，其黏结单元不仅可以传递力，还可以传递力矩，如图2 左下角所示。

图2 模拟海冰冻结效应的平行黏结模型Fig. 2 Parallel bond model for simulating the freezing effect of sea ice

2 锚泊系统作用下浮式结构的冰载荷分析

冰区采油平台或者采油船等浮式结构通常采用锚泊系统在冰区作业，同时，周边会有破冰船协助破冰。Kulluk 石油钻井浮式平台隶属于壳牌公司，其曾长期服役于阿拉斯加等北极地区，在冰区积累了大量的冰情及冰载荷实测资料[5]。

2.1 浮式结构运动求解

Kulluk 浮式平台的整体结构为双圆锥体，向下的倾斜角为31.4°，该倾斜角的设计使得海冰在此易发生弯曲破碎而降低冰载荷；其上甲板处直径为81 m，吃水为11.5 m，水线处直径70 m，排水量28 000 t。平台现场工作图及对应的有限元模型如图3 所示，平台结构外表面由三角形单元构成。冰区作业中，浮式结构会受到海水浮力、海流拖曳力及海浪的作用，本文在采用离散元法模拟浮式结构与海冰相互作用时，将只考虑浮力与拖曳力对浮式结构的影响，而忽略海浪的作用。浮式结构在整体坐标轴方向平动和沿局部坐标轴转动的动力方程可以表示为：

图3 Kulluk 浮式平台Fig. 3 Kulluk floating platform

浮式结构的锚泊系统可以采用张紧式锚链将锚泊点与浮式平台连接起来[25]，具体形式如图4(a)所示。锚链系统采用12 根锚链将平台进行约束，锚链的具体分布形式如图4(b) 所示。图中：D，d分别为Kulluk 平台双圆锥结构上、下表面的外直径；L为锚链长度；h为水深；xG，yG分别为横、纵坐标轴；θ 为锚链间夹角。

图4 锚链布置图Fig. 4 Layout of the mooring lines

张紧式锚泊系统采用线性弹簧模拟单一锚链的恢复力T：

2.2 对比分析实测数据

根据在Kulluk 浮式平台现场实测的数据，可知平台周围主要以破碎冰为主。利用Voronoi算法生成浮冰，浮冰区域面积为 300 m×200 m，破碎冰平均面积为 200 m2，破碎冰的运动由速度控制在0.5 m/s 的固定边界推进。离散单元其他模拟参数如表1 所示。

表1 Kulluk 浮式平台与碎冰作用时的离散单元计算参数Table 1 Computational parameters of DEM when Kulluk floating platform interacts with pack ice

Kulluk 浮式平台在冰厚为0.8 m 的碎冰区作业时，纵向、横向和垂向的冰载荷F时程曲线如图5 所示。从该冰力时程曲线上可以看出，平台的纵向、垂向冰力更连续且数值较大。其中，纵向冰力对平台漂移产生作用，垂向冰载荷是影响结构升沉运动的主要因素，而在横向，由于计算冰域与结构的对称性，导致其冰载荷均值接近0，大体成对称性分布。

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

在Kulluk 浮式平台报告中，冰载荷数据为作用于结构上的总体冰载荷，即图5 中3 个方向冰载荷的合力。另外，报告中未有效分清海冰冰况的密集度、冰型和冰块尺寸等信息，仅将其复杂的工况简化成了冰厚−冰载荷限值形式（式（22）和式（23）)[5]，其对应冰载荷的2 个限值方程为：

限值Fα对应于平台周围海冰清理不良时的冰载荷上限值，此时海冰冰况严重，包括冰脊、重叠冰等复杂冰型；Fβ对应于平台周围海冰清理较好时的冰载荷上限值，此时海冰冰况较轻，平台附近会有破冰船协助破冰，有很多开阔的水域。图6示出了采用离散元法模拟的1.1 m 厚碎冰与Kulluk浮式平台相互作用的过程，其中图6(a)显示的是碎冰在平台前发生堆积，而平台两侧的碎冰则继续向前运动并在平台后方形成一段开阔的水域(图6(b))。从中可以看出，采用离散元法模拟的海冰冰情介于实测报告中所提的清理良好与清理不良这2 种工况之间。

图6 离散元法模拟中破碎冰与Kulluk 浮式平台间的相互作用Fig. 6 Interaction between pack ice and Kulluk floating platform simulated with DEM

图7 所示为不同冰厚(hi=0.5～2.0 m) 下平台总体冰载荷时程均值（Fi）与实测报告中数据（以吨计位）的对比。从中可以看出，由离散元法模拟得到的结果位于实测冰载荷范围内，说明离散元冰载荷模型可以较好地适用于浮式结构的冰载荷分析。离散元数值结果显示，在冰厚较小时黏结单元较易失效，大块海冰碎裂成小尺寸冰块。此时，海冰不宜在平台迎冰面发生堆积转而向平台两侧滑移，进而导致冰力降低，可见，平台冰载荷受冰厚的影响较大。海冰管理作业可以及时、有效地减弱平台周边的冰情，因此在实际平台作业过程中，可以根据观测的冰厚信息及时调配破冰船进行海冰管理作业，破碎其周边海冰以提升结构安全性[5,26]。

图7 离散元法模拟结果与实测数据的对比[26]Fig. 7 Comparison between the DEM results and field test data[26]

3 结 论

本文利用离散元法建立了用于浮式结构冰载荷分析的数值预报模型，并通过对比Kulluk 浮式平台的实测数据，对该模型进行了验证性分析，得到以下主要结论及未来需要研究的方向：

1） 基于CUDA-C 并行算法和Voronoi 分割算法图形处理技术，可以较好地模拟冰区浮式结构−冰的相互作用过程，且该数值预报模型的有效性已通过了Kulluk 浮式平台多年实测数据的检验，因此本文提出的数值模型可以为冰区浮式平台冰载荷快速预报提供可行的研究手段。

2） 数值分析结果及实测数据均表明，冰厚对浮式结构冰载荷的影响较大，将直接影响平台的作业安全性，因此，可通过海冰管理等方式对浮式结构的冰载荷进行控制。然而，国内有关冰区工程管理的经验多集中于固定式平台结构的风险分析与生命周期管理，如何建立科学调配破冰船队、气象预报及作业单元的管理体系成为浮式平台安全运营的关键。

3） 冰区浮式结构的锚泊系统模拟技术大多采用线性弹簧来模拟张紧式锚泊系统，然而锚链自身的大变形和非线性特征也会对浮式结构运动响应特性带来较大的影响，因此，下一步将采用更加符合真实锚链力学特性的力学模型。