基于S-ALE 方法的圆柱体垂直出水破冰研究1)

汪春辉 王嘉安 王 超 郭春雨 朱广元

(哈尔滨工程大学船舶工程学院，哈尔滨 150001)

引言

在极区，平整冰与海洋结构物间的相互作用形式主要为挤压破坏和弯曲破坏，冰层受结构物作用产生断裂破碎、裂纹扩展和滑移清除过程[1-4].在某些特殊情况，结构物需要上浮穿透冰层破出水面，对冰层施以自下而上的垂向载荷，而冰层垂直贯穿破裂过程十分复杂.

学界对冰与结构物垂直相互作用的研究还不充分，国内外相关研究主要以试验和理论方法为主.由于试验条件要求苛刻及数值技术的发展，一些具备完善冰作用机理的数值方法成为解决这些问题的好方法.并且海冰是一种复杂的材料，有许多因素会影响冰的物理特性，例如，温度孔隙率和盐度、速率依赖性.由于冰可发生大而复杂的变形，模拟冰与结构的相互作用是一个巨大的挑战.

目前，有限元方法(FEM)、离散元方法(DEM)、光滑粒子法(SPH)和近场动力学(PD)等数值模拟方法已被用于模拟冰与结构物间的相互作用，本文采用有限元方法开展圆柱体垂直出水破冰研究.针对水下结构物在冰层下方上浮至破出水面的过程的研究在世界范围内仍处于起步阶段，国内外公开资料极其稀少.因此利用有限元方法分析结构物与冰层垂直相互作用进展可以参考船与冰相互接触的研究资料[5-7].

针对冰与结构物耦合作用有限元数值分析方法，国际上已有大量的研究探索.Ehlers 和Kujala[8]在LS-DYNA 有限元软件中使用弹塑性应变率材料模拟冰的弯曲行为，并与4 点弯曲试验结果进行对比.Polach 和Franzvon[9]应用LS-DYNA 有限元软件采用MAT_DAMAGE_3 材料模拟冰的力学特性，并基于损伤力学对数学模型进行分析，将模拟结果与冰池试验结果对比.Gagnon[10]提出了一种弹性线性硬化材料，可模拟冰力学特征，但不能模拟冰的破坏断裂行为.Wang 和Derradji-Aouat[11]基于罚函数法及拉格朗日欧拉法模拟碎冰对系泊结构的作用过程，以评估系泊结构的整体载荷.Ranta 等[12]通过FEM-DEM 耦合方法构建冰与结构物相互作用二维数值模型，利用多元回归线性数学模型和变量消元法对冰载荷特征值分析，探究冰载荷影响因素.国内方面，郭春雨等[13]基于LS-DYNA 有限元软件模拟了极地船和浮冰的耦合作用过程，实现了冰区船舶在碎冰区域的航行阻力预报.Wang 等[14]应用FEM方法与DEM-CFD 方法模拟了冰区船舶在冰缘区的航行过程，并将数值模拟结果与模型试验对比，分析了碎冰载荷的主要影响因素.黄焱等[15]基于FEM方法构建并分析了固定约束及摩擦约束边界条件下的冰温度膨胀力，发现固定约束会造成高度的应力集中.任奕舟和邹早建[16]、王峰[17]基于FEM 方法，通过粘聚力单元和弹塑性冰本构模型构建船舶与海冰碰撞场景的数值仿真，验证了粘聚元方法模拟海冰的准确性与可行性.Chen 等[18]应用FEM-SPH 方法建立了锥体和破冰船与平整冰相互作用的数值模型，预测了破冰阻力及碎冰堆积过程.

本文应用 LS-DYNA 有限元软件采用结构化-任意拉格朗日欧拉 (S-ALE)流固耦合方法及罚函数接触算法，建立了冰、水、结构物耦合作用数值模拟方法，自主搭建了圆柱体垂直贯穿冰层试验台架.通过冰梁结构4 点弯曲及楔形体入水数值模拟，验证了弹塑性应变率本构材料表征冰力学特性的准确性和流固耦合方法的可行性，对比了无水环境下圆柱体垂直贯穿冰层和圆柱体垂直出水破冰过程中的冰载荷、冰层挠度变化、冰层应力响应及断裂现象，以期为冰层与海洋结构物垂向相互作用研究提供参考.

1 理论基础

LS-DYNA 有限元软件被广泛应用于航天、航海和武器装备等工程制造领域.由于LS-DYNA 有限元软件包在结构计算、材料损伤破坏等领域的优势，越来越多的学者将其应用于船冰碰撞、冰桨切削等冰工程方面的研究中，运用Lagrange 格式的显式中心差分法计算大变形的结构物与冰层耦合作用问题.本文采取结构化任意拉格朗日欧拉方法，基于罚函数流固耦合算法模拟圆柱体上浮垂直出水破冰中冰、水、结构物耦合作用的强非线性动态过程.

1.1 结构化-任意拉格朗日欧拉流固耦合方法

结构化任意拉格朗日欧拉方法原理上是一种基于统一连续介质的动量守恒方程的二阶精度的平流算法[19].通过接触算法获取固体之间的接触作用时的界面力，能够有效解决物体在空间上存在大位移同时本身存在大变形问题.在S-ALE 计算中，定义了材料域和参考域，材料的速度为u，参考系的移动速度为v，速度差u-v表示为w，由此可对S-ALE 计算方法的控制方程进行推导，如下所示

式中J为雅可比矩阵，表示两个域的相对体积微分.

S-ALE 公式中的材料时间导数写成公式

式中，br表示b的参考域函数.

S-ALE 方法的质量、动量以及能量守恒方程如下

1.2 罚函数接触算法

罚函数接触算法将可能发生接触作用的两个表面分别定义为主表面和从表面，发生穿透时，类似在主表面和从表面之间放置一系列法向弹簧，弹簧的作用是为了限制穿透，弹簧的力被称为界面接触力.

罚函数法的接触力由下面计算公式得出

式中k为单元接触面刚度，δ为穿透量.

单元接触面刚度k由如下公式决定

式中pf为单元接触面刚度的惩罚因子，K是接触单元体积弹性模量，A是接触段面积，V是主段体积.

在数学方法上，用欧拉变量描述流体的动态，用拉格朗日变量描述结构的运动边界，通过分布节点力和插值速度来表示流场和结构物的交互作用.该方法适用于广泛的流固耦合应用，界面接触力为估算求得且可能存在数值振荡[20].

1.3 冰本构模型

冰的特性取决于温度、应变率和载荷方向，冰在受力变形破坏过程中，受加载速率的影响，会呈现韧脆转换现象[21].本文选用冰本构材料为弹塑性应变率本构材料*MAT _PLASTICITY_COMPRESSION_TENSION，该材料为一种各向同性的弹塑性材料，能够模拟不同的压缩和拉伸响应.此外，它还考虑了应变率敏感性[22-24].

该材料模型需要定义两条用于拉伸与压缩行为的屈服应力对有效塑性应变的曲线，以及一条应变率对压缩屈服应力比例系数(CYSF)的曲线，用于表征应变率对压缩强度的影响.假设弯曲强度与应变率无关，拉伸和压缩下的屈服应力(σy)与有效塑性应变(εp)的曲线均由以下方程推导得出

式中σ0为初始屈服应力.

假定压缩和拉伸过程中冰材料的塑性硬化模量(Ep)相同.通过Sánchez 等[25]给出的冰材料压缩强度应变率敏感性拟合函数可以计算得出应变率()为10-3s-1时冰材料的初始压缩屈服强度()

为了定义应变率对压缩屈服应力比例系数的曲线，不同应变率下的压缩强度由应变率为10-3s-1时的初始压缩强度进一步归一化得出，本文采用的海水冰材料应变率为10-3s-1时的初始压缩强度为5.8 MPa.海水冰应变率和相应的压缩屈服应力比例系数如表1 所示.对于拉伸行为，实验结果显示脆性破坏及应变率对冰拉伸强度的影响可以忽略不计[26].

表1 海水冰应变率和压缩屈服应力比例系数Table 1 Strain rates and compressive yield stress scale factors of sea ice

该材料模型采用等效塑性应变破坏准则，当单元的等效塑性应变达到定义的临界值(εf)时，该单元将被删除.可以通过调整εf值以获得适当的弯曲响应.

2 数值计算模型

2.1 无水环境下圆柱体垂直贯穿冰层有限元模型

无水环境下圆柱体垂直贯穿冰层数值模型包含以下两个部分:冰层和圆柱结构体，均采用拉格朗日算法进行描述.模型的几何尺寸见下图1，模型冰厚为3 cm，长宽均为50 cm;圆柱体直径为3 cm，高为5 cm.冰层边界条件采用四边约束，圆柱体置于冰层正中心下1 cm 处.

图1 数值模型尺寸Fig.1 Size of numerical model

本文将圆柱体定义为刚体，冰层本构选用弹塑性应变率材料模型，根据数值模拟需求将冰材料分为淡水冰及海水冰两类.主要材料参数如表2所示[24-26].

表2 主要材料参数Table 2 Main material parameters

圆柱体与冰层的相互作用以及冰单元间相互作用通过侵蚀接触算法实现，冰层单元不断发生破坏并删除失效单元形成裂纹，一旦网格单元失效且接触表面发生变化，通过侵蚀接触可重新建立表面形状[21].

2.2 圆柱体垂直出水破冰有限元模型

圆柱体垂直出水破冰数值模型设置与文2.1 节开展的圆柱体垂直贯穿冰层试验工况保持一致.包含4 部分:空气域、水域、圆柱结构体和冰层，其中，空气域和水域采用欧拉算法进行描述，圆柱结构体和冰层仍采用拉格朗日算法进行描述.空气域和水域的几何尺寸及网格划分如图2，空气域和水域的采用渐进式网格划分，即应用*ALE_STRUCTURED_MESH 和*ALE_STRUCTURED_MESH_CONTROL_OINTS 关键字生成S-ALE 结构化网格.水和空气通过*MAT_NULL 材料来描述，水的状态方程为OS-GRUNEISEN 来实现，空气的控制方程为EOSINER-POLYNOMIAL，空气与水的相关模型参数设置如表3，圆柱体结构及冰层结构参数设置与无水环境的数值模型参数设置保持一致.

表3 空气、水相关参数设置Table 3 Setting of air and water related parameters

图2 圆柱体垂直出水破冰数值模型Fig.2 Numerical model of vertical upward water of cylinder breaking through level ice

针对圆柱体垂直出水破冰流固耦合问题，通过INITIAL_HYDROSTATIC_ALE 关键字实现流体区域静水压力分层，使用*CONSTRAINED_AGRANGE_IN_SOLID 关键字实现冰层、圆柱体和流体区域间的流固耦合作用.

3 计算方法验证

3.1 冰本构模型验证

Sodhi[27]指出冰层在垂向载荷作用下的破坏方式主要以弯曲破坏为主.通过对冰梁的4 点弯曲模型试验进行数值模拟，验证本文所用应变率弹塑性材料模拟冰弯曲破坏过程的可靠性.以Kujala 等[28]冰梁弯曲破坏模型试验为依据，模型几何尺寸如图3，上方的支撑架刚性固定，下方的支撑架以2.756 mm/s匀速向上移动.冰材料参数与表2 中的海水冰材料一致.

图3 冰梁四点弯曲数值模拟几何模型Fig.3 Geometric model for numerical simulation of four-point bending of ice beam

冰梁受到弯曲载荷后在加载支撑处形成断裂裂纹如图4(a)，通过对比文献中基于SPH 方法的断裂现象即图5(b)发现，本文模拟的冰梁4 点弯曲断裂现象与试验[28]及SPH 方法[23]计算结果吻合度极高.

冰梁4 点弯曲过程的弯曲载荷-时间曲线如图5，试验及数值计算结果的曲线趋势以及峰值点一致性良好.结合图4 冰梁4 点弯曲破坏断裂现象可以说明本文所用应变率弹塑性材料可以准确模拟冰层弯曲破裂问题.

图5 冰梁4 点弯曲数值验证模拟结果Fig.5 Numerical verification of simulation results of four-point bending of ice beam

3.2 圆柱体垂直贯穿冰层试验

基于自主开发的试验台架、以冻结模型冰为试验对象、开展无水环境下圆柱体垂直贯穿冰层试验，探究该动态破冰过程的现象及破冰载荷.试验台架如图6，主要仪器包括伺服电动缸、数据采集与分析处理系统和接触式压力传感器等，试验台架主要用来固定模型冰以及搭载伺服电动缸

图6 试验用主要仪器Fig.6 Main instruments for test

试验使用冻结型模型冰作为模型材料，将煮沸后的纯净水倒入模具中冷却凝固至定型，冻结24 h定型后取出进行尺寸修正，再将修正过的模型冰放于低温试验箱中储存24 h 后进行试验，低温试验箱环境温度恒定为-25 °C.为测量制备模型冰的力学性能参数，分别进行了单轴压缩试验和巴西圆盘试验测量其抗压强度及拉伸强度.单轴压缩试验和巴西圆盘试验压头加载速率均为4 mm/s，与圆柱体垂直贯穿冰层试验加载速率保持一致.

单轴压缩试验试验选用长方体模型冰试样，冰试样尺寸为40 mm×40 mm×100 mm，尽可能减少冰体试样本身对试验结果的影响，测得模型冰试样平均抗压强度为2.41 MPa.单轴压缩试验典型破坏过程如图7.

巴西圆盘法是一种间接测量抗拉强度的方法，能够有效避开直拉法试验过程中存在的冰试件夹持和加载困难的缺点[29].试验用圆盘状冰体试样直径为10 cm 厚度为5 cm，测得模型冰试样平均拉伸强度为0.39 MPa.巴西圆盘试验典型破坏过程如图8.

Sodhi[27]发现冰层的垂向贯穿破裂过程中的大挠度变形仅发生于施加载荷附近区域，通常小于冰厚度的10～ 20 倍，且固定约束(四边约束)下冰载荷稳定性优于自由约束下冰载荷稳定性，但固定约束(四边约束)下冰载荷均值要大于自由约束下冰载荷均值.为获得可靠的试验结果及规律性结论，选取符合尺度效应的冰层尺寸(50 cm×50 cm×3 cm)的模型冰板并采用固定约束(四边约束)进行多组试验，如图9.将模型冰板置于试验台架上，铝型材扣置于模型冰上方约束冰层运动.试验用圆柱体试样直径3 cm，高5 cm，压载速率为4 mm/s，圆柱体试样匀速上升与冰层接触，通过数据采集系统和高速摄像机，记录模型冰板断裂破坏过程中冰载荷的时历曲线以及圆柱体试样与模型冰接触作用时冰层破裂现象.

图9 试验用冻结模型冰Fig.9 Frozen model ice used in the test

圆柱体试样与冰层接触后，冰层试样不会瞬间破裂，如图10(a)，随着圆柱体压头的上压，冰层中间部分受压发生向上弯曲，冰层四周向上翘曲，如图10(b).随着圆柱体试样的进一步上压，裂纹从冰层中部产生，瞬间扩散至冰层四周，并在四边约束的作用下，冰层先会沿约束内侧产生周向断裂，最终被顶起形成碎冰堆，如图10(c).

图10 试验现象Fig.10 Experimental phenomenon

图10 试验现象(续)Fig.10 Experimental phenomenon (continued)

由于在冰层的制取过程中，冰的形状、厚度存在一定不可避免的差异，无法保证每次试验模型完全相同，导致冰载荷极值与作用时间存在一定的差异.选取4 组3 cm 厚度冰层进行试验重复性验证，其试验结果如图11，结果表明整体试验数据一致性良好.

图11 冰层垂直贯穿破裂试验重复性验证Fig.11 Repeatability verification of vertical penetr-ation test of level ice

3.3 圆柱体垂直贯穿冰层数值模型验证

网格质量影响数值计算精度，过多的网格往往会极大的延长数值计算时间，而数值计算采用的模型计算单元需要保证结果准确性的同时提高计算效率.由于圆柱体模型为刚体材料，网格质量影响相对较小，所以为验证网格合理性，分别对冰层数值模型划分大小为2.5 mm，3.5 mm 以及4 mm 的正方形网格，进行网格无关性分析.

如图12 所示，图中3 条曲线总体来看趋势一致且峰值点近似相同，呈振荡上升趋势且冰载荷特征稳定，2.5 mm 网格单元曲线与3.5 mm 网格单元曲线拟合度良好，4 mm 网格单元曲线载荷卸载时间相对提前.所以选取冰层网格单元大小为3.5 mm，既可以保证数值精度又可以提高计算效率.

图12 不同网格尺寸下的冰载荷时历曲线Fig.12 Time history curve of ice load under different element sizes

开展无水环境下圆柱体垂直破冰数值仿真以验证数值方法准确性，数值计算模型见1.2 节，冰层材料选用表2 中淡水冰材料，计算工况与上述试验一致，图13 为数值模拟结果.从图13(a)应力云图中可以看出，当圆柱体与冰层接触而冰层并未发生断裂现象时，冰层发生弹性变形，应力集中在冰层中心，由于冰层采用四边约束应力区域呈环状.随着圆柱体与冰层的持续作用，部分冰层单元达到失效极限，出现径向裂纹向四周扩展，裂纹相互作用进一步形成楔子，如图13(b)所示.当圆柱体继续上升，由于应力集中于冰层中央，导致冰层中部破碎形成碎冰堆积，冰层产生周向裂纹，随后冰层向中心塌陷，如图13(c)所示，对比图10(c)可以看出冰层向中心塌陷整体现象以及破坏范围均呈现良好的一致性，且均存在约束下的径向裂缝在其形成后不再进一步扩展，并且由径向裂缝形成的大块楔形物随圆柱体进一步加载而随之推动.

图13 圆柱体垂直破冰数值模拟冰层破裂过程Fig.13 Numerical simulation of level ice rupture process by cylinder vertical ice breaking

图13 圆柱体垂直破冰数值模拟冰层破裂过程(续)Fig.13 Numerical simulation of level ice rupture process by cylinder vertical ice breaking (continued)

图14 为圆柱体垂直破冰时数值模拟与试验结果的冰载荷时历曲线.在圆柱体接触冰层且匀速上升过程中，数值模拟冰载荷数据呈震荡上升的趋势，伴随着冰层裂纹产生和扩展，整体趋势、载荷峰值且突破时间均与试验值一致.

图14 圆柱体垂直破冰载荷数值与试验对比曲线Fig.14 Numerical and experimental comparison curve of vertical ice breaking load of cylinder

3.4 基于S-ALE 方法的楔形体入水数值模拟

为验证使用S-ALE 方法分析流固耦合问题的准确性与可行性，应用此方法模拟了二维楔形体入水砰击过程，并与陈光茂等[30]基于SPH 方法模拟二维楔形体入水问题的数值结果及Zhao 和Faltinsen[31]的实验结果进行对比.

依据文献[30]建立二维楔形体入水数值模型，水池长2.4 m，高1.3 m，以此保证入水砰击的过程中壁面反射不会对楔形体受力产生干扰.楔形体形状为等腰三角形，顶边长度为0.5 m，斜升角(即斜边与水平方向的夹角)为30°.在楔形体入水前的瞬间，它的垂向速度为6.15 m/s，随后自由入水，计算使用的二维楔形体为实心刚体，重量为 241 kg.

图15 为S-ALE 方法模拟的楔形体入水过程自由表面变形，与文献[30]试验现象一致，通过图16入水时间0.015 8 s 时压力分布对比图可知S-ALE 方法数值模拟结果与文献[30]中应用SPH 方法模拟结果相同，且图17 中入水垂向速度随时间变化对比曲线吻合性良好，证明了S-ALE 方法处理流固耦合问题的准确性与合理性.

图15 楔形体入水过程自由表面变形Fig.15 Deformation of free surface of wedge during water entry

图16 入水时间0.015 8 s 时压力分布对比图Fig.16 Comparison diagram of pressure distribution when water entry time is 0.015 8 s

图16 入水时间0.015 8 s 时压力分布对比图(续)Fig.16 Comparison diagram of pressure distribution when water entry time is 0.015 8 s (continued)

图17 入水垂向速度随时间变化图Fig.17 Variation diagram of vertical velocity of water inflow with time

4 圆柱体垂直出水破冰过程计算结果分析

圆柱体垂直出水破冰数值模拟有限元模型如2.2 节所示，模型冰材料选取及模拟工况均与无水环境下圆柱体垂直贯穿冰层一致.

4.1 冰载荷分析

圆柱体垂直出水破冰和无水环境下圆柱体垂直贯穿冰层冰载荷时历对比曲线如图18，出水破冰时圆柱体所受冰载荷整体趋势与无水破冰时类似，冰载荷峰值基本一致;出水破冰时冰层与圆柱体相互作用时间比无水破冰时的作用时间长，原因是在圆柱体与冰层接触之前的近场逼近过程中，结构体与冰层之间的水由于受到结构物与被约束的冰层的双重挤压而预先产生的一个高压力场，该压力场使船-冰之间产生一个对冰层存在支撑作用的阻力，即水垫效应[32].水垫效应作用在冰层，使的冰层弹性变形阶段更长，且冰层失效后未塌陷，仍然与圆柱体接触产生作用力;出水破冰时圆柱体上浮过程中水会上升挤压冰面，高压力场导致圆柱体尚未与冰层接触前发生弯曲变形从而产生挠度，所以有水环境下圆柱体与冰层的冰载荷要略晚于无水环境出现(0.08 s)，且有水情况下的冰载荷变化趋势相对平缓.

图18 有水和无水环境下圆柱体上浮破冰冰载荷时历曲线Fig.18 Time history curve of ice breaking load on floating cylinder in water and waterless environment

4.2 冰层挠度变化分析

圆柱体垂直出水破冰过程可分为水下上浮阶段和出水破冰阶段两部分.为研究圆柱体出水破冰过程中冰层挠度变化趋势，分别选取冰层上表面中心点A和Y轴四等分点B和X轴四等分点C观测挠度变化曲线，冰层挠度测量点A，B及C如图19 所示.

图19 冰层挠度测量点Fig.19 Measuring point of level ice deflection

圆柱体水下上浮阶段未与冰层接触时A，B，C点的挠度变化曲线如图20.可以看出A，B，C，3 点挠度变化趋势及斜率基本一致，挠度变化曲线均出现负值，呈现震荡后趋于稳定;A点挠度变化随圆柱体上浮过程整体大于B和C两点的挠度变化，说明水下上浮阶段圆柱体运动过程中引起水上升挤压冰面，导致圆柱体尚未与冰层接触前发生变形从而产生挠度变化，且应力集中在冰面中心点;B和C两点挠度变化曲线几近吻合，说明水下上浮阶段圆柱体对冰层的影响范围以冰面中心点向四周扩散且逐渐衰减，与冰层放置于水面之上且四周约束的工况设定一致.

图20 圆柱体水下上浮阶段冰层各点挠度Fig.20 Deflection of each point of level ice during underwater floating of cylinder

由图20 知，圆柱体出水破冰过程中冰层A点挠度随上浮时间变化最为明显，因此选取A点对比无水破冰和出水破冰时冰层的挠度变化.

图21 为无水环境和有水环境下圆柱体上浮破冰过程中冰层A点处的挠度变化曲线.由图21 可知两条挠度变化曲线整体趋势一致，随着圆柱体逐渐上升挤压冰层，冰层挠度逐渐增大;出水破冰过程冰层挠度变化峰值为1.48 mm，高于无水破冰时冰层挠度变化峰值0.81 mm;圆柱体突破冰层的时间也有明显差异，有水环境下的冰层突破时间远大于无水环境下的冰层突破时间;在2.5～ 2.58 s 间有一段平滑阶段，由于结构物与冰层挤压水形成的高压力场，使得在圆柱体与冰层接触之前冰层就存在初始挠度，即圆柱体在2.5～ 2.58 s 间未接触到冰层.图21 与图18 曲线表现一致，均证明了有水环境下结构物-冰层耦合作用时存在水垫效应.

图21 有水环境和无水环境下圆柱体上浮破冰过程中冰层A 点挠度变化Fig.21 Deflection change of point A of level ice in the process of cylinder floating and breaking ice in water environment and waterless environment

4.3 冰层应力响应及断裂现象分析

圆柱体出水破冰过程中冰层未发生断裂时von Mises 等效应力云图如图22.在圆柱体水下上浮阶段，圆柱体运动排开水挤压冰层，对冰层产生影响，应力主要集中于冰层中心即结构物正上方，但也存在应力并不集中的情况，由于结构物处于短暂加速状态会使应力值不稳定，如图22(a)到图22(b)过程;在圆柱体处于匀速上升阶段，应力完全集中在冰层中心，且冰层中心应力随着结构物运动时间越来越大，如图22(c)到图22(d)过程.在结构物出水破冰阶段，冰层应力极值点随着圆柱体的运动迅速增大，当t=2.58 s 时结构物接触冰层，应力极值为542.9 kPa，如图22(e);当t=2.94 s 时结构物即将突破冰层，应力极值达到最大值，冰层突破应力为2509 kPa，如图22(f).圆柱体出水破冰过程中的应力响应云图与冰层挠度变化趋势表现一致，无水环境下圆柱体突破冰层应力极值为2381 kPa，小于有水环境下冰层的突破应力极值.

图22 冰层未发生断裂时von Mises 等效应力云图Fig.22 von Mises equivalent stress diagram when the level ice does not break

图22 冰层未发生断裂时von Mises 等效应力云图 (续)Fig.22 von Mises equivalent stress diagram when the level ice does not break (continued)

图23 为圆柱体出水破冰过程中冰层发生断裂时的von Mises 等效应力云图及破裂现象.圆柱体出水破冰过程中冰层断裂现象大体上整体上与无水破冰相似，可以分为3 个阶段.第1 阶段，冰层受到圆柱体挤压达到突破载荷超过冰的抗拉强度，浮冰板生成十字径向裂纹，如图23(a)所示，此时发生轻微卸载，径向裂纹周围的应力小于圆柱体的突破应力，随着圆柱体运动十字径向裂纹持续扩展，径向裂纹周围的应力也逐渐增大，如图23(b)所示;第2 阶段，冰层在出现十字径向裂纹之后，随着圆柱体继续上浮，冰层上萌生更多径向裂纹形成楔子，此时冰层仍未到达失效极限，如图23(c)所示;第3 阶段，随着圆柱体继续上升，达到冰层的失效极限，由于在冰层的上、下表面已经存在大量裂纹，圆柱体突破冰层造成冰层失效，部分碎冰在水浮力的作用下沿圆柱体堆积形成碎冰堆，如图23(d)所示.

图23 冰层发生断裂时von Mises 等效应力云图及现象Fig.23 von Mises equivalent stress diagram and phenomenon when level ice breaks

图23 冰层发生断裂时von Mises 等效应力云图及现象(续)Fig.23 von Mises equivalent stress diagram and phenomenon when level ice breaks (continued)

5 结论

本文针对结构物与冰层相互作用强非线性问题，应用LS-DYNA 有限元软件建立了基于S-ALE流固耦合方法及罚函数接触算法的冰-水-结构物耦合作用数值模拟方法，研究了圆柱体垂直出水破冰过程，并对比分析了无水环境及有水环境下圆柱体垂直贯穿冰层破裂的共同点与区别，得到如下结论.

(1)通过冰梁结构4 点弯曲数值模拟，对比文献中的试验结果及数值模拟现象，验证了弹塑性应变率材料的冰力学特性.进行了无水环境下圆柱体垂直贯穿冰层模型试验及数值模拟，对比破冰载荷和冰层破裂现象，验证了有限元方法研究结构物-冰层耦合作用问题的准确性.基于S-ALE 方法模拟了楔形体入水流固耦合过程，验证了流固耦合算法的可行性.

(2)开展了圆柱体垂直出水破冰数值模拟，通过与无水破冰对比可知:有水环境下结构物-冰层间作用存在水垫效应.冰层突破载荷极值大小基本不会因有水环境或无水环境而产生明显变化.由于水对冰层的支撑作用，有水环境下的结构物突破冰层的冰载荷持续时间显著长于无水环境下冰载荷持续时间.

(3)描述了圆柱体垂直出水过程中的冰层断裂现象，与无水破冰相比裂纹萌生到扩展现象相似，最终碎冰失效形式不同.无水环境下碎冰会迅速塌陷，而有水环境下的碎冰会沿圆柱体周围堆积.有水环境下冰层弹性变形阶段更长，挠度变化大于无水环境下的挠度变化.

综上所述，由于在有水、无水环境下的结构物突破冰层冰载荷极值大小基本一致，研究冰层突破载荷极值及结构物强度相关问题时均可在无水环境下研究提高计算效率.而由于水垫效应的存在有水环境下结构物垂直出水贯穿冰层破裂过程与无水环境下结构物垂直贯穿冰层破裂过程存在明显区别，需要根据研究目标选择合适的数值模拟方法.