基于ABAQUS的海冰單元開發及冰載荷直接計算法

龔榆峰，張正藝，2，3，劉敬喜，2，3，董問，解德，2，3

1華中科技大學船舶與海洋工程學院，湖北武漢430074

2船舶和海洋水動力湖北省重點實驗室，湖北武漢430074

3高新船舶與深海開發裝備協同創新中心，上海200240

基于ABAQUS的海冰單元開發及冰載荷直接計算法

龔榆峰1，張正藝1，2，3，劉敬喜1，2，3，董問1，解德1，2，3

1華中科技大學船舶與海洋工程學院，湖北武漢430074

2船舶和海洋水動力湖北省重點實驗室，湖北武漢430074

3高新船舶與深海開發裝備協同創新中心，上海200240

［目的］隨著全球氣溫的不斷升高，北極冰層逐漸融化，開發利用北極航線和北極資源受到國際社會的廣泛關注，各國掀起了新一輪的極地船舶建造高潮，而在極地船舶的設計過程中，如何確定船舶冰載荷顯得尤為重要。近年來，采用數值方法計算極地結構受到的冰載荷逐漸成為國內外學者關注的重點，而目前的通用商業有限元軟件中缺少特定用于模擬海冰力學特性的相關海冰單元，為相關數值研究工作的開展增加了困難。針對這一問題，［方法］結合Reinicke和Remer提出的各向同性海冰失效準則，在ABAQUS二次開發平臺UEL上開發用于模擬海冰失效的新型單元。在所開發的單元中引入彈性基礎用以模擬浮力和重力的作用。采用典型的斜坡和圓錐結構算例并結合相關規范公式對所開發的新型單元進行驗證。［結果］對比二者之間的計算結果發現，直接計算法與規范公式的計算結果吻合較好。［結論］所提出的直接計算法能夠為冰區結構設計提供工具和一定的參考。

冰載荷；直接計算法；ABAQUS；沙漏控制；自定義單元

0 引 言

隨著全球氣溫的升高，北極地區的海冰逐漸融化，人類對極地船舶與海洋結構物的需求也隨之增加。在進行極地船舶與海洋結構物的設計時，如何合理、有效確定結構物與海冰相互作用過程中受到的冰力，是一個長期被國內外相關研究人員關注的問題。然而海冰與結構物的相互作用是一個相當復雜的過程，在目前的研究中，研究人員通常將海冰與不同結構物相互作用過程中海冰的破壞模式歸納為壓碎、彎折和劈裂等幾類。

當海冰與斜坡或圓錐結構物作用時海冰常常會出現彎折破壞。國內外相關研究人員針對海冰與海洋結構物的彎折破壞進行了大量研究工作，其中 Bercha和 Danys［1］，Croasdale等［2］將冰排理想化為放置于彈性基礎上的線彈性板，從而提出了計算冰排和斜坡、圓錐結構相互作用時冰力大小的理論公式。Ralston［3］則基于塑性上限理論提出了計算冰排和圓錐結構相互作用時冰力大小的理論公式。李峰等［4］采用彈性基礎梁理論分析了冰與斜面結構作用時的彎曲破壞。法國船級社（Bureau Veritas，BV）公布的針對船舶和海洋結構物上冰力大小的計算指南（BV_565NI）［5］和美國石油協會（American Petroleum Institute，API）公布的規范API-RP-2N［6］中將海冰理想化為作用于彈性基礎上的線彈性材料，提出了用于計算冰排和斜坡、圓錐結構相互作用時冰力大小的相關公式。

對于海冰彎折破壞的數值計算研究，Sand和Horrigmoe［7］，Derradji-Aouat［8］等采用商用有限元軟件ANASYS計算了冰排與斜坡結構相互作用時結構受到的冰力，在其研究中冰排被理想化為作用于彈性或彈塑性基礎上的線彈性板。黃焱等［9］在LS-DYNA中模擬了圓環形冰排彎曲破壞以及冰排與錐體結構相互作用的過程。宋安等［10］采用縮尺的冰模型試驗研究了冰排作用于斜面結構時的冰力大小，并將結果同相關理論公式進行了對比。

綜上所述，國內外研究人員已對冰排與固定傾斜結構作用時結構受到的冰載荷進行了大量的理論、試驗和數值研究。然而在已開展的冰排與固定傾斜結構作用時冰載荷的數值研究中，海冰本構關系的模擬大多基于商業有限元軟件已有的材料模型和失效準則，使得難以對海冰的力學特性進行合理描述。為了更好地評估冰排與固定傾斜結構作用時結構受到的冰載荷，應在數值分析中引入適用于模擬海冰力學特性的特定材料模型和失效準則。因此，進一步開展冰排與固定傾斜結構作用時的冰載荷直接計算法研究，并基于現有商業有限元軟件開發具有一定通用性的海冰單元以描述海冰的力學行為與考慮浮力的作用很有必要，這將進一步提高采用直接計算法確定結構冰載荷的效率，并使得所編制的程序具有一定的通用性，具有一定的工程實踐意義。

本文將基于通用商業有限元軟件ABAQUS的自定義單元子程序創建集成了Reinicke和Remer［11］提出的各向同性海冰失效準則和彈性基礎的8節點六面體單元用于模擬海冰。同時為了對所開發的單元進行驗證，將采用開發的單元計算冰排與斜坡結構相互作用過程中的冰力大小，并將所獲得的計算結果與BV和API給出的相關規范公式進行對比。

1 單元及失效準則介紹

1.1 所開發單元基本原理

8節點六面體縮減積分單元示意圖如圖1所示。

圖1 單元示意圖Fig.1 Schematic of element

根據有限元理論［12］，單元內任意一點的位移u可以表示為

式中：u，v，w為沿坐標軸x，y，z軸的位移分量；x，y，z為總體坐標系的3個坐標軸。單元內任意一點的位移u可以由單元的形函數和單元各個節點的位移來表示：

式中：Na為基于節點位置事先給定的函數，又稱單元的形函數；ua為單元內各個節點的位移。

單元任意一點的應變可用引起應變的3個分量來表示：

式中：S為單元應力矩陣；Bi為單元應變矩陣。

所以Bi可以由式（4）給出。

根據應力—應變關系，采用矩陣D可以將單元中的6個應力分量和6個應變分量聯系起來，單元中應力σ可以表示為

最后，由虛功原理可以得到單元剛度矩陣K，

式中，J為雅克比行列式。采用高斯積分時，式（6）可以簡化為

式中，H為權函數。由于本文所開發的單元為完全縮減積分單元，只有一個積分點，即式（7）中i=1。對應的積分點在局部坐標系下的坐標為（0，0，0），即式（7）中ξi=0，ηi=0，ζi=0，對應于積分點數量的權函數H=8.0。

為了考慮海冰受到的浮力，在單元中引入彈性基礎用于考慮海冰的浮力作用。彈性基礎反力R與單元底部節點的垂向位移w和彈性基礎中彈簧的剛度ks之間的關系如下［13］：

由節點力平衡方程可以得到彈性基礎反力和垂向位移之間的關系：

式中，K0為考慮彈性基礎而引入的單元剛度附加矩陣。

最后，可以得到引入彈性基礎后的單元平衡方程：

式中，Ke為不考慮彈性基礎時的單元剛度矩陣。

1.2 Reinicke和Remer各向同性海冰失效準則

Reinicke和 Remer［11］在總結前人研究的基礎上提出了用于描述各向同性海冰失效行為的失效準則。

式中：k1，k2和k3為材料常數，通過試驗確定；I1為應力張量第1不變量；J2為應力偏量第2不變量。當f(I1，J2，k1，k2，k3)=0 時，材料失效。失效后的材料彈性模量和泊松比分別變為Ec=60 MPa，vc=0.33。本文中使用的材料常數的取值參考文獻［9］中由模型試驗擬合得到的結果，k1=0.674（MPa）-2，k2=0.853（MPa）-1，k3=0.014（MPa）-2。

1.3 沙漏控制

本文采用啞單元方法并結合ABAQUS用戶自定義輸出變量子程序UVARM來引入單元的沙漏控制，實現所開發單元的結果可視化。啞單元是指與所開發單元使用相同節點的一組C3D8R實體單元，其彈性模量定義為一個非常小的值，如1.0×10-20Pa，同時可以在啞單元中直接定義單元的沙漏剛度，從而引入單元沙漏控制。通過UVARM子程序，將所開發單元的應力、應變和損傷變量等單元狀態變量賦值給啞單元上的自定義輸出變量，從而可以在ABAQUS后處理中實現計算結果的可視化。具體流程如圖2所示。

圖2 單元計算流程圖Fig.2 Flow chart of element computation

2 算例驗證

為了驗證所開發的單元，首先采用所開發的單元模擬正方體冰體單軸拉伸和壓縮情況，將計算得到的單元應力—應變關系與解析解結果進行對比。而后，采用所開發的單元計算3組算例：端部受垂向載荷的半無限長冰梁、半無限大冰排與斜坡結構、錐臺結構的相互作用。并將計算結果與 BV［5］以及 Croasdale［14］，Ralston［3］給出的建議公式進行比較，其中BV給出的一端受垂向載荷作用的半無限長冰排能承受的最大外力F的計算公式如式（12）所示，半無限長冰排失效位置計算公式如式（13）所示。計算中采用的材料參數及Reinicke和Remer失效準則參數如表1所示。

式中：B為冰排的寬度；h為冰排的厚度；σf為冰排的拉伸強度；E為冰排的彈性模量；ρw為海水的密度；g為重力加速度。

表1 計算中使用材料參數Table 1 The material constants used in the simulation

2.1 單元基本特性驗證

為了對所開發單元的基本特性進行驗證，計算了一頂部受到強制位移、底面約束的正方體冰塊的拉伸和壓縮，并將計算得到的單元應力—應變關系與解析解結果進行對比。算例的示意圖如圖3所示。算例相關幾何參數如下：模型的長度l=1 m，寬度w=1 m，厚度h=1 m。模型共有1個單元，8個節點。其中：頂部節點為5，6，7，8；底部節點為1，2，3，4。頂部節點施加z方向的強制位移；底部節點中節點1約束x，y，z方向的位移，節點2約束y，z方向的位移，節點3約束z方向的位移，節點4約束x，z方向的位移。

圖3 正方體冰塊單軸拉伸和壓縮示意圖Fig.3 Schematic of ice cube under uniaxial stretching and compression

使用新型單元得到的計算結果與Reinicke和Remer屈服準則解析解的對比如圖4所示。將新型單元的計算結果與解析解的對比可以發現，該新型單元不僅能夠描述海冰的拉伸與壓縮強度不一致這一基本力學特性，而且計算得到的拉伸強度和壓縮強度與Reinicke和Remer屈服準則解析解吻合很好，驗證了本文提出的新型單元能夠很好地實現Reinicke和Remer提出的海冰本構關系，并描述了海冰的基本力學特性。

圖4 應力—應變曲線計算結果對比Fig.4 Comparison of stress-strain curve

2.2 半無限長冰梁

為了驗證所開發的單元，計算了半無限長冰梁在端部受到垂向位移作用時的反力，并將所得計算結果與BV給出的建議公式進行了比較。一端受垂向載荷的半無限長冰梁模型如圖5所示。模型幾何尺寸如下：冰梁長度l=400 m，冰梁寬度w=1 m。為了研究冰梁厚度對冰力的影響，選取了不同厚度的冰梁進行計算。冰梁的厚度分別為h=0.5，1.0，1.5，2.0，2.5和3.0 m。冰梁模型共采用2 800個新型單元和2 800個啞單元。同時，由于冰梁主要發生彎折失效，為了保證直接計算法的計算精度并描述冰梁沿厚度方向的失效發展過程，將冰梁的厚度方向劃分為7個單元。模型邊界條件設置為：在冰梁模型邊界一端約束x，y方向的位移，另一端施加垂向的位移。

圖5 半無限長冰梁示意圖Fig.5 Schematic of semi-infinite ice beam

采用開發的新型單元模擬一端受垂向載荷的半無限長冰梁，得到不同時刻下厚度為0.5和3.0 m的冰梁損傷分布，如圖6和圖7所示。其中藍色部分表示未損傷單元，紅色表示已出現損傷的單元。從圖6和圖7可以看出，不同厚度的一端受到垂向載荷的半無限長冰梁其損傷發展隨著厚度的不同并無明顯區別。由于冰梁端部受到向下的垂直位移作用，因而冰梁在變形的過程中，上表面受到拉伸而下表面為壓縮變形，由于采用的本構模型中海冰的拉伸強度要明顯小于其壓縮強度，所以在模型中冰梁的損傷均首先出現在上表面。而后隨著端部位移的不斷增加，冰梁的損傷區域由出現初始損傷的上表面逐漸沿著厚度方向，從冰梁的上表面向下表面擴展，直至整個厚度方向全部失效。同時，隨著冰梁厚度的不同，冰梁出現初始損傷的位置到加載端的距離不斷增加。但當冰梁出現初始失效后，隨著端部位移的增加，冰梁的損傷區域不僅沿冰梁的厚度方向擴展，也會沿長度方向不斷向加載端靠近。

圖6 厚度為0.5 m的冰梁損傷分布Fig.6 Failure process of ice beam with thickness h=0.5 m

圖7 厚度為3.0 m的冰梁損傷分布Fig.7 Failure process of ice beam with thickness h=3.0 m

一端受到垂向載荷的不同厚度半無限長冰梁的冰力和端部位移之間的關系曲線如圖8所示。從圖8中可以看出，隨著垂向位移的逐漸增加，端部受到的反力也線性增加，當冰梁中的應力狀態達到單元失效準則所定義的失效面時，冰梁出現失效。當單元出現失效時，直立結構受到的冰力隨位移的增加出現明顯的下降。但由于彈性基礎的存在，以及首先出現失效的只是冰梁上表面的部分單元，冰梁并未全部失效，故仍能繼續承受載荷。所以加載端部受到的反力隨位移的增加，出現下降后又會再次繼續上升。

圖8 不同厚度的半無限冰梁端部的反力和垂向位移關系曲線Fig.8 Vertical force-displacement curves of semi-infinite ice beam with different thickness

為了驗證采用新型單元獲得的計算結果，將采用直接分析法獲得的冰梁出現初始失效時的失效位置、端部承受的反力與BV［5］給出的理論公式獲得的解析解進行了對比，如圖9和圖10所示。從圖中的結果對比可以看出，通過新型單元獲得的失效位置和對應的端部反力與解析解吻合良好，誤差僅為5%左右。說明所開發的新型單元能夠準確模擬一端受垂向載荷的半無限長冰梁的彎折失效。

圖9 直接計算法與BV公式端部反力結果對比Fig.9 Comparison of vertical force between direct analysis method and the equation presented by Bureau Veritas

圖10 直接計算法與BV公式失效位置結果對比Fig.10 Comparison of failure position between direct analysis method and the equation presented by Bureau Veritas

2.3 冰排與斜坡結構的相互作用

為了驗證所開發的單元，計算了斜坡結構與冰排作用時受到的冰力，并將所得計算結果與Croasdale［14］提出的理論公式進行了對比。模型示意圖如圖11所示。

圖11 斜坡結構示意圖Fig.11 Schematic of sloping structure

模型幾何尺寸如下：冰排長度l=400 m，冰排寬度w=1 m，冰排厚度h=3 m。為了研究斜坡結構的傾角對斜坡結構受到的冰力的影響，選取了不同傾角的斜坡結構進行計算。斜坡結構的傾角分別為α=30°，35°，40°，45°，50°，55°和 60°。冰排模型共采用2 800個新型單元和2 800個啞單元，矩形直立結構共采用6個R3D4剛體單元。模型邊界條件設置為：在冰排模型邊界AB上約束x，y方向的位移，其他邊界不施加約束。同時直立結構的約束方式為：約束除x方向平動以外的所有自由度，同時施加x方向的平動位移。

采用開發的新型單元模擬冰排與斜坡結構的相互作用，得到不同時刻下與傾角α=30°和60°的斜坡結構作用時的冰排損傷分布圖，如圖12和圖13所示。其中藍色部分表示未損傷單元，紅色表示已出現損傷的單元。從圖12和圖13中可以看出，與傾角α=30°和60°的斜坡結構作用時，冰排的損傷模式并未有明顯的差別。由于冰排的拉伸強度要明顯小于冰排的壓縮強度，所以冰排與斜坡結構相互作用時，冰排的損傷均是首先出現在冰排受到拉伸的下表面處，然后損傷區域逐漸沿著冰排的軸向和厚度方向發展。

圖12 斜坡結構傾角α=30°時冰排損傷分布Fig.12 Failure process of sloping structure with slope angleα=30°

圖13 斜坡結構傾角α=60°時冰排損傷分布Fig.13 Failure process of sloping structure with slope angleα=60°

冰排與不同傾角的斜坡結構作用時，斜坡結構受到的水平和垂直冰力與冰排的水平位移之間的關系如圖14和圖15所示。從圖14和圖15中可以看出，隨著冰排水平位移的逐漸增加，斜坡結構受到的水平和垂直冰力也線性增加。當冰排中的應力狀態達到單元失效準則所定義的失效面時，冰排出現失效。斜坡結構受到的水平和垂直冰力隨冰排水平位移的增加而下降。隨著冰排水平位移的繼續增加，冰排沿著斜坡結構的斜面爬升，并且由于彈性基礎的存在，導致斜坡結構受到的冰力再次隨位移的增加而增加，使得斜坡結構的受力再次繼續上升。

圖14 斜坡結構受到的水平冰力與位移之間的關系Fig.14 Relationship between horizontal force and displacement of sloping structure

圖15 斜坡結構受到的垂直冰力與水平位移之間的關系Fig.15 Relationship between vertical force and displacement of sloping structure

對于不同傾角的斜坡結構，隨著斜坡結構傾角的增加，冰排出現初始失效時對應的斜坡結構受到的水平冰力也隨之增加。斜坡結構傾角為30°時對應的冰排失效時結構受到的水平冰力為5.75×104N，而隨著斜坡結構傾角的增加，當斜坡結構傾角為60°時對應的水平冰力達到2.02×105N，比傾角為30°時對應的水平冰力增加了351%。但冰排出現初始失效時對應的斜坡結構受到的垂直冰力基本不隨斜坡結構傾角的變化而變化。斜坡結構傾角為30°時對應的冰排失效時結構受到的垂直冰力為9.96×104N，而隨著斜坡結構傾角的增加，當斜坡結構傾角為60°時對應的垂直冰力為1.16×105N，相對傾角為30°時對應的垂直冰力沒有明顯的變化，只略微增加了16%。這主要是由于冰排在與斜坡結構作用時因為冰排發生彎曲變形而出現了彎折失效。而冰排的彎折失效主要由冰排受到的垂向分力引起，因為不同算例中冰排的幾何尺寸和材料屬性并未改變，所以其發生彎折失效所對應冰排受到的垂向分力基本保持不變。但隨著斜坡角度的增加，冰排受到的水平分力隨之增加，所以導致冰排出現初始失效時對應的斜坡結構受到的水平冰力隨斜坡結構傾角的增加而增加，而垂直冰力基本不隨斜坡結構傾角的變化而變化。

為了驗證采用新型單元獲得的計算結果，將采用直接分析法獲得的冰排出現初始失效時的斜坡結構受到的水平和垂直冰力與采用Croasdale［14］提出的理論公式獲得的解析解進行了對比，如圖16和圖17所示。從圖16中的結果對比可以看出，通過新型單元獲得的冰排出現初始失效時斜坡結構受到的水平和垂直冰力與解析解吻合良好。但隨著斜坡結構傾角的增加，直接計算法相對于Croasdale［14］的理論公式結果之間的誤差有逐漸增加的趨勢，在結構傾角為60°時，最大誤差為15%。引起誤差的原因主要為冰排與斜坡結構作用時，在Croasdale［14］的理論公式中只考慮了冰排的彎曲，而忽略了斜坡結構對冰排沿水平方向存在的面內力，而在采用新型單元進行數值模擬時，斜坡結構存在對冰排的面內力，從而導致采用新型單元得到的冰力與Croasdale［14］的理論公式給出的冰力存在一定的誤差。

圖16 直接計算法與Croasdale［14］公式水平方向冰力結果對比Fig.16 Comparison of horizontal force between direct analysis method and the equation presented by Croasdale［14］

圖17 直接計算法與Croasdale［14］公式垂向冰力結果對比Fig.17 Comparison of vertical force between direct analysis method and the equation presented by Croasdale［14］

2.4 冰排與圓錐結構的相互作用

計算圓錐結構與冰排作用時受到的冰力，并將所得計算結果與Ralston［3］給出的塑性極限分析解進行了比較，模型示意圖如圖18所示。計算圓錐結構的角度為30°～60°，間隔5°，圓錐結構的水線面直徑D=12 m，圓錐結構水線面以上部分高度ht=3 m，水線面以下部分高度hb=6 m；冰排的長度l=100 m，寬度w=50 m，厚度h=0.8 m。建立的有限元模型共有41 760個單元，其中厚度方向有5個單元。

圖18 圓錐結構示意圖Fig.18 Schematic of cone structure

采用開發的新型單元模擬冰排與錐體結構的相互作用，得到的不同時刻下與錐角α=30°和60°的錐體結構作用時的冰排損傷分布圖如圖19和圖20所示。其中藍色部分表示未損傷單元，紅色表示已出現損傷的單元。從圖19和圖20中可以看出，與錐角α=30°和60°的錐體結構作用時，冰排的損傷模式并未有明顯的區別。冰排與錐體結構相互作用時，隨著錐體結構水平位移的增加，冰排在上表面和下表面同時出現失效，并且隨著錐體結構水平位移的增加，失效區域沿著徑向發展，形成徑向裂紋。當錐體結構的水平位移繼續增加時，之前形成的徑向裂紋沿著冰排的周向發展，形成周向裂紋，由于冰排的拉伸強度要明顯小于冰排的壓縮強度，所以冰排的下表面處先形成周向裂紋，且失效區域明顯大于上表面，最終使得冰排的上表面和下表面處與錐體結構相鄰的區域大面積失效。

圖19 錐體結構錐角α=30°時冰排損傷分布Fig.19 Failure process of cone structure with cone angleα=30°

冰排與不同錐角的錐體結構作用時，錐體結構受到的水平和垂直冰力與錐體結構的水平位移之間的關系如圖21和圖22所示。從圖21和圖22中可以看出，隨著錐體結構水平位移的逐漸增加，錐體結構受到的水平和垂直冰力也線性增加。當冰排中的應力狀態達到單元失效準則所定義的失效面時，冰排出現失效。錐體結構受到的水平和垂直冰力隨冰排水平位移的增加而下降。隨著錐體結構水平位移的繼續增加，冰排沿著錐體結構的斜面爬升，并且由于彈性基礎的存在，導致錐體結構受到的冰力再次隨錐體結構水平位移的增加而增加，使得錐體結構的受力再次繼續上升。

圖20 錐體結構錐角α=60°時冰排損傷分布Fig.20 Failure process of cone structure with cone angleα=60°

圖21 錐體結構受到的水平冰力與水平位移之間的關系Fig.21 Relationship between horizontal force and displacement of cone structure

圖22 錐體結構受到的垂直冰力與水平位移之間的關系Fig.22 Relationship between vertical force and displacement of cone structure

對于不同錐角的錐體結構，隨著錐體結構錐角的增加，冰排出現初始失效時對應的錐體結構受到的水平冰力也隨之增加。錐體結構錐角為30°時對應的冰排失效時結構受到的水平冰力為7.91×105N，而隨著錐體結構傾角的增加，當錐體結構錐角為60°時對應的水平冰力達到了2.84×106N，比錐角為30°時對應的水平冰力增加了359%。但冰排出現初始失效時對應的錐體結構受到的垂直冰力基本不隨斜坡結構傾角的變化而變化。錐體結構錐角為30°時對應的冰排失效時結構受到的垂直冰力為1.4×106N，而隨著錐體結構傾角的增加，當錐體結構錐角為60°時對應的垂直冰力為1.74×106N，相對錐角為30°時對應的垂直冰力沒有明顯的變化，只略微增加了20%。但冰排出現初始失效時對應的錐體結構受到的垂直冰力基本不隨錐體結構錐角的變化而變化。導致這一現象的原因與冰排和斜坡結構作用時類似，主要是由于冰排在與錐體結構作用時冰排因彎曲變形而出現彎折失效。而冰排的彎折失效由冰排受到的垂向分力引起，因為不同算例中，冰排的幾何尺寸和材料屬性并未改變，所以其發生彎折失效所對應的冰排受到的垂向分力也基本保持不變。但隨著錐體結構錐角的增加，冰排受到的水平分力也隨之增加，所以導致冰排出現初始失效時對應的錐體結構受到的水平冰力隨錐體結構錐角的增加而增加，但垂直冰力基本不隨錐體結構錐角的變化而變化。

為了驗證采用新型單元獲得的計算結果，將采用直接分析法獲得的冰排出現初始失效時的錐體結構受到的水平和垂直冰力與采用Ralston［3］提出的理論公式獲得的解析解進行了對比，如圖23和圖24所示。從圖中的結果對比可以看出，通過新型單元獲得的冰排出現初始失效時的錐體結構受到的水平冰力與解析解吻合良好，最大誤差在14%左右。但通過新型單元獲得的冰排出現初始失效時的錐體結構受到的垂直冰力與解析解的最大誤差在30%左右。引起誤差的原因主要有：

在所開發的新型單元中，冰被當作彈脆性材料來處理，而在Ralston［3］的理論公式中冰被當作理想彈塑性材料。由于彈塑性材料在達到屈服點后仍能承受一定的載荷，而彈脆性材料則不行，所以導致Ralston［3］給出的理論公式預測的冰力要比采用新型單元計算得到的冰力大。

圖23 直接計算法與Ralston［3］公式水平方向冰力結果對比Fig.23 Comparison of horizontal force between direct analysis method and the equation presented by Ralston［3］

圖24 直接計算法與Ralston［3］公式垂向冰力結果對比Fig.24 Comparison of vertical force between direct analysis method and the equation presented by Ralston［3］

在Ralston［3］的理論公式中，冰排與錐體結構之間的接觸被理想化為完全接觸，而在采用新型單元進行數值模擬時，從前述冰排與錐體結構相互作用時冰排的損傷分布圖中可以看出，冰排與錐體結構并非完全接觸。

冰排與錐體結構作用時，在Ralston［3］的理論公式中只考慮了冰排的彎曲，而忽略了錐體結構對冰排沿水平方向存在的面內力，而在采用新型單元進行數值模擬時，錐體結構存在對冰排的面內力，從而導致采用新型單元得到的冰力與Ralston［3］的理論公式給出的冰力存在一定的誤差。

3 結 論

本文基于ABAQUS自定義單元子程序，結合Reinicke和Remer提出的各向同性海冰失效準則和彈性基礎假定，開發了用于模擬海冰失效的8節點六面體單元。同時采用不同算例和相關理論公式對所開發的單元進行了驗證，得出以下結論：

1）進行不同算例的計算后發現，使用直接計算法得到的結果與規范公式得到的結果吻合良好，說明了直接計算法的可行性；

2）在計算規范公式所沒有的冰情和結構型式時直接計算法可以給出可信的結果，同時由于理論方法只能給出結構的最大受力，而直接計算法能夠直觀地給出冰損傷的過程及對應的結構受力變化過程，因而直接計算法可以為冰區結構設計提供一定的參考。

［1］BERCHA F G，DANYS J V.Prediction of ice forces on conical offshore structures［C］//Proceedings of the 3rd International Symposium on Ice Problems.Hanover，NH：Dartmouth College，1975：18-21.

［2］CROASDALE K R，CAMMAERT A B，METGE M.A method for the calculation of sheet ice loads on sloping structures［C］//Proceedings of the 12th IAHR Sympo?sium on Ice Problems.Trondheim，Norway：Norwe?gian Institute of Technology，1994.

［3］RALSTON T D.Plastic limit analysis of sheet ice loads on conical structures［M］//International Union of Theo?retical and Applied Mechanics Symposium on Physics and Mechanics of Ice.Berlin Heidelberg：Springer，1980：289-308.

［4］李鋒，岳前進.冰在斜面結構上的縱橫彎曲破壞分析［J］.水利學報，2000，31（9）：44-47.LI F，YUE Q J.Failure of ice sheet under simultane?ously longitudinal and transverse load on slope struc?tures［J］.Journal of Hydraulic Engineering，2000，31（9）：44-47（in Chinese）.

［5］Bureau Veritas.Ice characteristics and ice/structure in?teractions guidance note NI 565 DT R00 E［R］.［S.l.］：Bureau Veritas，2010.

［6］American Petroleum Institute.Recommended practice for planning，designing and constructing structures and pipelines for arctic conditions［S］.2nd ed.Wash?ington，DC：API Publications，1995.

［7］SAND B，HORRIGMOE G.Simulation of ice forces on sloping structures［C］//The Eighth International Off?shore and Polar Engineering Conference.Montreal，Canada：The International Society of Offshore and Po?lar Engineers，1998.

［8］DERRADJI-AOUAT A.Explicit FEA and constitutive modelling of damage and fracture in polycrystalline ice-simulations of ice loads on offshore structures［C］//18th International Conference on Port and Ocean Engi?neering under Arctic Conditions（POAC'05）.New York，United States：POAC，2005.

［9］黃焱，袁光奇，馬健鈞.海冰彎曲破壞下的破壞準則數值模擬方法研究［J］.數學的實踐與認識，2014，44（23）：115-126.HUANG Y，YUAN G Q，MA J J.Methodology re?search on the numerical simulation of the bending fail?ure of sea ice［J］.Mathematics in Practice and Theory，2014，44（23）：115-126（in Chinese）.

［10］宋安，孫金亮，史慶增，等.作用在引航導堤堤頭冰力的模型試驗研究［J］.天津大學學報，2006，39（5）：537-542.SONG A，SUN J L，SHI Q Z，et al.Model test for ice force on the bank-head of the lead-navigating bank［J］.Journal of Tianjin University，2006，39（5）：537-542（in Chinese）.

［11］REINICKE K M，REMER R.A procedure for the de?termination of ice forces-Illustrated for polycrystal?line ice［C］//Proceedings of the 4th IAHR Symposium on Ice Problems.Lule?，Sweden：［s.n.］，1978.

［12］ZIENKIEWICZ O C，TAYLOR R L.有限元方法：第1卷：基本原理［M］.曾攀，譯.5版.北京：清華大學出版社，2008.

［13］孔娟，張勇.彈性地基板單元在ABAQUS中的開發與 應 用［J］.燕 山 大 學 學 報 ，2010，34（4）：370-376.KONG J，ZHANG Y.Secondary development and ap?plication ofelastic foundation plate elementin ABAQUS［J］.Journal of Yanshan University，2010，34（4）：370-376（in Chinese）.

［14］CROASDALE K R.Forces on fixed rigid structures.Working group on ice forces on structures-a state-of-the-art-report［R］.Hanover，New Hamp?shire：CRREL，1980.

Development of sea ice element and direct analysis method of predicting ice loads based on ABAQUS

GONG Yufeng1，ZHANG Zhengyi1，2，3，LIU Jingxi1，2，3，DONG Wen1，XIE De1，2，3
1 School of Naval Architecture and Ocean Engineering，Huazhong University of Science and Technology，Wuhan 430074，China
2 Hubei Key Laboratory of Naval Architecture&Ocean Engineering Hydrodynamics，Wuhan 430074，China
3 Collaborative Innovation Center for Advanced Ship and Deep-Sea Exploration，Shanghai 200240，China

With the global temperature continuing to rise,the Arctic ice is melting.The development of the Arctic route and exploitation of Arctic resources are drawing the wide attention of international society.Many countries are starting a new round of the construction of polar ships.In recent years,domestic and foreign researchers have begun to focus on the numerical study of ice loads on structure.However,the current commercial finite element software lacks a specialized element used for simulating the mechanical properties of sea ice.This lack increases the difficulty of the numerical study of ice loads on structure.In this paper,based on Reinicke and Remer's elliptic failure criteria and Winkler's elastic foundation,an element has been developed to simulate the bending failure of isotropic granular ice.The ice loads of typical marine structures,such as slope structure and conical structure,are calculated through the direct analysis method with the developed element.A comparison between the results of the direct analysis method and analytical method show good agreement.In brief,the proposed direct analysis method may provide tools and certain references for structural design in an ice environment.

ice load；direct analysis method；ABAQUS；hourglass control；user-defined element

U661.4

：ADOI：10.3969/j.issn.1673-3185.2017.03.011

http://kns.cnki.net/kcms/detail/42.1755.TJ.20170512.1251.020.html期刊網址：www.ship-research.com

龔榆峰，張正藝，劉敬喜，等.基于ABAQUS的海冰單元開發及冰載荷直接計算法［J］.中國艦船研究，2017，12（3）：75-85.

GONG Y F，ZHANG Z Y，LIU J X，et al.Development of sea ice element and direct analysis method of predicting ice loads based on ABAQUS［J］.Chinese Journal of Ship Research，2017，12（3）：75-85.

2016-10-12< class="emphasis_bold">網絡出版時間

時間：2017-5-12 12:51

國家自然科學基金資助項目（51579110）；國家高技術研究發展計劃資助項目（2012AA112601）；華中科技大學自主創新研究基金資助項目（2015TS004）

龔榆峰，男，1988年生，博士生。研究方向：船體結構。E-mail：D201277373@hust.edu.cn

劉敬喜（通信作者），男，1975年生，博士，副教授。研究方向：船體結構。E-mail：Liu_Jing_Xi@hust.edu.cn