冰脊压剪试验及其对直立结构冰载荷的离散元分析1)

朱红日 季顺迎

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

冰脊广泛分布于渤海、波罗的海、极区等海域，通常由大面积浮冰相互挤压、剪切进而发生破碎、堆叠而成[1-2]。典型的冰脊一般包括上部的脊帆、水面附近的冻结层及水下的龙骨三个部分，厚度可达5∼30 m[3]。冰脊作用在海洋结构造成的冰载荷远大于平整冰，很多情况下是船舶与海洋工程结构的极限设计载荷。因此，对冰脊的力学性质及其对海洋工程结构冰载荷的研究具有重要的工程意义。

现场测量与室内模型试验一直是研究冰脊的主要手段。在波弗特海的Molikpaq 平台[4]、加拿大的联邦大桥[5]及波的尼亚湾的Norstrmsgrund 灯塔[6]等海洋工程结构上开展了大量的冰脊载荷现场测量；此外，在现场及实验室内也开展了一系列的直剪[7]、双轴压缩[8]、压剪[7,9-10]等力学性质试验以研究冰脊特别是龙骨部分的力学特性。文献[11-12] 基于不同的冰脊破坏模式发展了相应的冰载荷计算模型，但目前不同的方法存在很大离散性。

近年来，有限元和离散元(discrete element method，DEM)等数值方法能计算复杂工况的优势日益显现，成为冰脊载荷研究的重要途径。Heinonen[10]和 Serr[7]等建立了冰脊龙骨的本构模型，采用有限元法对冰脊的力学特性进行了系统的数值计算。Hopkins[13]采用二维块体离散元方法模拟了冰脊的形成过程，发现冰脊在形成过程中存在巨大的能量耗散；Polojrvi 等[14-15]发展了三维块体离散元方法并计算了部分冻结及未冻结条件下冰脊龙骨的压剪试验，分析了加载速度对压剪试验结果的影响；Yulmetov 等[16]利用由球体组合而成的冰块单元计算了冰脊与锥体结构的相互作用过程以及海冰堆积现象。在工程计算中，冰脊内部冰块的尺寸一般远小于冰脊尺寸，冰块的形状对冰脊整体力学特性的影响可以忽略，由此可将冰块简化为球体单元从而兼顾计算效率及准确性。

为此，本文采用具有粘结−破坏功能的球体单元构造冰脊模型，通过计算现场压剪试验过程以分析计算模型的合理性及主要计算参数对冰脊力学特性的影响；在此基础上进一步计算冰脊与圆形直立结构的相互作用，并通过与相关规范进行对比以验证该离散元方法的可靠性。

1 冰脊的离散元方法

粘结球体离散单元方法由Cundall 和Potyondy于2004 年提出[17]，目前已广泛用于岩石、混凝土、玻璃等破坏特性的数值计算，其在模拟材料脆性破坏及力学行为方面具有显著优势。冰脊由大量的碎冰块堆积冻结而成，可采用该离散元方法进行构造和力学性质的数值计算。

1.1 冰脊结构的离散元构造

冰脊自上而下分为脊帆、冻结层、龙骨三个部分，通常将脊帆及龙骨理想化为梯形结构，如图1(a)和图1(b) 所示[18]。冻结层部分由于已完全冻结，可采用平整冰的离散元模型由相同粒径的球体单元按密六方排列构造[19]。龙骨与脊帆由部分冻结或未冻结的冰块堆积而成，其内部结构较为复杂，可由不同粒径的球体单元随机排列构造，如图 1(c) 所示。本文采用一种球体几何随机排列方法构造龙骨和脊帆的离散元模型[20]。在冰脊离散元方法中，根据冰块间的冻结状态确定单元接触对间是否存在粘结力。接触力的计算将单元间的法向力简化为弹簧和阻尼器，并将切向力简化为弹簧、阻尼器和滑动器，根据颗粒间的相对运动状态进行计算。当颗粒接触对处于粘结状态时，采用平行粘结模型[21]模拟冰块间冻结作用，如图2 所示。

图1 冰脊的几何结构及其离散元模型

图2 冰脊球体单元间的平行粘结模型

在平行粘结模型中，粘结单元间存在一个虚拟的弹性圆盘，由此传递两个单元间的轴向力Fn、剪力Fs、弯矩Mt和扭矩Mb。圆盘的最大拉应力σmax和剪应力τmax根据梁单元模型进行计算，可表示为[17]

式中，A为圆盘的截面积，I为圆盘的惯性矩，J为圆盘的极惯性矩，为圆盘的半径。

当σmax或τmax超过其相应的强度极限时，粘结面发生破坏，其中剪切破坏强度σt基于摩尔−库伦准则确定，即

式中，σc表示切向粘结强度，µp是颗粒间内摩擦系数。海冰单元间的粘结强度是离散元计算中的重要参数，其与海冰的温度、盐度、加载速率等因素密切相关。

基于上述理论构造而成的冰脊模型中，脊帆部分在冰脊中占比很小，其离散元参数对冰脊载荷的影响可以忽略不计。而冻结层与平整冰的力学性质较为接近，其离散元参数可直接参考平整冰的校准方法[18]。因此，本文主要针对龙骨部分离散元参数进行研究。龙骨部分的结构十分庞大且位于水下，无法采用试验手段直接测量其力学参数，目前主要通过压剪实验等方法间接研究其力学特性。这里主要参考Heinonen 于2004 年在波的尼亚湾开展的现场压剪试验[10]进行离散元计算(图3)，分析离散元参数对龙骨力学性质的影响并进行计算参数校准。

图3 冰脊压剪试验[10] 及离散元计算模型

1.2 冰脊压剪过程的离散元分析

在 Heinonen 的冰脊龙骨压剪试验中，首先去除龙骨最深处上方附近的脊帆及冻结层部分，然后用一个圆形压板缓慢竖直作用于龙骨上部，直至压板下方的冰块与冰脊彻底脱离，如图3(a) 所示[10]。试验中主要记录压板载荷及运动时程数据，同时在水底放置摄像头记录龙骨底部的变形及破坏现象。图3(b) 为参照其现场试验建立的离散元计算模型，其中龙骨部分由于尺寸远大于压盘，其可以直接近似为长方体，且其边界可设为固定约束。该计算中的主要参数列于表1 中。龙骨的单元粒径取0.2∼0.4 m，由此生成的龙骨计算模型孔隙率为0.31。

图4 冰脊压剪试验中离散元数值及现场测量的载荷与位移关系

表1 冰脊压剪试验中主要计算参数

通过压剪试验的离散元计算得到的压板载荷随位移的变化曲线如图4 所示。对比数值结果与试验数据可以发现，无论在离散元计算还是现场试验中压板载荷都经历了快速线性增长至最大值、较平缓下降、保持在与压板下碎冰堆所受浮力平衡的位置等几个阶段。在ABC的线性加载阶段，可以由力和位移关系的斜率确定冰脊的刚度。龙骨在计算与试验中的整体刚度基本一致，峰值大小也非常接近。

为进一步分析龙骨在压剪试验中的破坏过程，图5 给出了四个不同时刻的冰脊变形状态。左侧为不同时刻的龙骨破坏状态、中间及右侧分别为竖直、水平方向的位移云图。从龙骨的破坏状态可以发现，在u=10.0 mm 的加载初期，龙骨呈局部压缩变形，并主要集中在压板正下方很小范围内；随着压板的向下运动，当u=15.0 mm 时，龙骨内部的压缩变形区域不断增大，受挤压作用的碎冰块不断水平向外扩散引起周围冰块移动形成一个楔形变形区域，楔形部分整体下移并向两侧扩张在边缘处形成带状剪切变形区域；当u=25.0 mm 时，龙骨楔形部分两侧的剪切变形达到极限，龙骨发生整体破坏，形成明显的剪切面，压板载荷开始下降；在u=100.0 mm 的卸载过程中，由于楔形部分受压膨胀与周围龙骨部分发生摩擦，故卸载过程较加载更为平缓。上述过程说明在数值计算中，龙骨呈现出剪切破坏和压缩破坏的混合破坏模式。这一现象与现场试验结果相一致。由于条件限制，现场试验并未取得龙骨破坏过程的内部图像，但通过对比计算结果与室内模型试验结果，可以发现其主要现象基本一致[22-23]。

图5 不同位移下龙骨内部的破坏状态过程及位移云图

1.3 离散元参数对压剪强度的影响

在离散元数值计算中，计算参数对冰脊的压剪失效特性有显著影响。这里选择对冰脊龙骨力学性质影响较大的单元弹性模量Ep、单元间粘结强度σc以及内摩擦系数µp三个参数进行研究。通过不同参数组合下的数值压剪试验，研究以上参数对压剪试验的影响。为分析海冰单元弹性模量Ep的影响，首先在表1 所列计算参数的基础上，取σc=50 kPa，Ep分别取 0.3，0.5，1.0，1.5，2.0 GPa 进行冰脊压剪过程的离散元模拟，由此得到的最大压力如图6 所示。计算结果表明，Ep在 0.3∼1.0 GPa 范围内，Fmax几乎保持不变；而当Ep>1.0 GPa 后，Fmax随着Ep的增大而线性下降。这是由于当Ep<1.0 GPa 时，龙骨在试验中保持以全局性剪切破坏为主的破坏形式，Fmax主要受强度参数σc等影响，与Ep关系不大。Ep主要影响龙骨的整体刚度，与其呈线性正相关关系，如图7 所示；当Ep>1.0 GPa 时，随着Ep的不断增大，龙骨的破坏形式逐渐转变为以连续的局部压缩破坏为主，压板载荷峰值变得密集，最大值Fmax不断减小，如图8 所示。

图6 压剪破坏过程中压板最大载荷与弹性模量Ep 的对应关系

图7 离散元模拟的冰脊整体刚度与弹性模量Ep 的对应关系

图8 不同弹性模量Ep 下压板载荷随竖直位移的变化过程

图9 离散元模拟的最大载荷与粘接强度σc 的对应关系

为分析海冰单元间粘接强度σc的影响，当Ep= 0.5 GPa，σc分别取 60，120，200，300 kPa时，离散元计算得到的压板载荷最大值Fmax如图9所示。可以看出，在所取σc值范围内，Fmax一直随σc的增大而线性增长。但从不同粘结强度下压板载荷达到最大值时龙骨内部裂纹分布情况可以发现 (如图 10 所示)，当σc<120 kPa 时，龙骨保持现场试验中常见的全局性剪切破坏为主的破坏形式；当σc>120 kPa 时，随着σc的增大，压板下方自底部向上方扩展的竖向裂纹逐渐明显；当σc=300 kPa时，龙骨截面只有竖直裂纹，没有明显的剪切带。这说明随着σc的增大，龙骨的破坏形式逐渐由剪切破坏为主转变为弯曲破坏为主。这也与一些室内模型试验中观测到的现象相一致[24]。

图10 不同σc 下离散元模拟的冰脊的破坏模式

为分析海冰单元间摩擦系数µp的影响，取Ep= 0.5 GPa 和σc= 50 kPa，当µp分别为0.2，0.4，0.6，0.8 时，压板载荷最大值时龙骨内部裂纹分布如图11 所示。从中可以发现，这些龙骨在压板作用下均以整体剪切破坏为主，其楔形断裂部分两侧均分布有明显的剪切带，且两侧剪切带与数值方向夹角θ随µp的增大而线性增大。计算结果还进一步表明，压板载荷的最大值Fmax不仅随σc的增大而线性增大，也与µp呈线性正相关关系，如图12 所示。这表明冰脊在压剪过程中表现出明显的摩尔−库伦材料特征。

从以上离散元数值模拟结果可以看出，单元弹性模量Ep、单元间粘结强度σc以及内摩擦系数µp均对龙骨在压剪试验中的力学行为有重要影响。这些参数不仅影响其压剪强度的大小，还决定其破坏形态甚至破坏模式。基于对这些参数影响的分析和现场试验结果，可以选择合理的离散元参数组合对冰脊进行数值计算以得到可靠的结果。

图11 不同µp 下冰脊的失效模式

图12 不同摩擦系数µp 下离散元模拟的最大载荷及失效夹角

2 冰脊与圆形直立结构相互作用的离散元分析

2.1 冰脊与圆形直立结构相互作用的计算模型

冰脊龙骨部分采用表 1 的离散元参数，将冻结层与平整冰简化为统一厚度，其他主要参数列于表2 中。冰脊的冻结层由球体单元按密六方排列形式构造，如图13 所示，其厚度hc为

式中，hc为冻结层厚度，d为单元直径，ne为单元层数。

表2 冰脊与直立结构相互作用中离散元计算的主要参数

图13 由球体单元密六方排列而成的冰脊冻结层

图14 和图15 为结构与冰脊作用过程图，图中颗粒颜色代表速度的大小，而图16 为冰脊对直立结构在水平方向上的冰载荷。从中可以发现，冰脊在整个作用过程中以结构前的局部挤压破坏为主。随着圆形直立结构不断贯入冰脊内部，结构前的碎冰堆积深度不断增大，导致冰载荷也随之增大；当冰脊龙骨发生整体剪切破坏后，冰载荷与结构前碎冰堆积深度开始减小直至龙骨部分冰载荷可以忽略不计，而冻结层则始终保持挤压破坏。由图16 还发现直立结构冰载荷最大值为9.53 MN，冻结层冰载荷最大值为6.50 MN。由此可见，冰脊中的冻结层是对结构作用的主要承载部分，而龙骨的冰载荷相对较小。

图14 直立结构与冰脊相互作用离散元模拟

图15 直立结构与冰脊作用过程的侧视图

图16 离散元模拟的冰脊对直立结构的冰载荷

2.2 冰脊对圆形直立结构冰载荷的验证

为验证计算结果的准确性，这里选择将离散元结果与基于 ISO 19906 规范[18]的计算结果进行对比。ISO 19906 将冰脊载荷Fr分为冻结层载荷Fc与龙骨载荷Fk两部分，即

式中，Fc按照平整冰经验公式计算

式中，CR按照规范推荐取 1.8 MPa，h1取 1 m，hc为冻结层厚度，D为圆柱结构直径，n=−0.5+hc/5，m为−0.16。Fk的计算基于改进的Dolgopolov 公式，即

式中，D为圆柱结构直径；hk为龙骨深度；η为龙骨部分孔隙率；ρw为海水密度；ρice为海冰密度；g为重力加速度；c为宏观上冰脊龙骨表现出的粘聚力，这里取12 kPa；ϕ为龙骨内摩擦角，结合压剪试验结果并基于Meyerhof 理论求得[25-26]

式中，Dind为压盘直径；Ku为被动冰压力系数，一般为常数0.95；Qu为压剪载荷最大值756 kN；R为压剪试验中龙骨完全破坏后的残余压力120 kN。

由此可计算冰脊总载荷Fr、冻结层载荷Fc及龙骨载荷Fk。这里将离散元和ISO 19906 规范结果列于表 3。从中可以发现，离散元计算结果中的Fc相对ISO 19906 规范较小，而Fk则相对较大。这主要是由于ISO 19906 规范对冻结层冰载荷的计算是基于大量平整冰的试验数据并取其较大数值，其相对比较保守以保障结构安全；而在龙骨冰载荷的计算中，离散元方法未具体考虑碎冰块的随机排列和冰块尺寸影响，而是按龙骨中的冰块充分冻结进行计算，从而导致其计算结果略高于ISO 19906 规范。此外，在离散元模拟中，单元间细观计算参数选取的合理性以及冰载荷动力过程的随机性也是影响冰载荷大小的重要因素。但从直立结构冰载荷的总体情况看，离散元计算结果与ISO 19906 规范比较一致，从而说明了离散元方法在计算冰脊对直立结构冰载荷中的适用性。

表3 冰脊对圆柱结构载荷的离散元计算与 ISO 19906 规范对比

3 结论

本文基于粘结球体单元构建了冰脊离散元模型，通过颗粒排列方式与离散元参数的变化表现冰脊冻结层、龙骨等不同部分的力学特征。利用该模型计算分析了一系列冰脊龙骨压剪试验，并与现场实测数据进行了对比。计算结果表明：在合理的离散元参数下该方法能较为准确地模拟冰脊的力学行为；单元弹性模量Ep、单元间的粘结强度σc以及内摩擦系数µp不仅影响压板载荷的大小，还决定龙骨受压后的破坏形式。在一定参数范围内龙骨主要以全局性剪切破坏为主，压板载荷会随σc和µp的增大而线性增长，且剪切面与垂直方向的夹角也会随µp的增大而增大。最后计算了冰脊与圆柱形直立结构的相互作用过程，其计算结果与ISO 19906 规范相近，由此证明了该方法在工程计算中的可靠性。离散元方法能模拟冰脊与直立结构作用中的堆积效果与破坏模式，在更复杂的工况计算中将具有良好的应用前景。