基于UMAT的冰—结构相互作用数值仿真

吕保达，刘敬喜，解　德，龚榆峰

基于UMAT的冰—结构相互作用数值仿真

吕保达，刘敬喜，解德，龚榆峰

华中科技大学船舶与海洋工程学院，湖北武汉430074

基于ABAQUS用户自定义子程序UMAT，采用已有的冰的弹性失效准则，针对冰与结构的相互作用进行数值模拟。首先，通过简单算例验证本文开发的UMAT的正确性；然后，采用开发的UMAT建立单个冰单元受刚性板挤压下的力学模型，得到力—位移曲线；最后，通过建立圆柱与冰相互作用的仿真分析模型，探讨冰板厚度对结构反力的影响。研究结果表明，采用用户自定义程序UMAT使得模拟更为复杂的冰力学特性成为可能，可为以后开展大规模的冰—结构相互作用研究奠定基础。

冰载荷；屈服面；弹性失效准则；冰—结构相互作用

网络出版地址：http://www.cnki.net/kcms/detail/42.1755.TJ.20150128.1201.002.html

期刊网址：www.ship-research.com

引用格式：吕保达，刘敬喜，解德，等.基于UMAT的冰—结构相互作用数值仿真［J］.中国舰船研究，2015，10（1）：39-45. LV Baoda，LIU Jingxi，XIE De，et al.Numerical simulation of ice-structure interaction based on UMAT［J］.Chinese Journal of Ship Research，2015，10（1）：39-45.

0　引言

北极具有丰富的渔业和矿物资源亟待开采，一旦开通北极航线，将会带来大量商用运输船舶与海洋平台及相关辅助工程船舶的需求。由于北极地区常年处于低温、多冰的环境，现有的常规船舶难以适应如此恶劣的环境，也很难应用于北极地区的矿产资源开采工作。因此，对现有船舶进行冰区加强改造，或者建造适用于极地作业的破冰船舶将成为北极开发必不可少的一个关键环节。

因此，开展冰与海洋工程结构相互作用的仿真分析研究，探讨冰及海洋工程结构物的损伤破坏机理，建立冰载荷作用下船舶及海洋工程结构的损伤评估方法，提出冰对船舶与海洋工程结构挤压及撞击载荷的计算方法，对预防因冰载荷而造成的海洋工程结构灾难性事故、提高海洋工程结构的安全性和可靠性具有重要的理论价值和社会意义。

国内外对冰与结构物的相互作用研究，称之为结构型冰力学。结构型冰力学主要探讨冰与结构物相互作用的力学机理。Sand［1］系统地研究了冰的力学性能，对相同失效准则下冰的力学性能进行了比较分析，并利用有限元仿真技术研究了不同形式的冰与海洋工程（简称“海工”）结构物的相互作用。Liu［2］对船舶与冰山碰撞的外部动力学和内部动力学问题进行了研究，并利用有限元软件进行了仿真分析。Karna等［3］利用子结构技术和有限元技术对浮冰挤压竖直海工结构的试验进行了仿真分析。Lu等［4］利用试验和解析方法研究了宽斜坡结构与平整冰的相互作用。Riahi等［5］利用试验和数值方法分析了冰与铝板相互作用时的力学行为。

国内针对冰与结构相互作用机理的研究起步较晚，某大学于1987年建成了我国第一座长5.5 m、宽1.9 m、深0.8 m的小型冰池实验室，开展了海冰与结构物相互作用的模型试验研究。史庆增等［6-10］利用物理模拟手段对冰与单桩、平台和斜面结构物的相互作用进行了研究。徐继祖［11-12］和王翎羽等［13］开展了海冰与窄体结构的动力作用研究。

因冰本身的力学特性不仅复杂，且差异很大，故给仿真模拟带来了很大的困难。而大多数仿真分析软件又都缺少对冰的材料属性的模拟，有关冰的仿真分析主要是借用其他类似材料的属性进行，这在很大程度上限制了分析的准确性。本文将采用用户自定义程序UMAT，开发针对冰弹性失效的材料模型程序，以期能够更好地利用相关试验研究成果，从而提高冰与其他结构相互作用数值模拟的准确性。

1　冰的力学特性研究

总体而言，冰是一种非均匀的各向异性的非线性材料。海冰的力学特性可以分为脆性与韧性2种。研究海冰的脆性主要考虑海冰的压缩、拉伸、弯曲和剪切等强度，主要通过海冰的弹性模量进行研究，这些力学特性通常与冰粒的大小、晶体形式、温度和应变率等因素有关，而研究海冰的韧性则还需考虑时间因素。

到目前为止，在海冰的所有力学特性中，弹性是被研究得最为广泛的一种特性。研究海冰的弹性力学行为主要采用静态和动态2种方法。研究结果显示，通过静态方法获得的冰的弹性模量在0.3～10 GPa范围内变化，采用动态方法测量所得的弹性模量在6～10 GPa范围内变化。冰在静态条件下的屈服强度小于动态条件下的屈服强度。在静态条件下，由于冰的黏弹性，测量冰的总变形的静态分析方法并不能准确得到冰对瞬时弹性变形的抵抗力，应力越大，冰的黏弹性和蠕变行为的影响就会越明显。因此，通过测量应变法得到的冰的弹性模量会随着应力的增大而减小。与此不同的是，动态测量方法主要受振动在冰试件中的传播速度和试件的固有频率的影响，可以忽略黏弹性和蠕变行为的影响。

1.1冰的弹性失效准则分析

海冰主要存在2种晶体结构形式：C轴随机分布粒状冰和C轴平行分布的柱状冰。粒状冰为各向同性，柱状冰在垂直C轴的平面内为各向同性，在平行C轴的平面内为正交各向异性。

要准确描述冰的材料特性目前仍然比较困难，本文提出将冰的应力是否满足失效曲面作为评判冰是否失效的准则。对于各向同性粒状冰，Reinicke和Remer［14］提出了一个椭圆的失效曲面方程：

式中：I1为第1应力张量不变量；J2为第2应力偏张量不变量；k1，k2，k3为材料常数，须通过试验获得。该准则的第1应力张量不变量I1是线性项，能够预测压缩和拉伸2种不同模式的失效。而2阶偏应力不变量则确保在静水力平面上的强度按照椭圆屈服面准则变化。图1（a）给出了在静水力平面上的屈服面形状，图1（b）给出了屈服面在偏张量平面上的投影。需注意的是，该准则在偏平面上的形状为圆形，类似于Von Mises失效准则，剪切变形能控制着失效。与其他准则相比，该准则能够更好地与各向同性的冰的试验数据相吻合。

图1　式（1）定义的三参数屈服面Fig.1　Three-parameter yield surface defined by Eq.（1）

Horrigmoe和Zeng［15］利用式（1）给出的Reinicke和Remer的失效准则，对冰盖的压缩试验进行了仿真分析，获得了准则中的材料常数，其中，k1=0.674 MPa-2，k2=0.853 MPa-2，k3=0.014 MPa-2。

Horrigmoe和Zeng［15］针对正交各向异性材料提出了如下失效准则：

式中：σx，σy，σz分别为x，y，z方向的正应力；τxy，τxz，τyz分别为xy，xz，yz平面内的切应力；由12组独立试验确定的材料常数a1，…，a12分别是：

1）3组xy，yz，xz平面内的剪切试验；

2）3组x，y，z方向上的简单压缩试验；

3）3组x，y，z方向上的简单拉伸试验；

4）3组x，y，z方向上的受限压缩试验。

对于横观各向同性，z轴作为旋转对称轴，独立材料常数由12个减少为7个：a1=a2，a5=a6，a7=a8，a10=2a1-a4，a11=a12。

则式（2）变为

对于横观各向同性，xz和yz平面内的剪切强度Sxz，Syz与xy平面内的剪切强度Sxy接近，若假定三者相等，则有

则式（3）可以写为

式中：a1=0.768，a3=0.108，a5=-0.117，a7=2.533，a9=0.672。

对于Horrigmoe和Zeng的失效准则（4），式（4），如果

则Horrigmoe和Zeng的失效准则（4），式（4），与Reinicke和Remer的失效准则（1），式（1），完全相同。

1.2冰压碎研究的弹性方法

在利用弹性方法研究冰受到挤压之后的响应时，冰的总应变为弹性应变εeij和压碎后的应变之和，即

在冰受到拉伸时，冰的总应变为

在本文的研究中，冰的应力达到失效准则后失效，冰在失效前后均保持弹性，弹性模量由失效前弹性模量E变为失效后弹性模量Ec。在失效前，对于横观各向同性冰，应力与应变的关系为

式中，

其中：Ex，Ey，Ez分别为x，y，z方向的弹性模量；νxy，νzx为泊松比；Gxy，Gzx为剪切模量；νxy，νyz，νxz为剪切应变。

对于各向同性冰，应力与应变的关系变为

当冰中的应力达到失效准则时，弹性模量由失效前弹性模量E变为失效后弹性模量Ec，泊松比由ν变为νc，应力与应变的关系为

由于失效后的弹性模量远小于失效前的弹性模量，所以，失效后，单元的刚度和失效单元的应力σ=Ecε迅速下降，如图2所示。图中，Tg为冰的拉伸强度，Cg为冰的压缩强度。

图2　冰的单轴拉伸压缩应力—应变曲线Fig.2　Uniaxial stress-strain curve of tension and compression

2　算例分析

2.1用户自定义材料子程序UMAT

在有限元分析软件ABAQUS中没有失效准则（1）和（2）的材料模型，本文通过ABAQUS用户自定义材料子程序UMAT，利用失效准则（1）模拟C轴随机分布的各向同性粒状冰，利用失效准则（4）模拟C轴平行分布的横观各向同性柱状冰。

ABAQUS中的用户自定义材料子程序UMAT可以定义材料的本构关系，使用ABAQUS材料库中没有包含的材料进行计算，扩充程序功能。

在本文定义的材料子程序中，分别利用失效准则，以状态变量的形式对冰的状态进行判断，当状态变量大于1时，材料失效。

在UMAT中，利用式（7）与式（8）定义材料的本构模型，即应力增量对应变增量的变化率。当材料达到失效状态后，通过减小弹性模量的方式处理材料的失效，即弹性模量由失效前的E变为失效后的Ec。

根据1.1节中Reinicke和Remer提出的失效准则（1）及1.2中的弹性方法，针对各向同性冰粒状冰，编制用户自定义材料UMAT，取参数如下：

失效之前，冰的弹性模量为6 GPa，失效之后为20 MPa，失效前、后泊松比不变。

根据Horrigmoe和Zeng的失效准则（4），编制用户自定义材料子程序UMAT，模拟横观各向同性的柱状冰，取参数如下：

失效之前，平行于C轴方向的弹性模量为8.2 GPa，横向弹性模量为6 GPa，失效之后，均为20 MPa，失效前、后冰的泊松比不变。

2.2单个单元受挤压响应验证

为了验证UMAT及算法的正确性，首先建立边长为1 mm的单个立方单元体，利用失效准则（1）对冰体进行仿真分析，如图3所示。

图3　单个冰单元受压仿真分析模型Fig.3　Ice cube modeled with a single C3D8 element

分别对单元上面4个节点施加z方向的位移，模拟各向同性冰受拉伸和压缩2种情况，得到z方向的应力—应变曲线，并与文献［1］中提供的单元应力—应变曲线进行比较，如图4所示。

从图4可以看出，本文计算得到的单元应力—应变曲线与文献［1］中计算所得的应力—应变曲线吻合良好，说明了本文开发的用户自定义UMAT的正确性。

为进一步分析各向同性冰在受到挤压时的力学行为，建立冰与刚性平板相互作用的有限元模型，对冰体受刚性平板挤压后的力学行为进行仿真分析，单个单元的模型如图5所示。

图4　单个单元z方向应力—应变曲线比较Fig.4　Stress-strain curves of uniaxial tensile and compressive in the z direction of the single element and the curves in the reference［1］

图5　单个单元受刚性板挤压仿真分析模型Fig.5　Single C3D8 element pressed by a rigid plate

利用失效准则（1）进行计算分析，得到刚性板所受的反力如图6所示。

图6　刚性板所受反力Fig.6　The reactive force on the rigid plate

从图6可以看出，刚性板在挤压冰时受到的反力的变化趋势与图4中单个单元施加位移时应力—应变曲线的变化趋势相一致。

对于C轴平行分布的横观各向同性柱状冰，建立相同模型，利用刚性板对冰体在平行C轴方向与垂直C轴方向进行挤压，得到刚性板反力如图7所示。

图7　刚性板所受反力Fig.7　The reactive force on the rigid plate

从图7可以看出，反力变化趋势与图2一致。在利用刚性板在平行C轴方向对冰单元挤压时的最大反力是8.21 N，垂直C轴方向时是3.88 N，平行C轴方向明显高于垂直C轴方向。

2.3圆柱挤压冰板仿真分析

在北极环境中，冰盖是冰的一种主要形式，对冰盖受挤压的力学行为进行仿真分析对于冰区船舶与海洋工程结构物的设计建造具有重要意义。

建立面积为1 000 mm2×1 000 mm2，厚度分别为10，20，30和40 mm的冰板，刚性圆柱曲面半径为100 mm，对冰板进行挤压模拟，如图8所示。网格大小为10 mm3×10 mm3×10 mm3，分析不同厚度冰板对刚性圆柱曲面挤压时的力学行为。

图8　冰与刚性圆柱曲面挤压时的仿真模型Fig.8　Ice compressed by a rigid cylindrical surface

图9给出了当刚性圆柱全面位移为不同厚度的冰板受到刚性圆柱面挤压后的失效云图，冰载荷受到圆柱挤压时为简单的拉压破坏。不同厚度的冰板失效区域大致相同。图10为刚性圆柱曲面上的反力曲线。从图中可以看出，随着冰板厚度的增加，反力呈线性增加趋势，说明冰的厚度对海上结构物的影响极为重要。在设计和建造冰区的海上结构物时，需要考虑所处区域冰厚度的影响，以免海上结构物受到损伤。

图9　厚度为10，20，30和40 mm的冰的失效云图Fig.9　Ice failure contours for different thickness of 10，20，30 and 40 mm

图10　不同厚度的冰对刚性圆柱曲面反力Fig.10　The reactive force on the rigid cylindrical surface

3　结语

本文基于ABAQUS的自定义用户程序UMAT，利用Reinicke和Remer提出的失效准则，开发出了冰材料UMAT程序，对冰在受到挤压时的力学行为进行仿真，获得了冰对结构的反力。与弹性方法的拉伸时所受反力的变化趋势相比，两者在变化趋势上相一致，说明本文提出的方法可用于冰的力学行为的仿真分析。

由于冰材料的复杂性，对冰本身材料的力学行为的仿真模拟难度很大。本文提出的方法可以简化材料仿真模拟的难度，并通过与其他方法相结合，可以进一步提高仿真分析的准确性，以更好地模拟冰载荷，促进冰区船舶与海洋工程结构物设计水平的提升。

［1］SAND B.Nonlinear finite element simulations of ice forces on offshore structures［D/OL］.Lulea，Sweden：Lulea University of Technology，2008.［2014-01-15］. http：//epubl.ltu.se/1402-1544/2008/39/index.html.

［2］LIU Z H.Analytical and numerical analysis of iceberg collisions with ship structures［D］.Norway：Norwegian University of Science and Technology，2011：235.

［3］KARNA T，KAMESAKI K，TSUKUDA H.A numeri⁃cal model for dynamic ice-structure interaction［J］. Computers and Structures，1999，72（4/5）：645-658.

［4］LU W J，LUBBAD R，HØYLAND K，et al.Physical model and theoretical model study of level ice and wide sloping structure interactions［J］.Cold Regions Sci⁃ence and Technology，2014，101：40-72.

［5］RIAHI M M，MARCEAU T D，LAFORTE C，et al. The experimental/numerical study to predict mechani⁃cal behavior at the ice/aluminum interface［J］.Cold Re⁃gions Science and Technology，2011，65（1）：191-202.

［6］史庆增.作用于直立圆柱上的横向冰力［J］.海洋学报，1994，16（4）：126-129.

［7］史庆增，宋安.海冰静力作用的特点及几种典型结构的冰力模型试验［J］.海洋学报，1994，16（6）：133-141.

［8］史庆增，王永安.辽东湾孤立桩柱上冰力的概率分布［J］.海洋学报，1995，17（4）：130-136.

［9］史庆增，徐继祖，宋安.海冰作用力的模拟实验［J］.海洋工程，1991，9（l）：16-22. SHI Qingzeng，XU Jizu，SONG An.The modeling test of the sea ice forces［J］.The Ocean Engineering，1991，9（1）：16-22.

［10］史庆增.带有斜面的结构物上冰力的试验研究［C］//第六届全国近海工程学术会论文集（下册），1992：496-500.

［11］徐继祖.海冰引起的结构振动［J］.海洋工程，1986， 4（2）：42-47.

［12］徐继祖.对渤海海冰工程中几个问题的认识［J］.中国海上油气（工程），1993，5（4）：13-17. XU Jizu.Some views on Bohai sea ice engineering［J］. China Offshore Oil and Gas Engineering，1993，5（4）：13-17.

［13］王翎羽，徐继祖.冰与结构动力相互作用的理论分析模型［J］.海洋学报，1993，15（3）：140-146.

［14］REINICKE，K M，REMER R.A procedure for the determination of ice forces-illustrated for polycrys⁃talline ice［C］//Proceedings of the 5th IAHR Sympo⁃sium on Ice.Lulea，Sweden，1978：217-238.

［15］HORRIGMOE G，ZENG L F.Modelling ductile be⁃havior of columnar ice using computational plasticity［C］//Proceedings of the 12th IAHR Symposium on Ice.Trondheim，Norway，1994，1：282-291.

［责任编辑：喻菁］

Numerical Simulation of Ice-Structure Interaction Based on UMAT

LV Baoda，LIU Jingxi，XIE De，GONG Yufeng
School of Naval Architecture and Ocean Engineering，Huazhong University of Science and Technology，Wuhan 430074，China

In this paper，the failure criterion of ice is adopted to define the user-defined material（UMAT）of Abaqus.The UMAT is first proved to be correct through a simple simulation example.Then，the simula⁃tion of a single ice element pressed by a rigid plate is performed，and the reaction curve of force-displace⁃ment is obtained through simulation.Finally，by simulating the interaction between the ice and a rigid cylin⁃der surface，the relation between the thickness of ice and the reaction force is analyzed，which is shown to be linear.The results presented in this paper indicate that it is feasible to simulate the mechanical proper⁃ties of ice by the user-defined material，and this finding provides valuable references to future researches concerning large-scale ice-structure interaction.

ice load；yield surface；elastic failure criterion；ice-structure interaction

U661.4

A

10.3969/j.issn.1673-3185.2015.01.006

2014-06-14

网络出版时间：2015-1-28 12:01

国家高技术研究发展计划（863计划）资助项目（2012AA112601）

吕保达，男，1990年生，硕士生。研究方向：船体结构。E⁃mail：lvbaoda@163.com

刘敬喜（通信作者），男，1975年生，博士，副教授。研究方向：船体结构。E⁃mail：liu_jing_xi@mail.hust.edu.cn