冰脊壓剪試驗及其對直立結構冰載荷的離散元分析1)

朱紅日 季順迎

(大連理工大學工業裝備結構分析國家重點試驗室，遼寧大連116024)

冰脊廣泛分布于渤海、波羅的海、極區等海域，通常由大面積浮冰相互擠壓、剪切進而發生破碎、堆疊而成[1-2]。典型的冰脊一般包括上部的脊帆、水面附近的凍結層及水下的龍骨三個部分，厚度可達5～30 m[3]。冰脊作用在海洋結構造成的冰載荷遠大于平整冰，很多情況下是船舶與海洋工程結構的極限設計載荷。因此，對冰脊的力學性質及其對海洋工程結構冰載荷的研究具有重要的工程意義。

現場測量與室內模型試驗一直是研究冰脊的主要手段。在波弗特海的Molikpaq 平臺[4]、加拿大的聯邦大橋[5]及波的尼亞灣的Norstrmsgrund 燈塔[6]等海洋工程結構上開展了大量的冰脊載荷現場測量；此外，在現場及實驗室內也開展了一系列的直剪[7]、雙軸壓縮[8]、壓剪[7,9-10]等力學性質試驗以研究冰脊特別是龍骨部分的力學特性。文獻[11-12] 基于不同的冰脊破壞模式發展了相應的冰載荷計算模型，但目前不同的方法存在很大離散性。

近年來，有限元和離散元(discrete element method，DEM)等數值方法能計算復雜工況的優勢日益顯現，成為冰脊載荷研究的重要途徑。Heinonen[10]和 Serr[7]等建立了冰脊龍骨的本構模型，采用有限元法對冰脊的力學特性進行了系統的數值計算。Hopkins[13]采用二維塊體離散元方法模擬了冰脊的形成過程，發現冰脊在形成過程中存在巨大的能量耗散；Polojrvi 等[14-15]發展了三維塊體離散元方法并計算了部分凍結及未凍結條件下冰脊龍骨的壓剪試驗，分析了加載速度對壓剪試驗結果的影響；Yulmetov 等[16]利用由球體組合而成的冰塊單元計算了冰脊與錐體結構的相互作用過程以及海冰堆積現象。在工程計算中，冰脊內部冰塊的尺寸一般遠小于冰脊尺寸，冰塊的形狀對冰脊整體力學特性的影響可以忽略，由此可將冰塊簡化為球體單元從而兼顧計算效率及準確性。

為此，本文采用具有粘結?破壞功能的球體單元構造冰脊模型，通過計算現場壓剪試驗過程以分析計算模型的合理性及主要計算參數對冰脊力學特性的影響；在此基礎上進一步計算冰脊與圓形直立結構的相互作用，并通過與相關規范進行對比以驗證該離散元方法的可靠性。

1 冰脊的離散元方法

粘結球體離散單元方法由Cundall 和Potyondy于2004 年提出[17]，目前已廣泛用于巖石、混凝土、玻璃等破壞特性的數值計算，其在模擬材料脆性破壞及力學行為方面具有顯著優勢。冰脊由大量的碎冰塊堆積凍結而成，可采用該離散元方法進行構造和力學性質的數值計算。

1.1 冰脊結構的離散元構造

冰脊自上而下分為脊帆、凍結層、龍骨三個部分，通常將脊帆及龍骨理想化為梯形結構，如圖1(a)和圖1(b) 所示[18]。凍結層部分由于已完全凍結，可采用平整冰的離散元模型由相同粒徑的球體單元按密六方排列構造[19]。龍骨與脊帆由部分凍結或未凍結的冰塊堆積而成，其內部結構較為復雜，可由不同粒徑的球體單元隨機排列構造，如圖 1(c) 所示。本文采用一種球體幾何隨機排列方法構造龍骨和脊帆的離散元模型[20]。在冰脊離散元方法中，根據冰塊間的凍結狀態確定單元接觸對間是否存在粘結力。接觸力的計算將單元間的法向力簡化為彈簧和阻尼器，并將切向力簡化為彈簧、阻尼器和滑動器，根據顆粒間的相對運動狀態進行計算。當顆粒接觸對處于粘結狀態時，采用平行粘結模型[21]模擬冰塊間凍結作用，如圖2 所示。

圖1 冰脊的幾何結構及其離散元模型

圖2 冰脊球體單元間的平行粘結模型

在平行粘結模型中，粘結單元間存在一個虛擬的彈性圓盤，由此傳遞兩個單元間的軸向力Fn、剪力Fs、彎矩Mt和扭矩Mb。圓盤的最大拉應力σmax和剪應力τmax根據梁單元模型進行計算，可表示為[17]

式中，A為圓盤的截面積，I為圓盤的慣性矩，J為圓盤的極慣性矩，為圓盤的半徑。

當σmax或τmax超過其相應的強度極限時，粘結面發生破壞，其中剪切破壞強度σt基于摩爾?庫倫準則確定，即

式中，σc表示切向粘結強度，μp是顆粒間內摩擦系數。海冰單元間的粘結強度是離散元計算中的重要參數，其與海冰的溫度、鹽度、加載速率等因素密切相關。

基于上述理論構造而成的冰脊模型中，脊帆部分在冰脊中占比很小，其離散元參數對冰脊載荷的影響可以忽略不計。而凍結層與平整冰的力學性質較為接近，其離散元參數可直接參考平整冰的校準方法[18]。因此，本文主要針對龍骨部分離散元參數進行研究。龍骨部分的結構十分龐大且位于水下，無法采用試驗手段直接測量其力學參數，目前主要通過壓剪實驗等方法間接研究其力學特性。這里主要參考Heinonen 于2004 年在波的尼亞灣開展的現場壓剪試驗[10]進行離散元計算(圖3)，分析離散元參數對龍骨力學性質的影響并進行計算參數校準。

圖3 冰脊壓剪試驗[10] 及離散元計算模型

1.2 冰脊壓剪過程的離散元分析

在 Heinonen 的冰脊龍骨壓剪試驗中，首先去除龍骨最深處上方附近的脊帆及凍結層部分，然后用一個圓形壓板緩慢豎直作用于龍骨上部，直至壓板下方的冰塊與冰脊徹底脫離，如圖3(a) 所示[10]。試驗中主要記錄壓板載荷及運動時程數據，同時在水底放置攝像頭記錄龍骨底部的變形及破壞現象。圖3(b) 為參照其現場試驗建立的離散元計算模型，其中龍骨部分由于尺寸遠大于壓盤，其可以直接近似為長方體，且其邊界可設為固定約束。該計算中的主要參數列于表1 中。龍骨的單元粒徑取0.2～0.4 m，由此生成的龍骨計算模型孔隙率為0.31。

圖4 冰脊壓剪試驗中離散元數值及現場測量的載荷與位移關系

表1 冰脊壓剪試驗中主要計算參數

通過壓剪試驗的離散元計算得到的壓板載荷隨位移的變化曲線如圖4 所示。對比數值結果與試驗數據可以發現，無論在離散元計算還是現場試驗中壓板載荷都經歷了快速線性增長至最大值、較平緩下降、保持在與壓板下碎冰堆所受浮力平衡的位置等幾個階段。在ABC的線性加載階段，可以由力和位移關系的斜率確定冰脊的剛度。龍骨在計算與試驗中的整體剛度基本一致，峰值大小也非常接近。

為進一步分析龍骨在壓剪試驗中的破壞過程，圖5 給出了四個不同時刻的冰脊變形狀態。左側為不同時刻的龍骨破壞狀態、中間及右側分別為豎直、水平方向的位移云圖。從龍骨的破壞狀態可以發現，在u=10.0 mm 的加載初期，龍骨呈局部壓縮變形，并主要集中在壓板正下方很小范圍內；隨著壓板的向下運動，當u=15.0 mm 時，龍骨內部的壓縮變形區域不斷增大，受擠壓作用的碎冰塊不斷水平向外擴散引起周圍冰塊移動形成一個楔形變形區域，楔形部分整體下移并向兩側擴張在邊緣處形成帶狀剪切變形區域；當u=25.0 mm 時，龍骨楔形部分兩側的剪切變形達到極限，龍骨發生整體破壞，形成明顯的剪切面，壓板載荷開始下降；在u=100.0 mm 的卸載過程中，由于楔形部分受壓膨脹與周圍龍骨部分發生摩擦，故卸載過程較加載更為平緩。上述過程說明在數值計算中，龍骨呈現出剪切破壞和壓縮破壞的混合破壞模式。這一現象與現場試驗結果相一致。由于條件限制，現場試驗并未取得龍骨破壞過程的內部圖像，但通過對比計算結果與室內模型試驗結果，可以發現其主要現象基本一致[22-23]。

圖5 不同位移下龍骨內部的破壞狀態過程及位移云圖

1.3 離散元參數對壓剪強度的影響

在離散元數值計算中，計算參數對冰脊的壓剪失效特性有顯著影響。這里選擇對冰脊龍骨力學性質影響較大的單元彈性模量Ep、單元間粘結強度σc以及內摩擦系數μp三個參數進行研究。通過不同參數組合下的數值壓剪試驗，研究以上參數對壓剪試驗的影響。為分析海冰單元彈性模量Ep的影響，首先在表1 所列計算參數的基礎上，取σc=50 kPa，Ep分別取 0.3，0.5，1.0，1.5，2.0 GPa 進行冰脊壓剪過程的離散元模擬，由此得到的最大壓力如圖6 所示。計算結果表明，Ep在 0.3～1.0 GPa 范圍內，Fmax幾乎保持不變；而當Ep>1.0 GPa 后，Fmax隨著Ep的增大而線性下降。這是由于當Ep<1.0 GPa 時，龍骨在試驗中保持以全局性剪切破壞為主的破壞形式，Fmax主要受強度參數σc等影響，與Ep關系不大。Ep主要影響龍骨的整體剛度，與其呈線性正相關關系，如圖7 所示；當Ep>1.0 GPa 時，隨著Ep的不斷增大，龍骨的破壞形式逐漸轉變為以連續的局部壓縮破壞為主，壓板載荷峰值變得密集，最大值Fmax不斷減小，如圖8 所示。

圖6 壓剪破壞過程中壓板最大載荷與彈性模量Ep 的對應關系

圖7 離散元模擬的冰脊整體剛度與彈性模量Ep 的對應關系

圖8 不同彈性模量Ep 下壓板載荷隨豎直位移的變化過程

圖9 離散元模擬的最大載荷與粘接強度σc 的對應關系

為分析海冰單元間粘接強度σc的影響，當Ep= 0.5 GPa，σc分別取 60，120，200，300 kPa時，離散元計算得到的壓板載荷最大值Fmax如圖9所示。可以看出，在所取σc值范圍內，Fmax一直隨σc的增大而線性增長。但從不同粘結強度下壓板載荷達到最大值時龍骨內部裂紋分布情況可以發現 (如圖 10 所示)，當σc<120 kPa 時，龍骨保持現場試驗中常見的全局性剪切破壞為主的破壞形式；當σc>120 kPa 時，隨著σc的增大，壓板下方自底部向上方擴展的豎向裂紋逐漸明顯；當σc=300 kPa時，龍骨截面只有豎直裂紋，沒有明顯的剪切帶。這說明隨著σc的增大，龍骨的破壞形式逐漸由剪切破壞為主轉變為彎曲破壞為主。這也與一些室內模型試驗中觀測到的現象相一致[24]。

圖10 不同σc 下離散元模擬的冰脊的破壞模式

為分析海冰單元間摩擦系數μp的影響，取Ep= 0.5 GPa 和σc= 50 kPa，當μp分別為0.2，0.4，0.6，0.8 時，壓板載荷最大值時龍骨內部裂紋分布如圖11 所示。從中可以發現，這些龍骨在壓板作用下均以整體剪切破壞為主，其楔形斷裂部分兩側均分布有明顯的剪切帶，且兩側剪切帶與數值方向夾角θ隨μp的增大而線性增大。計算結果還進一步表明，壓板載荷的最大值Fmax不僅隨σc的增大而線性增大，也與μp呈線性正相關關系，如圖12 所示。這表明冰脊在壓剪過程中表現出明顯的摩爾?庫倫材料特征。

從以上離散元數值模擬結果可以看出，單元彈性模量Ep、單元間粘結強度σc以及內摩擦系數μp均對龍骨在壓剪試驗中的力學行為有重要影響。這些參數不僅影響其壓剪強度的大小，還決定其破壞形態甚至破壞模式。基于對這些參數影響的分析和現場試驗結果，可以選擇合理的離散元參數組合對冰脊進行數值計算以得到可靠的結果。

圖11 不同μp 下冰脊的失效模式

圖12 不同摩擦系數μp 下離散元模擬的最大載荷及失效夾角

2 冰脊與圓形直立結構相互作用的離散元分析

2.1 冰脊與圓形直立結構相互作用的計算模型

冰脊龍骨部分采用表 1 的離散元參數，將凍結層與平整冰簡化為統一厚度，其他主要參數列于表2 中。冰脊的凍結層由球體單元按密六方排列形式構造，如圖13 所示，其厚度hc為

式中，hc為凍結層厚度，d為單元直徑，ne為單元層數。

表2 冰脊與直立結構相互作用中離散元計算的主要參數

圖13 由球體單元密六方排列而成的冰脊凍結層

圖14 和圖15 為結構與冰脊作用過程圖，圖中顆粒顏色代表速度的大小，而圖16 為冰脊對直立結構在水平方向上的冰載荷。從中可以發現，冰脊在整個作用過程中以結構前的局部擠壓破壞為主。隨著圓形直立結構不斷貫入冰脊內部，結構前的碎冰堆積深度不斷增大，導致冰載荷也隨之增大；當冰脊龍骨發生整體剪切破壞后，冰載荷與結構前碎冰堆積深度開始減小直至龍骨部分冰載荷可以忽略不計，而凍結層則始終保持擠壓破壞。由圖16 還發現直立結構冰載荷最大值為9.53 MN，凍結層冰載荷最大值為6.50 MN。由此可見，冰脊中的凍結層是對結構作用的主要承載部分，而龍骨的冰載荷相對較小。

圖14 直立結構與冰脊相互作用離散元模擬

圖15 直立結構與冰脊作用過程的側視圖

圖16 離散元模擬的冰脊對直立結構的冰載荷

2.2 冰脊對圓形直立結構冰載荷的驗證

為驗證計算結果的準確性，這里選擇將離散元結果與基于 ISO 19906 規范[18]的計算結果進行對比。ISO 19906 將冰脊載荷Fr分為凍結層載荷Fc與龍骨載荷Fk兩部分，即

式中，Fc按照平整冰經驗公式計算

式中，CR按照規范推薦取 1.8 MPa，h1取 1 m，hc為凍結層厚度，D為圓柱結構直徑，n=?0.5+hc/5，m為?0.16。Fk的計算基于改進的Dolgopolov 公式，即

式中，D為圓柱結構直徑；hk為龍骨深度；η為龍骨部分孔隙率；ρw為海水密度；ρice為海冰密度；g為重力加速度；c為宏觀上冰脊龍骨表現出的粘聚力，這里取12 kPa；?為龍骨內摩擦角，結合壓剪試驗結果并基于Meyerhof 理論求得[25-26]

式中，Dind為壓盤直徑；Ku為被動冰壓力系數，一般為常數0.95；Qu為壓剪載荷最大值756 kN；R為壓剪試驗中龍骨完全破壞后的殘余壓力120 kN。

由此可計算冰脊總載荷Fr、凍結層載荷Fc及龍骨載荷Fk。這里將離散元和ISO 19906 規范結果列于表 3。從中可以發現，離散元計算結果中的Fc相對ISO 19906 規范較小，而Fk則相對較大。這主要是由于ISO 19906 規范對凍結層冰載荷的計算是基于大量平整冰的試驗數據并取其較大數值，其相對比較保守以保障結構安全；而在龍骨冰載荷的計算中，離散元方法未具體考慮碎冰塊的隨機排列和冰塊尺寸影響，而是按龍骨中的冰塊充分凍結進行計算，從而導致其計算結果略高于ISO 19906 規范。此外，在離散元模擬中，單元間細觀計算參數選取的合理性以及冰載荷動力過程的隨機性也是影響冰載荷大小的重要因素。但從直立結構冰載荷的總體情況看，離散元計算結果與ISO 19906 規范比較一致，從而說明了離散元方法在計算冰脊對直立結構冰載荷中的適用性。

表3 冰脊對圓柱結構載荷的離散元計算與 ISO 19906 規范對比

3 結論

本文基于粘結球體單元構建了冰脊離散元模型，通過顆粒排列方式與離散元參數的變化表現冰脊凍結層、龍骨等不同部分的力學特征。利用該模型計算分析了一系列冰脊龍骨壓剪試驗，并與現場實測數據進行了對比。計算結果表明：在合理的離散元參數下該方法能較為準確地模擬冰脊的力學行為；單元彈性模量Ep、單元間的粘結強度σc以及內摩擦系數μp不僅影響壓板載荷的大小，還決定龍骨受壓后的破壞形式。在一定參數范圍內龍骨主要以全局性剪切破壞為主，壓板載荷會隨σc和μp的增大而線性增長，且剪切面與垂直方向的夾角也會隨μp的增大而增大。最后計算了冰脊與圓柱形直立結構的相互作用過程，其計算結果與ISO 19906 規范相近，由此證明了該方法在工程計算中的可靠性。離散元方法能模擬冰脊與直立結構作用中的堆積效果與破壞模式，在更復雜的工況計算中將具有良好的應用前景。