船-冰作用中海冰典型破壞模式的離散元數值模擬

雷建奇, 狄少丞*, 于加一, 王迎暉

(1.哈爾濱工程大學，哈爾濱 150001； 2.中國船舶重工集團有限公司第七○二研究所，無錫 214082)

1 引 言

近年來，極地資源開發以及極地航道開辟越來越受到人們的關注與重視[1]，冰區船舶作為進行這兩項工作的重要載體，如何保證其在冰區航行時的安全性顯得尤為重要。冰區船舶在平整冰區航行時，船體與冰排的相互作用給船體帶來的冰載荷對船舶的安全性與適航性有重要影響，船體所受冰載荷主要來源于冰排破碎[2]。對于船-冰相互作用，國內外學者們通過各種手段進行了廣泛的研究，但到目前為止人們對于這一問題的理解還十分有限[3]。海冰破碎時表現出離散的特性，同時不同尺度下海冰分布也十分離散[4,5]，因此使用離散元方法模擬海冰與結構物的作用具有一定優勢[6,7]，能較好地模擬船-冰作用中海冰從完整狀態到破碎狀態的過程[8]。

Shen等[9]最先在海冰動力學的研究中引入基于顆粒理論的離散元方法，在有限的假設下推導了理想二維破碎冰流場的碰撞流變學，并開發了與之對應的蒙特卡洛模擬方法。Lau等[10-11]使用三維離散元方法模擬了Canadian號破冰船與冰排的相互作用，之后進一步采用三維塊體離散元模擬了破冰船在冰區中的直線航行和回轉運動，分析了破冰船在不同航行狀態下受到的冰阻力。Karulin等[12]以一組具有粘彈性體相互作用的漂浮圓盤代表浮冰，計算不同漂浮速度下系泊船舶與浮冰的相互作用，同時將數值計算結果與冰池模型試驗進行對比，發現所建立的數值模型可用于浮冰中系泊油輪的理論分析。

國內的學者們在船-冰相互作用的離散元模擬方面也做了大量工作。蔡柯等[13]采用離散單元法以球形顆粒粘結成平整冰，同時考慮海流對海冰單元的浮力和拖曳力，模擬船舶在平整冰區中航行，分析船舶結構與海冰的相互作用，探討冰載荷分布規律。劉璐等[14]建立了基于擴展多面體的離散元方法，并以雪龍號為對象，研究了單船破冰條件下的冰荷載，將離散元計算結果與常用經驗公式進行對比。剛旭皓等[15]以具有粘結-破碎特性的球體單元、三維拓展圓盤和拓展多面體單元構造海冰的離散元模型，模擬船體在平整冰中的直航過程，觀察冰排破碎特性。

船-冰相互作用過程中，擠壓與彎曲破壞是冰排的兩種典型破壞模式[16，17]，不同破壞模式給船舶帶來的外載荷有較大差異。冰排的擠壓破壞與彎曲破壞受多種因素的影響，如冰厚、航速和船-冰作用部位等。本文通過離散元數值手段，建立海冰擠壓與彎曲破壞的數值模型，對這兩種典型破壞模式及影響因素進行分析。

2 海冰離散元模型

2.1 模型建立

采用具有粘結-破碎特性的離散元法對海冰的破壞特性進行分析[18]，通常采用單軸壓縮試驗和三點彎曲試驗對模型的細觀參數進行標定。首先，建立海冰單軸壓縮與三點彎曲數值模型，其中單軸壓縮數值試樣尺寸為70 mm×70 mm×175 mm，三點彎曲數值試樣尺寸為75 mm×75 mm×700 mm[19]，試樣中顆粒粒徑服從正態分布，建立的數值試樣如圖1所示，主要計算參數列入表1。

圖1 海冰壓縮試樣與彎曲試樣的離散元模型

表1 數值試驗的主要計算參數

海冰壓縮與彎曲試驗均采用恒定速率加載的方式，數值試驗中也采用同樣的加載方式。對于同尺寸的海冰試樣，文獻[20,21]進行了單軸壓縮與三點彎曲試驗研究，可得海冰壓縮強度平均取值為5.3 MPa，壓縮強度與彎曲強度的比值一般在2.5～3。

離散元模擬的海冰單軸壓縮與三點彎曲試驗中，海冰破碎狀態如圖2所示，計算得到的海冰壓縮與彎曲應力-應變和應力-撓度曲線如圖3所示，曲線中將直線段AB的斜率作為海冰試樣的彈性模量E，將壓縮應力的最大值作為海冰的壓縮強度σc，將彎曲應力的最大值作為海冰的彎曲強度σf。

圖2 離散單元法模擬的海冰單軸壓縮以及三點彎曲破裂情況

圖3 離散單元法模擬的海冰單軸壓縮以及三點彎曲試驗結果

2.2 細觀剛度值的標定

采用具有黏結-破碎特性的離散元法模擬海冰與海洋結構物相互作用時，需要對模型的細觀法向剛度和切向剛度兩個剛度參數以及法向黏結強度、切向黏結強度和摩擦系數三個強度參數進行標定。其中，細觀剛度系數對試樣的宏觀變形模量產生影響，細觀強度參數對試樣的宏觀強度產生影響[22]。

模型需要校核的細觀剛度參數有法向剛度系數(Kn)和切向法向剛度比(Ks/Kn)，校核時初始取值為1.60 MN/m，1.83 MN/m，2.00 MN/m，2.17 MN/m 和2.40 MN/m，以及0.10，0.33，0.50，0.67和0.90，細觀強度參數暫時定為法向粘結強度2.0 MPa、切向粘結強度2.0 MPa和摩擦角 12.5°。繪制出細觀剛度系數對彈性模量的影響規律，如圖4所示。

圖4 法向剛度和剛度比對彈性模量影響關系

選取海冰彈性模量2 GPa作為細觀剛度值標定的依據，從圖4可選定細觀剛度值為法向剛度Kn=1.8 MN/m，切向法向剛度比Ks/Kn=0.5。

2.3 細觀強度值的標定

模型中需要校核的細觀強度參數有法向粘結強度(Ten)、切向粘結強度(Coh)以及摩擦角(Fa)。其中，法向與切向粘結強度的校核取值分別為 1 MPa,2 MPa,2.5 MPa,3 MPa,3.5 MPa,4 MPa和5 MPa，摩擦角取值范圍為0°，12°，22°和25°。經計算得到壓縮強度等值線和壓縮強度與彎曲強度比值等值線，如圖5～圖8所示。

圖5 摩擦角0°時的等值線

圖6 摩擦角12°時的等值線

圖7 摩擦角22°時的等值線

圖8 摩擦角25°時的等值線

分析以上法向黏結強度與切向黏結強度對宏觀壓縮強度與壓縮彎曲強度比值的影響規律，結合試驗中壓縮強度和壓縮彎曲強度比值的參考值，最終選定的細觀強度參數為法向黏結強度3 MPa、切向黏結強度4 MPa以及摩擦角12°。

3 擠壓破壞過程離散元模擬

3.1 離散元模型建立

船冰作用過程中，冰排主要以船艏部位的彎曲破壞和船體舯部的擠壓破壞為主，同時伴隨著擠壓-彎曲共同作用及破壞模式。將船體與冰排作用部位簡化成具有不同傾角的平板，采用離散元數值方法，對具有不同傾角、不同作用方位及不同船-冰相對速度下的冰排破壞模式進行研究。分別構建了1.5 m×1.5 m×0.5 m,3 m×3 m×1 m和 6 m×6 m×2 m三組不同尺寸的冰排離散元模型，如圖9所示。經研究表明，不同尺寸的試樣中，在試樣尺寸與顆粒粒徑比值一致的情況下，試樣表現出的宏觀力學行為是一致的。因此，三組冰排在厚度方向上具有與單軸壓縮和三點彎曲試樣同樣的顆粒數量，以保證經標定后的細觀參數仍適用于較大尺寸冰排破壞過程的模擬。

圖9 長方體冰排離散元模型

將船-冰擠壓作用方式進一步劃分為側棱撞擊、頂角撞擊、下壓以及側舷刮擦四種作用方式，如圖10所示。冰和平板的相對速度設定為1節、2節和4節。將速度單位換算成m/s，A組中平板沿X軸和Y軸負方向的速度均為1.5 m/s；B組沿X軸、Y軸和Z軸負方向的速度均為1.2 m/s；C組沿Y軸和Z軸負方向的速度均為1.5 m/s；D組沿X軸負方向的速度為2.0 m/s，沿Y軸負方向的速度較小為0.1 m/s。

圖10 船-冰擠壓作用方式

3.2 作用方位對擠壓破壞過程的影響

選取冰厚為0.5 m，相對速度2節的工況，計算分析A，B，C，D四種作用方位下的冰排破壞過程，計算得到的冰力時程曲線如圖11所示，冰排破壞結果如圖12所示。

圖11 4種作用方位下的冰力時程曲線

圖12 4種加載方位下的冰排破壞方式

平板受力可反應出不同作用方位下船體的受力狀況。從冰力時程曲線可以看出，當平板以側棱撞擊(A工況)和側舷刮擦(D工況)方位與冰排作用時的冰力較大，均達到了250 kN，當以頂角撞擊(B工況)和下壓(C工況)方位作用時的冰力相對較小，均約為100 kN。

從冰排破壞方式來看，側棱撞擊(A工況)和側舷刮擦(D工況)方位下,冰排的破壞方式以水平面內的擠壓破壞為主;頂角撞擊(B工況)和下壓(C工況)方位下，冰排的破壞方式除擠壓外還存在彎曲破壞的成分，這就導致了A工況和D工況下平板的冰力較大。

3.3 冰厚對擠壓破壞的影響

分析相對速度2節四種作用方位、三種冰厚為0.5 m，1 m和2 m下冰排擠壓破壞過程，計算得到的冰力時程曲線如圖13和圖14所示。可以看出，不同作用方位下的平板受力均呈現出了隨冰厚增加而增大的趨勢，均體現為冰厚增加使得冰排與平板的接觸面積增大，進而使得冰力增加。就增幅而言，頂角撞擊(B工況)、下壓(C工況)及側舷刮擦(D工況)方位冰力增幅大于側棱撞擊(A工況)方位；對于B工況和C工況，隨著冰厚的增加，冰排破壞方式中彎曲破壞的比例逐漸降低，擠壓破壞的成分升高，最終導致冰力的增大。

圖13 A，B和C三種作用方位下，不同厚度冰排冰力比較

圖14 D組側舷刮擦作用方式下，不同厚度冰排冰力對比

3.4 相對速度對擠壓破壞的影響

分析三種相對速度1節、2節和4節、四種作用方位、冰厚1 m工況下的冰排擠壓破壞過程。由于平板速度不同，采用冰力-位移曲線來保證平板在冰排中的穿透量一致，計算得到的冰力-位移曲線如圖15所示。可以看出，側棱撞擊(A工況)和頂角撞擊(B工況)方位下，冰力的增大與速度呈正比關系；下壓(C工況)方位下，1節與2節速度對冰力的影響不明顯，4節速度冰力約為1節和2節速度冰力的4倍；側舷刮擦(D工況)方位下，速度改變對冰力的變化無影響。

圖15 三種加載速度下，不同加載方式的冰力對比

4 彎曲破壞過程離散元模擬

4.1 離散元模型建立

參照IACS極地規則的冰楔尺寸，建立了上底為0.5 m，厚度為0.5 m和1 m，夾角為30°，45°和60°的等腰梯形冰排模型，如圖16所示，通過平板勻速向下運動對冰排進行加載，計算分析不同楔形冰排在豎向載荷下的破壞模式，與擠壓破壞不同的是，楔形冰排還受到下方水體的浮力(冰密度為920 kg/m3，水密度為1035 kg/m3)。

圖16 冰楔與離散元模型

剛性固定梯形冰排底部三邊，以模擬實際工況中平整冰彎曲斷裂時的三側受約束情況。分別賦予平板0.15 m/s，0.3 m/s和0.6 m/s的加載速度，研究不同夾角、不同冰厚以及不同加載速度對彎曲破壞過程的影響。計算得到的冰楔夾角60°、加載速度0.15 m/s和冰厚0.5 m的冰排破壞后的情況如圖17所示，斷裂位置位于冰排的底部，并呈圓弧形。

圖17 冰排彎曲斷裂部位(加載速度0.15 m/s，冰厚0.5 m)

4.2 冰楔夾角對彎曲破壞冰力的影響

分析冰厚0.5 m和加載速度0.15 m/s下，不同夾角冰排的彎曲破壞過程，得到冰力時程曲線如 圖18 所示。可以看出，冰力曲線具有典型的彎曲破壞特征，隨著冰排夾角的增大，冰力變化不顯著。

圖18 不同夾角冰排破壞過程中的冰力時程曲線

4.3 冰排厚度對彎曲破壞冰力的影響

分析加載速度為0.15 m/s，三種冰楔夾角以及兩種冰厚下的冰排破壞冰力情況，冰力時程曲線如圖19所示。可以看出，隨著冰排厚度從0.5 m增大到1 m，冰排發生斷裂所用時間也發生了成倍的延長。對于夾角30°的楔形冰排，發生斷裂時的冰力隨厚度增加而增大4倍；45°與60°的楔形冰排，發生斷裂時的冰力隨厚度增加而增大約6倍。冰厚1 m時，不同夾角的冰排受豎向作用力時，冰排內部的破壞機制發生變化；夾角較小時，楔形冰排更趨向于懸臂梁，冰排內發生彎曲破壞的成份較多；夾角較大時，冰排整體受力狀態發生改變，使得冰排不再是簡單的受彎受力狀態，邊界約束的影響逐漸增強，最終導致了斷裂力的增大。

圖19 1 m和0.5 m冰厚下，不同夾角冰排的冰力時程曲線

4.4 加載速度對彎曲破壞冰力的影響

0.5 m厚度冰排，不同夾角楔形冰排在不同加載速度下的彎曲破壞冰力曲線如圖20所示。可以看出，當加載速度較低時，冰力曲線較為平穩，但對應0.3 m/s的冰力曲線在加載初期出現了第一次峰值，這是由于加載板對冰板產生了動態載荷，冰板內的動應力在固定邊界反射后造成冰板在豎直方向上的振動，與加載板發生相向運動時產生了一次冰力峰值。而加載速度達到0.6 m/s時，冰力曲線的震蕩幅度加大，這是由于更大的加載速度造成了更劇烈的冰排振動。可以看出，對于夾角30°的冰排，隨加載速度的增大，冰力增加不明顯；對于45°冰排，加載速度0.6 m/s時的冰力峰值為0.15 m/s時的1.5倍；對于60°冰排，加載速度0.6 m/s時的冰力峰值為0.15 m/s時的2倍。

圖20 冰排在不同加載速度下的冰力比較

5 結 語

本文采用具有粘結-破碎特性的離散元方法，分析了船-冰作用中的擠壓與彎曲兩種典型破壞模式，獲得了不同作用模式下的冰力特征，計算結果可為冰區航行船舶的設計和安全運營提供參考。

與理論求解不同，離散元黏結模型可對材料從變形、裂紋擴展直至完全破壞的全過程進行模擬，可以很好地對海冰材料的破壞過程進行再現。同時，采用海冰力學特性測量中普遍使用的現場單軸壓縮和三點彎曲試驗結果對離散元模型參數進行標定，使得該離散元方法可以更方便地用于船-冰相互作用的研究中。