基于黏聚单元法计算分析海洋平台锥形立柱冰载荷

詹开宇，曹留帅，万德成

(上海交通大学 船海计算水动力学研究中心(CMHL) 船舶海洋与建筑工程学院，上海 200240)

在寻求新资源和新能源的大环境下，越来越多的国家将目光聚集到了极地。南北极地区虽然常年被冰雪覆盖，但其辽阔冰原下覆盖的丰富资源是人类的一笔宝贵财富。近年来，人类对于海洋资源的开发逐步向寒区延伸和发展。在极地探索研究中，需要考虑冰载荷对海洋结构物的影响。冰载荷是寒区的主要环境载荷。极地的海冰以多种形式存在，如层冰、碎冰、冰脊和冰山等，其中层冰是主要存在形式之一，也是研究海冰与结构物相互作用时最常考虑的载荷形式。

关于结构物层冰载荷的测量与预报，国内外学者最早使用的方法是直接测量和试验模拟方法，并总结出了一系列用以估算层冰载荷的经验公式[1-3]。但是直接测量的难度较大，且设备的安装和维护成本较高；试验方法对试验环境的要求较高，且对于不同类型的海洋结构物和不同形式的海冰普适性不佳。近年来，学者们开始采用多种数值方法进行冰载荷模拟尝试，其中得到广泛应用的是离散单元法(discrete element method，简称DEM)和有限元法(finite element method，简称FEM)。离散单元法将海冰模拟成黏结在一起的颗粒，通过改变颗粒间作用形式和作用力的参数，能够较准确模拟处平整冰的破坏和堆积效果，在模拟海冰离散特性方面有较明显的优越性；有限元法的重点和优点在于可以较好地研究海冰结构的变形。季顺迎等[4]利用离散元法对海冰与直立结构相互作用进行了数值模拟，获得了不同柱径下的冰载荷和结构冰振响应。Kim等[5]采用新开发的有限元模型和模型试验方法，对一艘货船在浮冰条件下的破冰阻力进行了数值和试验研究，对比结果得出了一系列重要结论。

黏聚单元法是在有限元方法的基础上逐渐发展起来的，该方法在有限单元间插入了黏聚单元，因此在研究中可以同时考虑材料的变形和破坏。黏聚单元法可以较准确模拟出海冰的变形、破碎和堆积现象，如今已经在海冰材料模拟研究中得到了广泛应用。刘路平[6]使用黏聚单元法进行了平整冰和抗冰海洋平台相互作用的数值模拟研究，根据以往的试验结果对模拟结果进行了验证。王峰等[7]基于非线性有限元数值方法，引入黏聚单元模型并结合线性软化弹塑性本构模型，对平整冰与竖直固定圆柱体碰撞进行了数值模拟。

在海冰覆盖率较高的寒区，通常在海洋平台立柱的水线面附近设置抗冰结构，以增加平台的抗冰性能，常见的抗冰结构为抗冰锥体。由于海冰的抗剪强度远小于其抗压强度，因此在海冰与海洋平台相互作用时，带有特定锥角的抗冰锥体可以将海冰的主要破坏形式从挤压破坏转变为弯曲和剪切破坏，从而极大的减少平台所受载荷作用，过程如图1所示。抗冰锥体是一种非常有效的抗冰方式，对于抗冰锥体的研究是十分必要的，许多学者都对冰区的锥体结构进行了探索和研究[8-10]。

图1 层冰与锥形结构相互作用后的破坏与堆积Fig. 1 Failure and accumulation of ice after interaction with conical structure

基于黏聚单元法建立海冰模型，并进行层冰和带有抗冰锥角的固定平台立柱相互作用的数值模拟。该数值模拟将探究海冰的厚度等条件对结构物受力的影响。同时，对于带有抗冰锥体的立柱，锥体的角度对立柱的受力也存在一定的影响，也会改变海冰的破碎和堆积现象，这也是本文研究的另一重点。将数值模拟的计算结果与近似环境条件下的试验和观测结果、使用离散单元法模拟得到的结果以及基于ISO(international organization of standardization)冰力标准计算得到的结果进行对比验证，在一定范围内具有较好的一致性。研究表明，基于黏聚单元法建立的海冰模型在研究层冰与带有抗冰锥角的固定平台立柱相互作用时具有较好的适用性，对于后续有关寒区结构物的冰载荷模拟以及抗冰结构的设计具有重要的参考意义。

1 基于黏聚单元法的层冰模型建立

海冰是一种性质复杂的温度敏感性复合材料，其物理和力学性质与其本身的盐度、孔隙率，环境温度加载速率等都有密切关系，在不同条件下会表现出脆性和塑性两种特性。此外，在海冰与不同形式的结构物相互作用时，其破坏形式也是多种多样的，包括弯曲破坏，挤压破坏，剪切破坏等。因此在进行冰载荷的研究和数值模拟时，确定海冰的本构模型、材料参数和破坏形式，对得到准确的模拟结果十分关键。

有限元法是一种研究冰载荷的有效方法，可以模拟出海冰的变形和破坏，也能快速得到结构物的受力和响应情况，目前得到较广的应用。但该方法存在局限性，即经过碰撞破坏后的冰单元会失效并被删除，因此有限元法并不能很好地模拟出层冰破碎后的堆积现象以及碎冰对于结构物的后续作用，同时造成整个模拟过程的质量不守恒，使模拟结果与实际情况有偏差。为了解决这一问题，学者们提出了很多的解决方案，如FEM和SPH(Smoothed Particle Hydrodynamics)方法耦合，黏聚单元法等[11-13]。

黏聚单元法是在有限元法基础上优化和发展得到的一种方法。如图2所示，该方法将使用有限元方法建立的海冰单元离散化，并且在冰单元间插入黏聚单元。插入的黏聚单元与冰单元共用节点，且厚度为零。图3是构造的海冰的8节点六面体单元，水平和垂直方向的黏聚单元将海冰单元分隔。

图2 黏聚单元与冰实体单元Fig. 2 Cohesive element and bulk element

图3 海冰六面体单元模型与黏聚单元模型Fig. 3 Ice hexahedral element model and cohesive element model

在层冰和结构物相互作用的过程中，原先的海冰有限单元发生变形，通过与黏聚单元共用的节点传递力和位移；而插入的黏聚单元会根据海冰材料的本构模型和材料属性参数发生变形，当其位移达到设定的破坏临界点时，黏聚单元破坏，而与之相连的海冰单元将会脱落。在整个层冰的黏聚单元模型中，破坏仅在黏聚单元处发生。

黏聚单元的破坏是以给定的牵引力—位移曲线作为依据进行判定的。该曲线规定了黏聚单元所受的牵引力与相邻单元的分离位移之间的关系。牵引力—分离位移准则的重要参数有断裂应力、断裂韧性、最大分离位移、曲线形状等。根据上述参数的不同可以将该模型分为不同种类，目前常用的黏聚单元破坏准则模型有三种：线性软化模型；指数软化模型以及延性软化模型(又称梯形软化模型)[11]。图4为三种模型的牵引力—分离位移曲线。其中，曲线与坐标轴所围成的区域的面积即使黏聚单元破裂所需的能量。下文使用的是黏聚单元线性软化模型。

图4 黏聚单元破坏准则模型的牵引力—分离位移曲线Fig. 4 Common shapes of traction-separation law

2 数值模拟参数设置

2.1 结构物模型参数

结构物模型为带有抗冰锥角的海洋平台立柱，来源于渤海湾JZ20-2MUQ平台[12]。取一根立柱与海冰发生作用，未加锥角之前立柱的直径为1.7 m，锥体高度2.5 m，锥体最大直径为4 m。定义锥体斜面与水线面夹角θ为锥角大小，并保证水线面处(即立柱与层冰作用位置)锥体直径一致，为2.8 m。平台和锥形结构物如图5所示。

图5 渤海湾JZ20-2 MUQ平台与锥形结构物示意Fig. 5 JZ20-2 MUQ platform on Bohai Bay and the model of conical structure

碰撞过程为层冰以恒定速度撞击固定立柱。碰撞不考虑平台变形，因此选择立柱的材料为刚性材料。具体的材料属性如表1所示。

表1 锥形立柱材料属性

2.2 海冰模型参数

海冰的有限元单元实体模型采用8节点六面体单元建模，取18 m×40 m的层冰进行计算。研究了层冰厚度对于冰力大小的影响，冰厚h分别取0.2 m、0.25 m、0.3 m、0.35 m和0.4 m。采用各向同性的线性弹塑性材料模拟海冰实体单元的变形。在冰实体单元间插入厚度为零的黏聚单元项，黏聚单元的牵引力—分离位移曲线选择线性软化曲线模型。

海冰材料的基本参数受到许多因素的影响，在不同工况下很难有确定的数值。国内外学者对海冰进行了大量的力学试验，并考虑海冰的温度、盐度、孔隙率等对海冰性质的影响，给出了工程应用中所使用的海冰参数的推荐范围[13-14]。这里也根据该推荐范围选择了海冰弹塑性模型的材料参数，以满足本算例的工程需求。具体参数值见表2。

表2 海冰单元与黏聚单元模型参数

2.3 网格收敛性分析

网格尺度和时间步长对计算时间和数值结果有重要影响。在显示动力学分析中，时间步长通过最小单元尺寸确定，因此本研究中仅对网格收敛性进行分析。分别选取尺寸为0.075 m×0.075 m×0.01 m、0.05 m×0.05 m×0.01 m和0.025 m×0.025 m×0.01 m的网格，计算层冰以0.4 m/s的速度与锥形立柱的碰撞过程。模拟时长20 s并统计水平方向平均冰力大小、相对误差和计算所需时间，结果如表3所示。可以看到，随着网格尺寸的增加，计算结果逐渐收敛，同时计算时间迅速增加。权衡计算精度和计算时长，本研究选择尺寸为0.05 m×0.05 m×0.01 m的网格进行计算。

表3 网格尺寸对计算结果和计算时间的影响

2.4 接触模式

计算前设定层冰和结构物碰撞的接触模式。立柱和层冰之间的接触采用面面接触的侵蚀算法，在LS-DYNA中选择算法“CONTACT_ ERODING_ NODES_ TO_ SURFACE”，该算法可以定义海冰单元的破坏失效标准。层冰与立柱之间的摩擦力系数设定为0.2，破碎堆积的碎冰之间的接触采用单面接触的算法，在LS-DYNA中选择接触算法“CONTACT_ ERODING_ SINGLE_ SURFACE”，该算法可以简化碎冰之间的作用，节约计算时间。海冰间的摩擦系数设定为0.1。

2.5 计算域与边界条件设置

如图6所示的计算域中，海冰沿x轴负方向移动，z轴与海冰厚度方向平行，与水线面垂直，且定义向上为正方向。层冰与立柱碰撞的边界设定为自由边界，其余三边约束三个方向自由度，并设置为无反射边界条件，以模拟半无限大的层冰。本算例没有考虑水对碰撞过程的影响，仅在海冰下方增加海平面，并在垂直方向定义海水对海冰的浮力作用。

图6 计算域设置Fig. 6 Computational domain settings

3 计算结果与分析

3.1 层冰破坏与碎冰堆积

将使用黏聚单元法得到的模拟结果与国内外学者的实测、试验以及数值模拟等不同方法得到的结果进行对比，从而验证该方法的可行性。

图7给出了层冰以0.4 m/s的速度撞击固定锥形立柱的模拟结果，图中所示为模拟进行到t=10 s时的俯视图和侧视图。观察碰撞过程可知层冰与抗冰锥角斜面接触和碰撞情况。海冰与斜面接触时会在斜面上有一小段的爬升，同时海冰所受的主要载荷为弯曲载荷。当海冰的弯曲应力大于其抗弯强度时，海冰发生破坏，形成碎冰并在水线面附近堆积。随着层冰的不断推进，碎冰和层冰继续与立柱相互作用。模拟中还观察到碎冰在锥体斜面上爬升后滑落的现象。分析完整作用过程可知，层冰与固定锥形立柱相互作用时产生的冰载荷主要有两个部分，一部分为层冰对立柱的作用，另一部分为破碎后堆积的碎冰对立柱的继续作用。对比可知，本文模拟得到的一系列现象与渤海湾JZ20-2 MUQ平台的实测结果以及HSVA(Hamburg Ship Model Basin，汉堡船模试验水池)的模型试验结果基本吻合[15]，如图8所示。

图7 t=10 s时刻俯视图和侧视图Fig. 7 Top and side views at time t=10 s

图8 模拟结果与实测和试验结果对比Fig. 8 Simulation results compared with measured and experimental results

3.2 立柱所受冰载荷

图9为层冰运动速度为0.4 m/s时锥形立柱所受x方向水平载荷的时历曲线。从图中可知，在与层冰碰撞的过程中，立柱所受冰载荷呈波动变化且幅值较大。这是因为层冰与立柱碰撞过程中，层冰一直在进行“接触变形—破碎—堆积爬升—接触变形”这一循环过程，海冰受到压应力和弯曲应力作用，载荷上升，随后黏聚单元破坏，应力卸载造成载荷下降。同时由图9可知，碰撞前期冰载荷呈震荡上升趋势，随后逐渐趋于稳定，在一定范围内波动。这是因为在碰撞初期，随着层冰不断破坏，碎冰开始堆积，海冰与立柱的接触面积逐渐增大，冰载荷开始上升。当碎冰堆积达到一定程度时，接触面积不再有明显变化，层冰碰撞与碎冰堆积的过程达到稳定，冰载荷开始在一定范围内稳定波动。这一结果也与模型试验得到的冰载荷结果相吻合。

3.3 抗冰锥角大小对冰载荷的影响

3.3.1 ISO冰力标准

ISO-19906标准[16]是广泛应用于海洋结构物设计的行业标准，其中给出了不同类型的海冰与海洋结构物相互作用的冰载荷。使用冰力标准给出的层冰与锥形结构物相互作用过程中冰载荷计算公式计算得到了不同锥角条件下立柱所受最大冰载荷的大小，并将其与本文研究结果进行对比，以验证该黏聚单元模型在冰载荷研究中的可行性。ISO标准中层冰与锥形结构物相互作用的冰载荷计算公式主要基于海冰的塑性破坏和脆性破坏。海冰的塑性破坏通常发生在载荷加载速率较小的情况下，海冰在外力作用下的变形量超过其最大弹性变形量，产生明显的永久变形，进而破坏失效；当加载速率增加时，海冰破坏模式逐渐向脆性破坏转变；在大速率外载荷冲击下，海冰会发生无明显变形的脆性破坏。选取海冰的塑性破坏模式进行计算和对比。

塑形破坏模式中，锥形结构物所受冰力由两部分组成：层冰弯曲断裂和碎冰爬升对结构物的作用力。层冰破坏时的水平冰力HB和垂直冰力VB表示为：

(1)

VB=HB·hv

(2)

碎冰在锥面爬升时的水平冰力HR与垂直冰力VR表示为：

(3)

(4)

其中，f，gr和hv等系数由冰—结构物间摩擦系数、锥径、锥角大小和碎冰爬升高度等参数定义。

则在锥形结构上的总冰力可以表示为：

FH=HB+HR

(5)

FV=VB+VR

(6)

3.3.2 锥角影响

为了研究抗冰锥角θ的大小对冰载荷的影响，取锥角θ分别为40°，45°，50°，55°，60°和65°进行了一系列的模拟研究。为消除锥角改变时造成的锥径变化带来的影响，保持水线面处锥径不变，为2.8 m。计算了层冰以0.4 m/s的速度碰撞锥形立柱90 s时间内立柱所受水平方向和竖直方向的最大冰力与平均冰力大小，并将结果与ISO冰力标准计算得到的最大冰力以及狄少丞等[15]使用离散单元法得到的最大冰力结果进行对比。如图10和图11所示。

图10 水平方向冰载荷与θ大小的关系Fig 10. The relation between θ and ice force of horizontal direction

图11 竖直方向冰载荷与θ大小的关系Fig 11. The relation between θ and ice force of vertical direction

由计算结果可知，水平方向的最大冰力和平均冰力随着抗冰锥角的增大而增大；竖直方向的最大冰力在一定锥角范围内变化不大，而平均冰力随着锥角的增大而减小。与其他结果对比可知，使用黏聚单元法得到的结果与其他方法得到的结果吻合较好，在变化趋势上保持一致，但整体的冰载荷结果偏小。相较ISO冰力标准，造成这一结果的原因可能是黏聚单元和冰体单元参数的选择和实际海冰材料不完全一致；另外，实际海冰材料内部结构复杂，破碎形状随机，而使用黏聚单元法建立的海冰模型是规则的各向同性结构，这也是造成误差的主要原因。相较DEM，黏聚单元法在模拟海冰结构的变形和破坏时具有一定优势，由于黏聚单元是插入在海冰单元之间的零厚度片状单元，可以承受拉应力和剪应力但几乎无法承受压应力，而DEM中颗粒间的连接基于弹簧阻尼系统，可以更好地模拟海冰多个方向的受力，因此模拟结果更接近实际情况。这里使用的黏聚单元模型需要在这一方面进一步改进和优化。

考虑层冰与锥形结构物作用过程，锥角对冰载荷造成影响的主要原因可能是海冰与锥形结构物作用时，海冰的破坏形式发生了变化。海冰与锥体作用时的破坏形式主要由弯曲破坏和挤压破坏，两种破坏形式同时存在。当锥角较小时，海冰的破坏以弯曲破坏为主，随着锥角增大，弯曲和挤压同时发生且挤压破坏逐渐称为主要破坏形式，而海冰的抗压强度要明显大于抗弯强度。因此，随着锥角增大，x方向上的冰力逐渐增大；同时，随着锥角增大，碎冰在斜面上的爬升越来越困难，爬升高度降低，因此z方向的冰力呈减小趋势。

3.4 冰厚对冰载荷的影响

还研究了冰厚对冰载荷大小的影响。冰厚分别取0.2 m，0.25 m，0.3 m，0.35 m和0.4 m进行模拟。冰速取0.3 m/s，其余参数不变。图12给出了t=60 s时间内x方向最大冰力和平均冰力与冰厚的关系。可以明显看出，冰力大小与冰厚呈正相关，且增加速度接近线性。将计算结果与使用ISO冰力标准计算得到的结果相对比，具有较好的一致性。

图12 冰厚与冰载荷的关系Fig. 12 The relation between ice thickness and ice force

4 结 语

基于黏聚单元法建立了海冰模型，并模拟了带有锥形结构的立柱与层冰相互作用的动力过程，分析不同参数对碰撞过程和冰载荷的影响，得到的主要结论如下：

1) 黏聚单元法可以较准确地模拟层冰与锥形立柱的碰撞过程，该方法在模拟海冰的破碎、堆积以及碎冰在立柱上的爬升等方面有较好的表现，模拟结果与实测结果以及试验结果基本吻合。

2) 抗冰锥角的大小对冰载荷有明显影响。研究结果表明，锥角角度在40°左右时有较好的抗冰效果。随着锥角增大，水平方向的冰载荷有明显增加，而竖直方向冰载荷逐渐减小。造成这一结果的主要原因为锥角改变了海冰的破坏形式。该模拟结果与ISO冰力标准以及其他学者得到的结果吻合。

3) 研究了层冰冰厚对冰载荷影响。结果表明，锥形立柱所受层冰载荷与冰厚呈正相关，且增长速度接近线性。该结果与ISO冰力标准吻合较好。

在下一步工作中，将增加计算工况，进一步研究影响冰载荷的因素。此外，将对海冰的黏聚单元模型进行改进，同时优化网格，提高计算结果的准确性和可信度。