-
近 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}} $ 叠加计算,即$$ {{\mathit{\boldsymbol{F}}}}_{i}={{\mathit{\boldsymbol{F}}}}_{i}^{\mathrm{n}}+{{\mathit{\boldsymbol{F}}}}_{i}^{\mathrm{s}} $$ (1) 法向接触力
$ {{\mathit{\boldsymbol{F}}}}_{i}^{\mathrm{n}} $ 可直接由法向刚度系数$ {k}_{i}^{\mathrm{n}} $ 和法向重叠量$ {{\mathit{\boldsymbol{U}}}}^{\mathrm{n}} $ 计算得到:$$ {{\mathit{\boldsymbol{F}}}}_{i}^{n}=-{k}_{i}^{\mathrm{n}}{{\mathit{\boldsymbol{U}}}}^{\mathrm{n}} $$ (2) 式中:
${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}} $ 计算:$$ \Delta {{\mathit{\boldsymbol{F}}}}_{i}^{\mathrm{s}}=-{k}_{i}^{\mathrm{s}}\Delta {{\mathit{\boldsymbol{U}}}}_{i}^{\mathrm{s}} $$ (3) 式中:
$ {k}_{i}^{\mathrm{s}}{=\eta k}_{i}^{\mathrm{n}}, $ 为切向刚度系数,其中$ \eta $ 为切向刚度系数与法向系数之比;$ \Delta {{\mathit{\boldsymbol{U}}}}_{i}^{\mathrm{s}} $ 为切向重叠量的增量。同时,切向弹性接触力应满足Coulomb摩擦定律,即当单元间的切向力超过最大摩擦力时,切向力将变为最大摩擦力:
$$ {{\mathit{\boldsymbol{F}}}}_{i}^{\mathrm{s}}<{{\mathit{\boldsymbol{F}}}}_{i}^{\mathrm{s}}\left({{\mathit{\boldsymbol{F}}}}_{\mathrm{m}\mathrm{a}\mathrm{x}}^{\mathrm{s}}/\left|{{\mathit{\boldsymbol{F}}}}_{i}^{\mathrm{s}}\right|\right) $$ (4) 式中,
$ {{\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个黏结颗粒单元间设定一个弹性黏结圆盘,圆盘可以传递2个单元之间的力和力矩,即轴向力、剪力、弯矩和扭矩,并且力
$ {\bar{{\mathit{\boldsymbol{F}}}}}_{i} $ 和力矩$ {\bar{{\mathit{\boldsymbol{M}}}}}_{i} $ 都可以用各自的法向分量和切向分量来表示:$$ {\bar{{\mathit{\boldsymbol{F}}}}}_{i}={\bar{{\mathit{\boldsymbol{F}}}}}_{i}^{\mathrm{n}}+{\bar{{\mathit{\boldsymbol{F}}}}}_{i}^{\mathrm{s}} $$ (5) $$ {\bar{{\mathit{\boldsymbol{M}}}}}_{i}={\bar{{\mathit{\boldsymbol{M}}}}}_{i}^{\mathrm{n}}+{\bar{{\mathit{\boldsymbol{M}}}}}_{i}^{\mathrm{s}} $$ (6) 式中,
$ {\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}} $ 分别为力与力矩的法向分量和切向分量,其中法向分量可表示为:$$ {\bar{{\mathit{\boldsymbol{F}}}}}_{i}^{\mathrm{n}}={\bar{F}}^{\mathrm{n}}{{\mathit{\boldsymbol{n}}}}_{i} $$ (7) $$ {\bar{{\mathit{\boldsymbol{M}}}}}_{i}^{\mathrm{n}}={\bar{M}}^{\mathrm{n}}{{\mathit{\boldsymbol{n}}}}_{i} $$ (8) 式中:
$ {{\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 {\bar{{\mathit{\boldsymbol{F}}}}}_{i}^{\mathrm{n}}=-{k}_{i}^{\mathrm{n}}A\Delta {{\mathit{\boldsymbol{U}}}}^{\mathrm{n}} $$ (9) $$ \Delta {\bar{{\mathit{\boldsymbol{F}}}}}_{i}^{\mathrm{s}}=-{k}_{i}^{\mathrm{s}}A\Delta {{\mathit{\boldsymbol{U}}}}_{i}^{\mathrm{s}} $$ (10) 式中:
$ \Delta {{\mathit{\boldsymbol{U}}}}^{\mathrm{n}} $ 为法向位移增量;$ \Delta {{\mathit{\boldsymbol{U}}}}_{i}^{\mathrm{s}} $ 为切向位移增量;$ A $ 为黏结圆盘的横截面积。力矩增量可以表示为:
$$ \Delta {\bar{{\mathit{\boldsymbol{M}}}}}_{i}^{\mathrm{n}}=-{k}_{i}^{\mathrm{n}}J\Delta {{\mathit{\boldsymbol{\theta }}}}_{i}^{\mathrm{n}} $$ (11) $$ \Delta {\bar{{\mathit{\boldsymbol{M}}}}}_{i}^{\mathrm{s}}=-{k}_{i}^{\mathrm{s}}I\Delta {{\mathit{\boldsymbol{\theta }}}}_{i}^{\mathrm{s}} $$ (12) 式中,
$ \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可分别表示为:$$ \left\{ \begin{aligned} & A=\pi {\bar{R}}^{2}\\& J=\frac{1}{2}\pi {\bar{R}}^{4}\\& I=\frac{1}{4}\pi {\bar{R}}^{4} \end{aligned} \right. $$ (13) 式中,
$ \bar{R} $ 为黏结圆盘的半径。根据材料力学中的梁理论,其对应黏结单元圆盘周边上的最大拉应力
$ {\sigma }_{\mathrm{m}\mathrm{a}\mathrm{x}} $ 和剪切应力$ {\tau }_{\mathrm{m}\mathrm{a}\mathrm{x}} $ 可以表示为:$$ {\sigma }_{\mathrm{m}\mathrm{a}\mathrm{x}}=\frac{-{\bar{{\mathit{\boldsymbol{F}}}}}^{\mathrm{n}}}{A}+\frac{\left|{\bar{{\mathit{\boldsymbol{M}}}}}_{i}^{\mathrm{s}}\right|}{I}\bar{R} $$ (14) $$ {\tau }_{\mathrm{m}\mathrm{a}\mathrm{x}}=\frac{\left|{\bar{{\mathit{\boldsymbol{F}}}}}_{i}^{\mathrm{s}}\right|}{A}+\frac{\left|{\bar{{\mathit{\boldsymbol{M}}}}}^{\mathrm{n}}\right|}{J}\bar{R} $$ (15) 如果最大拉伸应力
$ {\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模拟中,将只考虑浮力与拖曳力对浮式结构的影响,而忽略海浪的作用。浮式结构在整体坐标轴方向平动和沿局部坐标轴转动的动力方程可以表示为:
$$ {{\mathit{\boldsymbol{F}}}}_{\mathrm{i}}+{{\mathit{\boldsymbol{F}}}}_{\mathrm{b}}+{{\mathit{\boldsymbol{F}}}}_{\mathrm{d}}+{{\mathit{\boldsymbol{F}}}}_{\mathrm{p}}=\left(M+A\right)\ddot{{\mathit{\boldsymbol{x}}}}+C\dot{{\mathit{\boldsymbol{x}}}}+K{\mathit{\boldsymbol{x}}} $$ (16) $$ {{\mathit{\boldsymbol{M}}}}_{\mathrm{i}}+{{\mathit{\boldsymbol{M}}}}_{\mathrm{b}}+{{\mathit{\boldsymbol{M}}}}_{\mathrm{d}}=I\dot{{\mathit{\boldsymbol{\omega }}}} $$ (17) 式中:
$ {{\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 }}}} $ 平台旋转运动的角加速度。浮力的计算采用阿基米德原理:
$$ {{\mathit{\boldsymbol{F}}}}_{\mathrm{b}}=-{\rho }_{\mathrm{w}}g{V}_{\mathrm{s}\mathrm{u}\mathrm{b}} $$ (18) 式中:
$ {\rho }_{\mathrm{w}} $ 为海水密度;$ g $ 为重力加速度;$ {V}_{\mathrm{s}\mathrm{u}\mathrm{b}} $ 为结构浸入水中的体积。浮式结构在水中所受的拖曳作用由2部分组成,分别为结构发生平动时的拖曳力
$ {{\mathit{\boldsymbol{F}}}}_{\mathrm{d}} $ 及旋转时的拖曳力矩$ {{\mathit{\boldsymbol{M}}}}_{\mathrm{d}} $ ,二者可表示为:$$ {{\mathit{\boldsymbol{F}}}}_{\mathrm{d}}={-C}_{\alpha }{\rho }_{\mathrm{w}}{V}_{\mathrm{s}\mathrm{u}\mathrm{b}}\left({{\mathit{\boldsymbol{v}}}}_{\mathrm{w}}-{{\mathit{\boldsymbol{v}}}}_{\mathrm{s}}\right) $$ (19) $$ {{\mathit{\boldsymbol{M}}}}_{\mathrm{d}}=-{C}_{\beta }{{\rho }_{\mathrm{w}}V}_{\mathrm{s}\mathrm{u}\mathrm{b}}{L}^{2}{{\mathit{\boldsymbol{\omega }}}}_{\mathrm{s}} $$ (20) 式中:
$ {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)所示。
张紧式锚泊系统采用线性弹簧模拟单一锚链的恢复力T:
$$ {\rm{T = }}{k_{\rm{m}}} \cdot \Delta {\rm{L}} $$ (21) 式中:
$ {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所示。表 1 Kulluk浮式平台与碎冰作用时的离散单元计算参数
Table 1. Computational parameters of DEM when Kulluk floating platform interacts with pack ice
参数 值 冰层大小(L×B)/m2 $ 300 \times 200 $ 海冰单元平均面积Si/m2 200 海冰密集度Ci/% 70 海冰密度ρi/(kg·m-3) 920 海水密度ρw/(kg·m-3) 1 035 冰速vi/(kg·m-1) 0.5 冰厚hi/m 0.5~2.0 颗粒间的弹性模量$ \bar{E} $/GPa 1 颗粒间法向(切向)黏结强度$ {\bar{\sigma }}_{\mathrm{n}}\left(\bar{\tau }\right) $/MPa 0.6 平台在冰厚为0.8 m的碎冰区作业时,纵、横向和垂向的冰载荷F时程曲线如图5所示。从该冰力时程曲线上可以看出,平台纵向和垂向的冰力更加连续且数值较大。其中,纵向冰力对平台漂移产生作用,垂向上的冰载荷是影响结构升沉运动的主要因素,而在横向方向上,因计算冰域与结构的对称性,导致其冰载荷均值接近0,大体成对称性分布。
在Kulluk浮式平台报告中,冰载荷数据为作用在结构上的总体冰载荷,即图5中3个方向冰载荷的合力。另外,报告中未有效地分清海冰冰况的密集度、冰型和冰块尺寸等信息,仅将其复杂的工况简化成了冰厚−冰载荷限值形式(式(22)和式(23))[5],其对应冰载荷的2个限值方程为:
$$ {F}_{\alpha }=87{h}_{\mathrm{i}}+91 $$ (22) $$ {F}_{\beta }=38{h}_{\mathrm{i}}+31 $$ (23) 式中,
$ {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]。
-
本文利用DEM建立了用于浮式结构冰载荷分析的数值预报模型,并通过对比Kulluk浮式平台实测数据,对该模型进行了验证性分析,得到了以下主要结论及未来需要研究的方向:
1)基于CUDA-C并行算法和Voronoi分割算法图形处理技术,可以较好地模拟冰区浮式结构-冰的相互作用过程,且该数值预报模型的有效性已通过了Kulluk浮式平台多年实测数据的检验,因此本文提出的数值模型可为冰区浮式平台冰载荷快速预报提供可行的研究手段。
2)数值分析结果及实测数据均表明,冰厚对浮式结构冰载荷的影响较大,将直接影响到平台的作业安全性,因此可通过海冰管理等方式对浮式结构的冰载荷进行控制。然而国内有关冰区工程管理的经验多集中于固定式平台结构的风险分析与生命周期管理,如何建立科学调配破冰船队、气象预报及作业单元的管理体系成为浮式平台安全运营的关键。
3)冰区浮式结构的锚泊系统模拟技术大多是采用线性弹簧来模拟张紧式锚泊系统,然而锚链自身的大变形和非线性特征也会较大地影响到浮式结构的运动响应特性,因此下一步将采用更加符合真实锚链力学特性的力学模型。
Numerical analysisof ice load on the floating platform in polar region based on high-performance discrete element method
-
摘要:
目的 冰载荷是极区浮式平台结构安全性分析及完整性管理的重要输入,建立快速有效的冰载荷预报技术可保障浮式平台分割的安全运营。 方法 基于CUDA-C并行处理技术,建立可用于海冰模拟的高性能离散元方法(DEM)。采用Voronoi分割算法生成平台作业时的碎冰域,其海冰由具有黏结-失效效应的球体单元组成,且黏结单元的破碎由失效准则控制。平台结构由三角单元构成,并考虑了其在浮力、锚链恢复力及冰载荷共同作用下的六自由度运动。 结果 计算了锚泊结构在海冰持续作用下的总体冰载荷,统计分析了冰厚的影响并与Kulluk实测结果进行了对比验证。 结论 极区浮式平台冰区运行时冰载荷受冰厚影响较大,应配套破冰服务以减弱其作业区域的冰情。 Abstract:Objectives Ice load is the important input for the floating platform in polar area in aspect of structural safety analysis and integrity management. And then the establishment of fast and effective forecasting technique of ice load can guarantee the safety operation of floating platform. Methods A high-performance discrete element method(DEM) used for the sea ice modeling based on CUDA-C parallel processing technique was established. The pack ice around the platform during the operation was generated by Voronoi algorithm and the ice sheet was constructed by the sphere elements considering the bond-break effect which was controlled by the failure criterion. The structure hull was constructed with triangular elements and treated as a six-degree of freedom system under the action of buoyancy, restoring force of mooring line, and ice load. The global ice load on the moored structure was calculated under the continuous action of ice. Results Moreover, the influence of ice thickness was studied statistically and the numerical results were validated with the field test results of Kulluk. Conclusions The influence of ice thickness on the ice load impacted on the floating platform during the operation in the ice-covered waters. And then the ice breaking service should be supplied to weaken the ice condition of the operation area. -
Key words:
- floating platform /
- ice load /
- discrete element method(DEM)) /
- parallel bond model
-
表 1 Kulluk浮式平台与碎冰作用时的离散单元计算参数
Table 1. Computational parameters of DEM when Kulluk floating platform interacts with pack ice
参数 值 冰层大小(L×B)/m2 $ 300 \times 200 $ 海冰单元平均面积Si/m2 200 海冰密集度Ci/% 70 海冰密度ρi/(kg·m-3) 920 海水密度ρw/(kg·m-3) 1 035 冰速vi/(kg·m-1) 0.5 冰厚hi/m 0.5~2.0 颗粒间的弹性模量 $ \bar{E} $ /GPa1 颗粒间法向(切向)黏结强度 $ {\bar{\sigma }}_{\mathrm{n}}\left(\bar{\tau }\right) $ /MPa0.6 -
[1] 张侠, 杨惠根, 王洛. 我国北极航道开拓的战略选择初探[J]. 极地研究, 2016, 28(2): 267–276. ZHANG X, YANG H G, WANG L. Strategic thinking on China's involvement in the development of Arctic sea routes[J]. Chinese Journal of Polar Research, 2016, 28(2): 267–276 (in Chinese). [2] EBINGER C K, ZAMBETAKIS E. The geopolitics of Arctic melt[J]. International Affairs, 2009, 85(6): 1215–1232. doi: 10.1111/j.1468-2346.2009.00858.x [3] YUE Q J, BI X J. Ice-induced Jacket structure vibrations in Bohai sea[J]. Journal of Cold Regions Engineering, 2000, 14(2): 81–92. doi: 10.1061/(ASCE)0887-381X(2000)14:2(81) [4] 李想, 李红霞, 黄一. 核电平台连接机构设计与运动响应分析[J]. 中国舰船研究, 2020, 15(1): 152–161. doi: 10.19693/j.issn.1673-3185.01786 LI X, LI H X, HUANG Y. Design of connecting mechanism and motion response analysis on nuclear power platform[J]. Chinese Journal of Ship Research, 2020, 15(1): 152–161 (in Chinese). doi: 10.19693/j.issn.1673-3185.01786 [5] WRIGHT B. Full scale experience with Kulluk Stationkeeping operations in pack ice (with reference to Grand Banks Developments)[R]. Canada: National Research Council, 2000: 52. [6] 董斌, 钱源, 李元泰, 等. 船体(平台)渤海冰区作业安全性分析[J]. 中国舰船研究, 2020, 15(1): 145–151, 169. doi: 10.19693/j.issn.1673-3185.01553 DONG B, QIAN Y, LI Y T, et al. Safety analysis of hull (platform) operation in Bohai sea ice area[J]. Chinese Journal of Ship Research, 2020, 15(1): 145–151, 169 (in Chinese). doi: 10.19693/j.issn.1673-3185.01553 [7] 王健伟, 邹早建. 基于非线性有限元法的船舶-冰层碰撞结构响应研究[J]. 振动与冲击, 2015, 34(23): 125–130. WANG J W, ZOU Z J. Ship's structural response during its collision with level ice based on nonlinear finite element method[J]. Journal of Vibration and Shock, 2015, 34(23): 125–130 (in Chinese). [8] 薛彦卓, 陆锡奎, 王庆, 等. 冰三点弯曲试验的近场动力学数值模拟[J]. 哈尔滨工程大学学报, 2018, 39(4): 607–613. XUE Y Z, LU X K, WANG Q, et al. Simulation of three-point bending test of ice based on peridynamic[J]. Journal of Harbin Engineering University, 2018, 39(4): 607–613 (in Chinese). [9] HOPKINS M A, THORNDIKE A S. Floe formation in Arctic sea ice[J]. Journal of Geophysical Research: Oceans, 2006, 111(C1): C11S23. [10] SUN S S, SHEN H H. Simulation of pancake ice load on a circular cylinder in a wave and current field[J]. Cold Regions Science and Technology, 2012, 78: 31–39. doi: 10.1016/j.coldregions.2012.02.003 [11] SHEN H H, HIBLER W D, LEPPÄRANTA M. On applying granular flow theory to a deforming broken ice field[J]. Acta Mechanica, 1986, 63(1/2/3/4): 143–160. [12] CUNDALL P A, STRACK O D L. A discrete numerical model for granular assemblies[J]. Géotechnique, 1979, 29(1): 47–65. [13] TIMCO G W, WEEKS W F. A review of the engineering properties of sea ice[J]. Cold Regions Science and Technology, 2010, 60(2): 107–129. doi: 10.1016/j.coldregions.2009.10.003 [14] POTYONDY D O, CUNDALL P A. A bonded-particle model for rock[J]. International Journal of Rock Mechanics and Mining Sciences, 2004, 41(8): 1329–1364. doi: 10.1016/j.ijrmms.2004.09.011 [15] Itasca Consulting Group. Inc. PFC3D (Particle Flow Code in 3 Dimensions) Version 4.0[M]. Minneapolis: ICG, 2008. [16] 阿比尔的, 郑颖人, 冯夏庭, 等. 平行黏结模型宏细观力学参数相关性研究[J]. 岩土力学, 2018, 39(4): 1289–1301. ABI E D, ZHENG Y R, FENG X T, et al. Relationship between particle micro and macro mechanical parameters of parallel-bond model[J]. Rock and Soil Mechanics, 2018, 39(4): 1289–1301 (in Chinese). [17] 龙雪, 宋础, 季顺迎, 等. 锥角对锥体结构抗冰性能影响的离散元分析[J]. 海洋工程, 2018, 36(6): 96–100. LONG X, SONG C, JI S Y, et al. Influence of cone angle on anti-icing performance of conical structure with numerical simulations of discrete element method[J]. The Ocean Engineering, 2018, 36(6): 96–100 (in Chinese). [18] 蔡柯, 季顺迎. 平整冰与船舶结构相互作用的离散元分析[J]. 船舶与海洋工程, 2016, 32(5): 5–14. CAI K, JI S Y. Analysis of interaction between level ice and ship hull based on discrete element method[J]. Naval Architecture and Ocean Engineering, 2016, 32(5): 5–14 (in Chinese). [19] 季顺迎, 狄少丞, 李正, 等. 海冰与直立结构相互作用的离散单元数值模拟[J]. 工程力学, 2013, 30(1): 463–469. doi: 10.6052/j.issn.1000-4750.2011.07.0417 JI S Y, DI S C, LI Z, et al. Discrete element modelling of interaction between sea ice and vertical offshore structures[J]. Engineering Mechanics, 2013, 30(1): 463–469 (in Chinese). doi: 10.6052/j.issn.1000-4750.2011.07.0417 [20] 狄少丞, 王庆, 薛彦卓, 等. 破冰船冰区操纵性能离散元分析[J]. 工程力学, 2018, 35(11): 249–256. doi: 10.6052/j.issn.1000-4750.2017.09.0698 DI S C, WANG Q, XUE Y Z, et al. Manoeuvrability analysis of an icebreaker based on discrete element method[J]. Engineering Mechanics, 2018, 35(11): 249–256 (in Chinese). doi: 10.6052/j.issn.1000-4750.2017.09.0698 [21] 王帅霖, 季顺迎. 锥体导管架海洋平台冰激振动的DEM-FEM耦合分析及高性能算法[J]. 海洋学报, 2017, 39(12): 98–108. WANG S L, JI S Y. Ice induced vibration of conical platform based on coupled DEM-FEM model with high efficiency algorithm[J]. Haiyang Xuebao, 2017, 39(12): 98–108 (in Chinese). [22] 朱红日, 季顺迎, 刘璐. 基于Voronoi切割算法的碎冰区构造及离散元分析[J]. 计算力学学报, 2019, 36(4): 454–463. doi: 10.7511/jslx20180410001 ZHU H R, JI S Y, LIU L. Construction of broken ice field with Voronoi tessellation algorithm and its DEM application[J]. Chinese Journal of Computational Mechanics, 2019, 36(4): 454–463 (in Chinese). doi: 10.7511/jslx20180410001 [23] POLOJÄRVI A, TUHKURI J. On modeling cohesive ridge keel punch through tests with a combined finite-discrete element method[J]. Cold Regions Science and Technology, 2013, 85: 191–205. doi: 10.1016/j.coldregions.2012.09.013 [24] 狄少丞. 基于GPU并行算法的海洋平台及船舶结构冰荷载的离散元分析[D]. 大连: 大连理工大学, 2015. DI S C. Discrete element simulation of ice load on offshore platform and ship hull based on GPU parallel algorithm[D]. Dalian: Dalian University of Technology, 2015 (in Chinese). [25] ZHOU L, SU B, RISKA K, et al. Numerical simulation of moored structure station keeping in level ice[J]. Cold Regions Science and Technology, 2012, 71: 54–66. doi: 10.1016/j.coldregions.2011.10.008 [26] LU W J, LUBBAD R, ALEKSEY S, et al. Parallel channels' fracturing mechanism during ice management operations. Part I: theory[J]. Cold Regions Science and Technology, 2018, 156: 102–116. doi: 10.1016/j.coldregions.2018.07.010 -