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

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

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

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

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

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

图  1  单元间法向与切向接触模型

Figure 1.  Normal and tangential contact model between elements

法向接触力$ {{\mathit{\boldsymbol{F}}}}_{i}^{\mathrm{n}} $可直接由法向刚度系数$ {k}_{i}^{\mathrm{n}} $和法向重叠量$ {{\mathit{\boldsymbol{U}}}}^{\mathrm{n}} $计算得到：

式中：${k}_{i}^{\mathrm{n}}=\dfrac{{R}^{\mathrm{A}}\cdot {R}^{\mathrm{B}}}{{R}^{\mathrm{A}}+{R}^{\mathrm{B}}}G\pi$，为法向刚度系数其中$ G $为颗粒的弹性模量，$ {R}^{\mathrm{A}},{R}^{\mathrm{B}} $分别为2个接触颗粒A和B的半径。

切向弹性接触力由其增量$ \Delta {{\mathit{\boldsymbol{F}}}}_{i}^{\mathrm{s}} $计算：

式中：$ {k}_{i}^{\mathrm{s}}{=\eta k}_{i}^{\mathrm{n}}， $为切向刚度系数，其中$ \eta $为切向刚度系数与法向系数之比；$ \Delta {{\mathit{\boldsymbol{U}}}}_{i}^{\mathrm{s}} $为切向重叠量的增量。

同时，切向弹性接触力应满足Coulomb摩擦定律，即当单元间的切向力超过最大摩擦力时，切向力将变为最大摩擦力：

式中，$ {{\mathit{\boldsymbol{F}}}}_{\mathrm{m}\mathrm{ax}}^{\mathrm{s}}$为最大摩擦力，${{\mathit{\boldsymbol{F}}}}_{\mathrm{max}}^{\mathrm{s}}=\mu \left|{{\mathit{\boldsymbol{F}}}}_{i}^{\mathrm{n}}\right| $，其中$ \mu $为单元间摩擦系数。

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

图  2  模拟海冰冻结效应的平行黏结模型

Figure 2.  Parallel bond model for simulating the freezing effect of sea ice

所谓平行黏结模型，是指在2个黏结颗粒单元间设定一个弹性黏结圆盘，圆盘可以传递2个单元之间的力和力矩，即轴向力、剪力、弯矩和扭矩，并且力$ {\bar{{\mathit{\boldsymbol{F}}}}}_{i} $和力矩$ {\bar{{\mathit{\boldsymbol{M}}}}}_{i} $都可以用各自的法向分量和切向分量来表示：

式中，$ {\bar{{\mathit{\boldsymbol{F}}}}}_{i}^{\mathrm{n}},{\bar{{\mathit{\boldsymbol{M}}}}}_{i}^{\mathrm{n}},{\bar{{\mathit{\boldsymbol{F}}}}}_{i}^{\mathrm{s}},{\bar{{\mathit{\boldsymbol{M}}}}}_{i}^{\mathrm{s}} $分别为力与力矩的法向分量和切向分量，其中法向分量可表示为：

式中：$ {{\mathit{\boldsymbol{n}}}}_{i} $为黏结单元法向量；$ {\bar{F}}^{\mathrm{n}} $和$ {\bar{M}}^{\mathrm{n}} $分别为$ {\bar{{\mathit{\boldsymbol{F}}}}}_{i}^{\mathrm{n}} $和$ {\bar{{\mathit{\boldsymbol{M}}}}}_{i}^{\mathrm{n}} $对应的标量。

当把平行黏结施加于两颗粒上时，颗粒之间相对位置和相对转角的改变会在接触点处产生力增量和力矩增量，其中力的增量可表示为：

式中：$ \Delta {{\mathit{\boldsymbol{U}}}}^{\mathrm{n}} $为法向位移增量；$ \Delta {{\mathit{\boldsymbol{U}}}}_{i}^{\mathrm{s}} $为切向位移增量；$ A $为黏结圆盘的横截面积。

力矩增量可以表示为：

式中，$ \Delta {{\mathit{\boldsymbol{\theta }}}}_{i}^{\mathrm{n}},\Delta {{\mathit{\boldsymbol{\theta }}}}_{i}^{\mathrm{s}} $分别为法向角度增量和切向角度增量，可由$ \Delta {{\mathit{\boldsymbol{\theta }}}}_{i}=\left({{\mathit{\boldsymbol{\omega }}}}_{i}^{B}-{{\mathit{\boldsymbol{\omega }}}}_{i}^{A}\right)\Delta t $计算，其中$ {{\mathit{\boldsymbol{\omega }}}}_{i}^{A} $和$ {{\mathit{\boldsymbol{\omega }}}}_{i}^{\mathrm{B}} $分别为颗粒A和B的旋转速度；J为圆盘横截面积的极惯性矩；I为圆盘横截面积的惯性矩。

黏结圆盘的横截面积$ A $、极惯性矩$ J $和惯性矩I可分别表示为：

式中，$ \bar{R} $为黏结圆盘的半径。

根据材料力学中的梁理论，其对应黏结单元圆盘周边上的最大拉应力$ {\sigma }_{\mathrm{m}\mathrm{a}\mathrm{x}} $和剪切应力$ {\tau }_{\mathrm{m}\mathrm{a}\mathrm{x}} $可以表示为：

如果最大拉伸应力$ {\sigma }_{\mathrm{m}\mathrm{a}\mathrm{x}} $超过了给定的法向黏结强度，或最大剪切应力$ {\tau }_{\mathrm{m}\mathrm{a}\mathrm{x}} $超过了给定的切向黏结强度，颗粒间的黏结将发生破坏，单元之间将失去黏结效应，进而表现为海冰材料中产生裂纹，最终导致海冰破碎。

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

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

图  3  Kulluk浮式平台

Figure 3.  Kulluk floating platform

式中：$ {{\mathit{\boldsymbol{F}}}}_{\mathrm{i}},{{\mathit{\boldsymbol{F}}}}_{\mathrm{b}},{{\mathit{\boldsymbol{F}}}}_{\mathrm{d}},{{\mathit{\boldsymbol{F}}}}_{\mathrm{p}} $分别为结构的冰力、海水浮力、海流拖曳力及锚链或推进器的作用力；$ {{\mathit{\boldsymbol{M}}}}_{\mathrm{i}}, {{\mathit{\boldsymbol{M}}}}_{\mathrm{b}},{{\mathit{\boldsymbol{M}}}}_{\mathrm{d}} $分别为海冰对结构产生的力矩、海水浮力矩和拖曳力产生的力矩；M，A，I，C，K分别为浮式结构的质量、附加质量、转动惯量、阻尼及刚度系数；$ \ddot{{\mathit{\boldsymbol{x}}}},\dot{{\mathit{\boldsymbol{x}}}},{\mathit{\boldsymbol{x}}} $分别为平台平移的加速度、速度和位移；$ \dot{{\mathit{\boldsymbol{\omega }}}} $平台旋转运动的角加速度。

浮力的计算采用阿基米德原理：

式中：$ {\rho }_{\mathrm{w}} $为海水密度；$ g $为重力加速度；$ {V}_{\mathrm{s}\mathrm{u}\mathrm{b}} $为结构浸入水中的体积。

浮式结构在水中所受的拖曳作用由2部分组成，分别为结构发生平动时的拖曳力$ {{\mathit{\boldsymbol{F}}}}_{\mathrm{d}} $及旋转时的拖曳力矩$ {{\mathit{\boldsymbol{M}}}}_{\mathrm{d}} $，二者可表示为：

式中：$ {C}_{\alpha },{C}_{\beta } $分别为对应于拖曳力与拖曳力矩的拖曳系数；$ {{\mathit{\boldsymbol{v}}}}_{\mathrm{w}} $为水流速度；$ {{\mathit{\boldsymbol{v}}}}_{\mathrm{s}} $为结构平移的速度；$ {{\mathit{\boldsymbol{\omega }}}}_{\mathrm{s}} $为结构转动的角速度；$ L $为浮体长度。

浮式结构的锚泊系统可以采用张紧式锚链将锚泊点与浮式平台连接起来^[25]，其具体形式如图4(a)所示。锚链系统采用12根锚链将平台进行约束，其锚链具体的分布形式如图4(b)所示。

图  4  锚链布置图

Figure 4.  Layout of the mooring lines

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

式中：$ {k}_{\mathrm{m}}=1.03\times {10}^{3}\;\mathrm{k}\mathrm{N}/\mathrm{m} $，为锚链的刚度；$ \Delta L $为锚泊点与平台连接点之间的相对变形量。

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

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

图  5  3个方向的船体冰载荷时程曲线

Figure 5.  Time history curves of ice loads in three directions

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

式中，$ {F}_{\alpha } $和$ {F}_{\beta } $分别为2种海冰情况下的限值；$ {h}_{\mathrm{i}} $为海冰冰厚。

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

图  6  DEM模拟中破碎冰与Kulluk浮式平台间的相互作用

Figure 6.  Interaction between pack ice and Kulluk floating platform simulated with DEM

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

图  7  DEM结果与实测数据的对比^[26]

Figure 7.  Comparison between the DEM results and field test data^[26]

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

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

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

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