基于ABAQUS的海冰單元開(kāi)發(fā)及冰載荷直接計(jì)算法

龔榆峰，張正藝，2，3，劉敬喜，2，3，董問(wèn)，解德，2，3

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

2船舶和海洋水動(dòng)力湖北省重點(diǎn)實(shí)驗(yàn)室，湖北武漢430074

3高新船舶與深海開(kāi)發(fā)裝備協(xié)同創(chuàng)新中心，上海200240

基于ABAQUS的海冰單元開(kāi)發(fā)及冰載荷直接計(jì)算法

龔榆峰1，張正藝1，2，3，劉敬喜1，2，3，董問(wèn)1，解德1，2，3

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

2船舶和海洋水動(dòng)力湖北省重點(diǎn)實(shí)驗(yàn)室，湖北武漢430074

3高新船舶與深海開(kāi)發(fā)裝備協(xié)同創(chuàng)新中心，上海200240

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

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

0 引 言

隨著全球氣溫的升高，北極地區(qū)的海冰逐漸融化，人類對(duì)極地船舶與海洋結(jié)構(gòu)物的需求也隨之增加。在進(jìn)行極地船舶與海洋結(jié)構(gòu)物的設(shè)計(jì)時(shí)，如何合理、有效確定結(jié)構(gòu)物與海冰相互作用過(guò)程中受到的冰力，是一個(gè)長(zhǎng)期被國(guó)內(nèi)外相關(guān)研究人員關(guān)注的問(wèn)題。然而海冰與結(jié)構(gòu)物的相互作用是一個(gè)相當(dāng)復(fù)雜的過(guò)程，在目前的研究中，研究人員通常將海冰與不同結(jié)構(gòu)物相互作用過(guò)程中海冰的破壞模式歸納為壓碎、彎折和劈裂等幾類。

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

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

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

本文將基于通用商業(yè)有限元軟件ABAQUS的自定義單元子程序創(chuàng)建集成了Reinicke和Remer［11］提出的各向同性海冰失效準(zhǔn)則和彈性基礎(chǔ)的8節(jié)點(diǎn)六面體單元用于模擬海冰。同時(shí)為了對(duì)所開(kāi)發(fā)的單元進(jìn)行驗(yàn)證，將采用開(kāi)發(fā)的單元計(jì)算冰排與斜坡結(jié)構(gòu)相互作用過(guò)程中的冰力大小，并將所獲得的計(jì)算結(jié)果與BV和API給出的相關(guān)規(guī)范公式進(jìn)行對(duì)比。

1 單元及失效準(zhǔn)則介紹

1.1 所開(kāi)發(fā)單元基本原理

8節(jié)點(diǎn)六面體縮減積分單元示意圖如圖1所示。

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

根據(jù)有限元理論［12］，單元內(nèi)任意一點(diǎn)的位移u可以表示為

式中：u，v，w為沿坐標(biāo)軸x，y，z軸的位移分量；x，y，z為總體坐標(biāo)系的3個(gè)坐標(biāo)軸。單元內(nèi)任意一點(diǎn)的位移u可以由單元的形函數(shù)和單元各個(gè)節(jié)點(diǎn)的位移來(lái)表示：

式中：Na為基于節(jié)點(diǎn)位置事先給定的函數(shù)，又稱單元的形函數(shù)；ua為單元內(nèi)各個(gè)節(jié)點(diǎn)的位移。

單元任意一點(diǎn)的應(yīng)變可用引起應(yīng)變的3個(gè)分量來(lái)表示：

式中：S為單元應(yīng)力矩陣；Bi為單元應(yīng)變矩陣。

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

根據(jù)應(yīng)力—應(yīng)變關(guān)系，采用矩陣D可以將單元中的6個(gè)應(yīng)力分量和6個(gè)應(yīng)變分量聯(lián)系起來(lái)，單元中應(yīng)力σ可以表示為

最后，由虛功原理可以得到單元?jiǎng)偠染仃嘖，

式中，J為雅克比行列式。采用高斯積分時(shí)，式（6）可以簡(jiǎn)化為

式中，H為權(quán)函數(shù)。由于本文所開(kāi)發(fā)的單元為完全縮減積分單元，只有一個(gè)積分點(diǎn)，即式（7）中i=1。對(duì)應(yīng)的積分點(diǎn)在局部坐標(biāo)系下的坐標(biāo)為（0，0，0），即式（7）中ξi=0，ηi=0，ζi=0，對(duì)應(yīng)于積分點(diǎn)數(shù)量的權(quán)函數(shù)H=8.0。

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

由節(jié)點(diǎn)力平衡方程可以得到彈性基礎(chǔ)反力和垂向位移之間的關(guān)系：

式中，K0為考慮彈性基礎(chǔ)而引入的單元?jiǎng)偠雀郊泳仃嚒?/p>

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

式中，Ke為不考慮彈性基礎(chǔ)時(shí)的單元?jiǎng)偠染仃嚒?/p>

1.2 Reinicke和Remer各向同性海冰失效準(zhǔn)則

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

式中：k1，k2和k3為材料常數(shù)，通過(guò)試驗(yàn)確定；I1為應(yīng)力張量第1不變量；J2為應(yīng)力偏量第2不變量。當(dāng)f(I1，J2，k1，k2，k3)=0 時(shí)，材料失效。失效后的材料彈性模量和泊松比分別變?yōu)镋c=60 MPa，vc=0.33。本文中使用的材料常數(shù)的取值參考文獻(xiàn)［9］中由模型試驗(yàn)擬合得到的結(jié)果，k1=0.674（MPa）-2，k2=0.853（MPa）-1，k3=0.014（MPa）-2。

1.3 沙漏控制

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

圖2 單元計(jì)算流程圖Fig.2 Flow chart of element computation

2 算例驗(yàn)證

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

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

表1 計(jì)算中使用材料參數(shù)Table 1 The material constants used in the simulation

2.1 單元基本特性驗(yàn)證

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

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

使用新型單元得到的計(jì)算結(jié)果與Reinicke和Remer屈服準(zhǔn)則解析解的對(duì)比如圖4所示。將新型單元的計(jì)算結(jié)果與解析解的對(duì)比可以發(fā)現(xiàn)，該新型單元不僅能夠描述海冰的拉伸與壓縮強(qiáng)度不一致這一基本力學(xué)特性，而且計(jì)算得到的拉伸強(qiáng)度和壓縮強(qiáng)度與Reinicke和Remer屈服準(zhǔn)則解析解吻合很好，驗(yàn)證了本文提出的新型單元能夠很好地實(shí)現(xiàn)Reinicke和Remer提出的海冰本構(gòu)關(guān)系，并描述了海冰的基本力學(xué)特性。

圖4 應(yīng)力—應(yīng)變曲線計(jì)算結(jié)果對(duì)比Fig.4 Comparison of stress-strain curve

2.2 半無(wú)限長(zhǎng)冰梁

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

圖5 半無(wú)限長(zhǎng)冰梁示意圖Fig.5 Schematic of semi-infinite ice beam

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

圖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

一端受到垂向載荷的不同厚度半無(wú)限長(zhǎng)冰梁的冰力和端部位移之間的關(guān)系曲線如圖8所示。從圖8中可以看出，隨著垂向位移的逐漸增加，端部受到的反力也線性增加，當(dāng)冰梁中的應(yīng)力狀態(tài)達(dá)到單元失效準(zhǔn)則所定義的失效面時(shí)，冰梁出現(xiàn)失效。當(dāng)單元出現(xiàn)失效時(shí)，直立結(jié)構(gòu)受到的冰力隨位移的增加出現(xiàn)明顯的下降。但由于彈性基礎(chǔ)的存在，以及首先出現(xiàn)失效的只是冰梁上表面的部分單元，冰梁并未全部失效，故仍能繼續(xù)承受載荷。所以加載端部受到的反力隨位移的增加，出現(xiàn)下降后又會(huì)再次繼續(xù)上升。

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

為了驗(yàn)證采用新型單元獲得的計(jì)算結(jié)果，將采用直接分析法獲得的冰梁出現(xiàn)初始失效時(shí)的失效位置、端部承受的反力與BV［5］給出的理論公式獲得的解析解進(jìn)行了對(duì)比，如圖9和圖10所示。從圖中的結(jié)果對(duì)比可以看出，通過(guò)新型單元獲得的失效位置和對(duì)應(yīng)的端部反力與解析解吻合良好，誤差僅為5%左右。說(shuō)明所開(kāi)發(fā)的新型單元能夠準(zhǔn)確模擬一端受垂向載荷的半無(wú)限長(zhǎng)冰梁的彎折失效。

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

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

2.3 冰排與斜坡結(jié)構(gòu)的相互作用

為了驗(yàn)證所開(kāi)發(fā)的單元，計(jì)算了斜坡結(jié)構(gòu)與冰排作用時(shí)受到的冰力，并將所得計(jì)算結(jié)果與Croasdale［14］提出的理論公式進(jìn)行了對(duì)比。模型示意圖如圖11所示。

圖11 斜坡結(jié)構(gòu)示意圖Fig.11 Schematic of sloping structure

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

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

圖12 斜坡結(jié)構(gòu)傾角α=30°時(shí)冰排損傷分布Fig.12 Failure process of sloping structure with slope angleα=30°

圖13 斜坡結(jié)構(gòu)傾角α=60°時(shí)冰排損傷分布Fig.13 Failure process of sloping structure with slope angleα=60°

冰排與不同傾角的斜坡結(jié)構(gòu)作用時(shí)，斜坡結(jié)構(gòu)受到的水平和垂直冰力與冰排的水平位移之間的關(guān)系如圖14和圖15所示。從圖14和圖15中可以看出，隨著冰排水平位移的逐漸增加，斜坡結(jié)構(gòu)受到的水平和垂直冰力也線性增加。當(dāng)冰排中的應(yīng)力狀態(tài)達(dá)到單元失效準(zhǔn)則所定義的失效面時(shí)，冰排出現(xiàn)失效。斜坡結(jié)構(gòu)受到的水平和垂直冰力隨冰排水平位移的增加而下降。隨著冰排水平位移的繼續(xù)增加，冰排沿著斜坡結(jié)構(gòu)的斜面爬升，并且由于彈性基礎(chǔ)的存在，導(dǎo)致斜坡結(jié)構(gòu)受到的冰力再次隨位移的增加而增加，使得斜坡結(jié)構(gòu)的受力再次繼續(xù)上升。

圖14 斜坡結(jié)構(gòu)受到的水平冰力與位移之間的關(guān)系Fig.14 Relationship between horizontal force and displacement of sloping structure

圖15 斜坡結(jié)構(gòu)受到的垂直冰力與水平位移之間的關(guān)系Fig.15 Relationship between vertical force and displacement of sloping structure

對(duì)于不同傾角的斜坡結(jié)構(gòu)，隨著斜坡結(jié)構(gòu)傾角的增加，冰排出現(xiàn)初始失效時(shí)對(duì)應(yīng)的斜坡結(jié)構(gòu)受到的水平冰力也隨之增加。斜坡結(jié)構(gòu)傾角為30°時(shí)對(duì)應(yīng)的冰排失效時(shí)結(jié)構(gòu)受到的水平冰力為5.75×104N，而隨著斜坡結(jié)構(gòu)傾角的增加，當(dāng)斜坡結(jié)構(gòu)傾角為60°時(shí)對(duì)應(yīng)的水平冰力達(dá)到2.02×105N，比傾角為30°時(shí)對(duì)應(yīng)的水平冰力增加了351%。但冰排出現(xiàn)初始失效時(shí)對(duì)應(yīng)的斜坡結(jié)構(gòu)受到的垂直冰力基本不隨斜坡結(jié)構(gòu)傾角的變化而變化。斜坡結(jié)構(gòu)傾角為30°時(shí)對(duì)應(yīng)的冰排失效時(shí)結(jié)構(gòu)受到的垂直冰力為9.96×104N，而隨著斜坡結(jié)構(gòu)傾角的增加，當(dāng)斜坡結(jié)構(gòu)傾角為60°時(shí)對(duì)應(yīng)的垂直冰力為1.16×105N，相對(duì)傾角為30°時(shí)對(duì)應(yīng)的垂直冰力沒(méi)有明顯的變化，只略微增加了16%。這主要是由于冰排在與斜坡結(jié)構(gòu)作用時(shí)因?yàn)楸虐l(fā)生彎曲變形而出現(xiàn)了彎折失效。而冰排的彎折失效主要由冰排受到的垂向分力引起，因?yàn)椴煌憷斜诺膸缀纬叽绾筒牧蠈傩圆⑽锤淖?，所以其發(fā)生彎折失效所對(duì)應(yīng)冰排受到的垂向分力基本保持不變。但隨著斜坡角度的增加，冰排受到的水平分力隨之增加，所以導(dǎo)致冰排出現(xiàn)初始失效時(shí)對(duì)應(yīng)的斜坡結(jié)構(gòu)受到的水平冰力隨斜坡結(jié)構(gòu)傾角的增加而增加，而垂直冰力基本不隨斜坡結(jié)構(gòu)傾角的變化而變化。

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

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

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

2.4 冰排與圓錐結(jié)構(gòu)的相互作用

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

圖18 圓錐結(jié)構(gòu)示意圖Fig.18 Schematic of cone structure

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

圖19 錐體結(jié)構(gòu)錐角α=30°時(shí)冰排損傷分布Fig.19 Failure process of cone structure with cone angleα=30°

冰排與不同錐角的錐體結(jié)構(gòu)作用時(shí)，錐體結(jié)構(gòu)受到的水平和垂直冰力與錐體結(jié)構(gòu)的水平位移之間的關(guān)系如圖21和圖22所示。從圖21和圖22中可以看出，隨著錐體結(jié)構(gòu)水平位移的逐漸增加，錐體結(jié)構(gòu)受到的水平和垂直冰力也線性增加。當(dāng)冰排中的應(yīng)力狀態(tài)達(dá)到單元失效準(zhǔn)則所定義的失效面時(shí)，冰排出現(xiàn)失效。錐體結(jié)構(gòu)受到的水平和垂直冰力隨冰排水平位移的增加而下降。隨著錐體結(jié)構(gòu)水平位移的繼續(xù)增加，冰排沿著錐體結(jié)構(gòu)的斜面爬升，并且由于彈性基礎(chǔ)的存在，導(dǎo)致錐體結(jié)構(gòu)受到的冰力再次隨錐體結(jié)構(gòu)水平位移的增加而增加，使得錐體結(jié)構(gòu)的受力再次繼續(xù)上升。

圖20 錐體結(jié)構(gòu)錐角α=60°時(shí)冰排損傷分布Fig.20 Failure process of cone structure with cone angleα=60°

圖21 錐體結(jié)構(gòu)受到的水平冰力與水平位移之間的關(guān)系Fig.21 Relationship between horizontal force and displacement of cone structure

圖22 錐體結(jié)構(gòu)受到的垂直冰力與水平位移之間的關(guān)系Fig.22 Relationship between vertical force and displacement of cone structure

對(duì)于不同錐角的錐體結(jié)構(gòu)，隨著錐體結(jié)構(gòu)錐角的增加，冰排出現(xiàn)初始失效時(shí)對(duì)應(yīng)的錐體結(jié)構(gòu)受到的水平冰力也隨之增加。錐體結(jié)構(gòu)錐角為30°時(shí)對(duì)應(yīng)的冰排失效時(shí)結(jié)構(gòu)受到的水平冰力為7.91×105N，而隨著錐體結(jié)構(gòu)傾角的增加，當(dāng)錐體結(jié)構(gòu)錐角為60°時(shí)對(duì)應(yīng)的水平冰力達(dá)到了2.84×106N，比錐角為30°時(shí)對(duì)應(yīng)的水平冰力增加了359%。但冰排出現(xiàn)初始失效時(shí)對(duì)應(yīng)的錐體結(jié)構(gòu)受到的垂直冰力基本不隨斜坡結(jié)構(gòu)傾角的變化而變化。錐體結(jié)構(gòu)錐角為30°時(shí)對(duì)應(yīng)的冰排失效時(shí)結(jié)構(gòu)受到的垂直冰力為1.4×106N，而隨著錐體結(jié)構(gòu)傾角的增加，當(dāng)錐體結(jié)構(gòu)錐角為60°時(shí)對(duì)應(yīng)的垂直冰力為1.74×106N，相對(duì)錐角為30°時(shí)對(duì)應(yīng)的垂直冰力沒(méi)有明顯的變化，只略微增加了20%。但冰排出現(xiàn)初始失效時(shí)對(duì)應(yīng)的錐體結(jié)構(gòu)受到的垂直冰力基本不隨錐體結(jié)構(gòu)錐角的變化而變化。導(dǎo)致這一現(xiàn)象的原因與冰排和斜坡結(jié)構(gòu)作用時(shí)類似，主要是由于冰排在與錐體結(jié)構(gòu)作用時(shí)冰排因彎曲變形而出現(xiàn)彎折失效。而冰排的彎折失效由冰排受到的垂向分力引起，因?yàn)椴煌憷?，冰排的幾何尺寸和材料屬性并未改變，所以其發(fā)生彎折失效所對(duì)應(yīng)的冰排受到的垂向分力也基本保持不變。但隨著錐體結(jié)構(gòu)錐角的增加，冰排受到的水平分力也隨之增加，所以導(dǎo)致冰排出現(xiàn)初始失效時(shí)對(duì)應(yīng)的錐體結(jié)構(gòu)受到的水平冰力隨錐體結(jié)構(gòu)錐角的增加而增加，但垂直冰力基本不隨錐體結(jié)構(gòu)錐角的變化而變化。

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

在所開(kāi)發(fā)的新型單元中，冰被當(dāng)作彈脆性材料來(lái)處理，而在Ralston［3］的理論公式中冰被當(dāng)作理想彈塑性材料。由于彈塑性材料在達(dá)到屈服點(diǎn)后仍能承受一定的載荷，而彈脆性材料則不行，所以導(dǎo)致Ralston［3］給出的理論公式預(yù)測(cè)的冰力要比采用新型單元計(jì)算得到的冰力大。

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

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

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

冰排與錐體結(jié)構(gòu)作用時(shí)，在Ralston［3］的理論公式中只考慮了冰排的彎曲，而忽略了錐體結(jié)構(gòu)對(duì)冰排沿水平方向存在的面內(nèi)力，而在采用新型單元進(jìn)行數(shù)值模擬時(shí)，錐體結(jié)構(gòu)存在對(duì)冰排的面內(nèi)力，從而導(dǎo)致采用新型單元得到的冰力與Ralston［3］的理論公式給出的冰力存在一定的誤差。

3 結(jié) 論

本文基于ABAQUS自定義單元子程序，結(jié)合Reinicke和Remer提出的各向同性海冰失效準(zhǔn)則和彈性基礎(chǔ)假定，開(kāi)發(fā)了用于模擬海冰失效的8節(jié)點(diǎn)六面體單元。同時(shí)采用不同算例和相關(guān)理論公式對(duì)所開(kāi)發(fā)的單元進(jìn)行了驗(yàn)證，得出以下結(jié)論：

1）進(jìn)行不同算例的計(jì)算后發(fā)現(xiàn)，使用直接計(jì)算法得到的結(jié)果與規(guī)范公式得到的結(jié)果吻合良好，說(shuō)明了直接計(jì)算法的可行性；

2）在計(jì)算規(guī)范公式所沒(méi)有的冰情和結(jié)構(gòu)型式時(shí)直接計(jì)算法可以給出可信的結(jié)果，同時(shí)由于理論方法只能給出結(jié)構(gòu)的最大受力，而直接計(jì)算法能夠直觀地給出冰損傷的過(guò)程及對(duì)應(yīng)的結(jié)構(gòu)受力變化過(guò)程，因而直接計(jì)算法可以為冰區(qū)結(jié)構(gòu)設(shè)計(jì)提供一定的參考。

［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ìn).冰在斜面結(jié)構(gòu)上的縱橫彎曲破壞分析［J］.水利學(xué)報(bào)，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］黃焱，袁光奇，馬健鈞.海冰彎曲破壞下的破壞準(zhǔn)則數(shù)值模擬方法研究［J］.數(shù)學(xué)的實(shí)踐與認(rèn)識(shí)，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］宋安，孫金亮，史慶增，等.作用在引航導(dǎo)堤堤頭冰力的模型試驗(yàn)研究［J］.天津大學(xué)學(xué)報(bào)，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版.北京：清華大學(xué)出版社，2008.

［13］孔娟，張勇.彈性地基板單元在ABAQUS中的開(kāi)發(fā)與 應(yīng) 用［J］.燕 山 大 學(xué) 學(xué) 報(bào) ，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期刊網(wǎng)址：www.ship-research.com

龔榆峰，張正藝，劉敬喜，等.基于ABAQUS的海冰單元開(kāi)發(fā)及冰載荷直接計(jì)算法［J］.中國(guó)艦船研究，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">網(wǎng)絡(luò)出版時(shí)間

時(shí)間：2017-5-12 12:51

國(guó)家自然科學(xué)基金資助項(xiàng)目（51579110）；國(guó)家高技術(shù)研究發(fā)展計(jì)劃資助項(xiàng)目（2012AA112601）；華中科技大學(xué)自主創(chuàng)新研究基金資助項(xiàng)目（2015TS004）

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

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