基于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