冰區航行中船體結構冰壓力分布特性的離散元分析

劉 璐，曹 晶，張志剛，季順迎

（1.大連理工大學工業裝備結構分析國家重點實驗室，遼寧大連116024；2.中國船級社上海規范研究所，上海200135）

0 引 言

隨著全球氣候變暖，人類在極區的自然資源勘探考察和開發逐漸增多。其中，北極航道的開辟大大減少了海運航程，也避免了航經不穩定地區帶來的安全性因素。然而，船舶在冰區的運行面臨破冰過程中的結構安全問題，急需研究冰載荷在船體結構表面的分布情況，從而識別船體結構的危險區域并分析結構的冰激疲勞問題，從而確定冰載荷監測的最優位置[1]。

自上世紀80 年代起，國外在極區開展了大量的船體結構冰壓力現場測量工作，積累了大量一年冰和多年冰與船體結構相互作用的相關數據[2]。為研究局部冰壓力的分布特性并從海冰的物理力學性能角度分析冰壓力成因，國外也開展了大量的冰壓力現場試驗[3]。然而，現場的試驗測量經濟成本較高，且現場試驗結果較難獲得重復性驗證，難以獲得更加廣泛的規律。因此，借助于海冰物理力學性質的相關研究成果實現對冰壓力的數值預測，可為深入探究冰壓力的形成機理和冰區船舶的工程設計提供必要的手段[4]。對冰-船相互作用的物理過程進行系統的總結和歸納，形成海冰斷裂過程的物理模型并確定模型參數，同時結合大量的實測數據結果，通過簡化的海冰斷裂模型發展半解析半數值方法預測船體結構冰載荷和冰壓力分布情況[5]。而采用有限元、黏聚單元、離散元等方法分析海冰與結構作用過程則更具有工程實用性，可深入分析不同海冰工況下海冰與結構相互作用過程中的斷裂機理[6]。

近些年來，國內逐步開展了破冰船冰載荷的現場監測和數值計算工作[7]。在海洋工程中使用較多的壓力盒測試法可以直觀地測得結構冰載荷，但其成本較高，不適合船-冰接觸面積較大的船體結構。根據國外的相關經驗，采用時域法中的影響系數矩陣法可通過船體結構上的應變測試結果反演冰載荷[8]。在寒區海洋工程領域，也開展了相應的海冰模型試驗工作[9]。在破冰船冰載荷的數值預測領域，國內開展了較多的工作[10]，在船-冰作用模型、冰壓力預測等方面發展迅速。值得注意的是，在數值模擬中大多采用商業軟件如ABAQUS、LS-DYNA，或將高度經驗化的數學模型引入數值方法中，不能較好地滿足數值預測的廣泛適用性和對復雜條件的兼容性要求。

離散單元法最早用于離散介質的動力過程模擬，也可用于碎冰的堆積成型過程模擬[11]。隨著其粘結模型的不斷發展，該方法也可用于連續介質的破壞過程，在巖土工程、地質體演變等領域應用廣泛[12]。離散元法在海冰領域的應用較早，在冰脊成型、海冰與結構相互作用等方面發展較快[13]。采用擴展圓盤單元可以模擬船體在碎冰區的航行過程，可對碎冰區船舶結構冰載荷的影響因素進行定量化的分析[14]。通過引入基于GPU（Graphics Processing Unit）的并行計算可實現大規模離散元計算，提高對海冰斷裂過程的分析精度和效率[15]。在海冰離散元中，海冰單元的斷裂完全取決于外力的傳遞過程和材料自身參數，體現出了更好的自適應性和健壯性。不同于前述數值模型專注于船體結構的受力分析，海冰離散元方法同時可分析海冰的動力和斷裂過程，包括堆積、阻塞等情況。因此，離散元方法是一種能夠滿足廣泛需求的船-冰相互作用分析及冰載荷預測的方法。

本文通過帶有粘結-破碎模型的海冰離散元方法模擬破冰船與平整冰作用過程，同時采用基于GPU 的大規模并行技術大幅提高計算效率和規模。針對船體結構在平整冰區直行破冰過程，分析了船體結構上的總體冰載荷及局部冰壓力，并通過規范驗證了冰壓力模擬的可靠性。

1 破冰船冰區航行的海冰離散元方法

在平整冰與船體結構的相互作用過程中，將海冰離散為均勻大小和質量的球體單元。海冰單元的受力考慮單元之間粘結力、接觸力以及海水對單元的浮力和拖曳力[16]，這里主要闡述粘結-破碎模型的實施方式。同時，通過對船體結構表面節點力的統計合理建立船體結構表面冰壓力的評估方法。

1.1 海冰單元的粘結-破碎模型

采用平行粘結模型對離散球體單元進行凍結[17]。平行粘結模型假定相鄰的兩個球體單元之間通過圓形截面梁進行粘結，該模型可以同時傳遞粘結力和力矩，如圖1所示。

圖1 球體單元之間的平行粘結模型Fig.1 Parallel bond model between spherical elements

粘結力考慮因法向和切向彈性變形而引起的彈性力模型，這里采用與彈性接觸力模型統一的力模型計算，根據梁模型計算法向和切向的力矩。將力和力矩分為拉力、剪力、扭矩和彎矩，即

式中，F 和M 分別代表力和力矩；下標b 代表粘結，上標n 和s 分別代表法向和切向分量。在粘結模型中考慮梁的拉伸、彎曲和扭轉作用，兩個球形單元之間的最大法向應力和剪切應力為

式中，A、J 和I分別為梁的截面積、極慣性矩和慣性矩，且有A = πR2，J = πR4/2，I = πR4/4，R 為球體單元的半徑。假設拉力為正，壓力為負，采用Mohr-Coulomb理論判斷粘結的斷裂失效，

式中：τb即為切向粘結強度；C 為材料粘聚力強度；σ 是法向應力；μ表示內摩擦系數，μ = tanθ，θ 表示曲線傾角，如圖2 所示。粘結的失效準則可以表示為

圖2 海冰離散元中的Mohr-Coulomb失效準則Fig.2 Mohr-Coulomb criterion in DEM for sea ice

式中，σb為法向拉伸強度。圖2中實線外側即為斷裂區域，即上式表示的失效準則。

在海冰離散元的實際應用中，如何根據已知海冰的強度等力學性質合理選擇離散元參數是數值模擬準確性的關鍵。通過單軸壓縮和三點彎曲試驗的離散元模擬，研究顆粒單元尺寸、粘結強度和內摩擦系數對離散元模擬結果的影響，進而可確定海冰宏觀強度與以上三者之間關系[18]。海冰的彎曲強度與顆粒單元的粘結強度緊密相關，可通過對海冰的三點彎曲試驗確定海冰宏觀強度σf與離散元粘結失效參數σb之間的關系，

式中，h 為海冰厚度，D 為顆粒單元直徑。另外，一般在具體的模擬中，令μ = 0.2 且σb= C[15]。以上關系式與渤海現場試驗中測得的海冰彎曲強度對比驗證結果良好[19]，可用于渤海海冰力學性質的模擬分析，是海冰離散元參數選取的重要依據。

1.2 船體結構的冰壓力計算

將船體的外殼結構劃分為若干三角形單元。通過海冰單元對船體結構的碰撞力插值到三角形單元的三個節點上，并考慮三角形單元的面積，即可得到每個節點上的冰壓力。

由于離散元方法的積分步長非常小，通常在10-9～10-5s 范圍內，直接統計在某個時間步上的瞬時冰壓力顯然是不合理的。這里通過統計平均的方式，將一段時間內的冰壓力平均值作為這段時間的冰壓力，那么第i個統計時間步某一節點上的冰壓力可表示為

式中，pj是第j 個時間步的冰壓力，r 和s 分別代表該段時間的起始和結束時間步數。式中表示的冰壓力即為瞬時冰壓力，它體現了船體結構在較短時間內的冰壓力情況。但是船體結構的瞬時冰壓力不能很好地體現航行過程中的統計特性。為了克服瞬時冰壓力在表述上的缺點，更好地體現船體結構上冰壓力較大且作用頻率較高的區域，即高壓和易疲勞區域，這里引入累計最大冰壓力和累計平均冰壓力。第N個統計時間步某一節點上的累計最大冰壓力和累計平均冰壓力可表示為

從船體結構上的所有節點來看，累計最大冰壓力體現了船體結構上的高壓危險位置，該區域可能會產生結構強度破壞；累計平均冰壓力體現了船體結構上冰壓力較大且作用頻率較高的位置，會產生較為持續且強度較高的冰壓力作用，容易產生材料的疲勞強度破壞。

2 破冰船在平整冰區的離散元模擬

采用海冰離散元模擬船體結構在平整冰區破冰時與海冰的相互作用過程，分析船舶破冰過程中冰壓力在船體結構表面的分布形式。船體結構模型采用某極地科學考察船，計算中將船體結構視為剛體，船體長、寬和高分別為167 m、22.6 m 和13.5 m，吃水深度為9 m。圖3 所示為該科學考察船船體表面劃分為三角形的結構形式，模擬中三角單元的平均尺寸為1.5 m。

計算中主要的計算參數列于表1 中。其中，海冰的彎曲強度是海冰物理力學參數的重要表征，直接影響到海冰的破壞模式，進而影響海冰與結構物相互作用過程中的冰載荷[19]。在海冰與結構相互作用的離散元模擬中，需考慮將海冰的宏觀強度轉化為離散元中球體單元間的微觀強度，確保模擬中海冰的宏觀強度與實際保持一致，該過程中還需綜合考慮單元的尺寸效應[20]。海冰受到海流的拖曳力模型參考Sun和Shen（2012）[16]的相關工作。

圖3 科學考察船的船體結構三角單元劃分Fig.3 Triangular meshes of ship hull for the scientific research ship

表1 船體結構破冰過程模擬的離散元參數Tab.1 DEM parameters in the simulation of icebreaking process of ship hull

2.1 船體結構破冰過程的離散元分析

連續的直行破冰是破冰船冰區航行中的主要破冰操作模式，這里采用離散元方法分析船體結構直行破冰過程中的總體冰載荷和局部冰壓力，主要計算參數列于表2中。在模擬中，平整冰區計算域的長×寬為600 m×200 m。計算域邊界和海冰單元之間采用彈簧約束，彈簧剛度與單元之間的接觸剛度相同。這里參考ISO19906，取極地海冰的彎曲強度σf=0.7 MPa[21]。當海冰厚度為1.5 m 時，兩層球體單元的總單元數約為41 萬，且球體單元粒徑為0.826 m。海冰彎曲強度為0.7 MPa 時，根據式（7）算得單元粘結強度為1.315 MPa。采用基于GPU 的并行技術，破冰船的直行破冰過程模擬如圖4 所示。從圖中可以看出，在船體結構的作用下海冰發生斷裂，在冰面上形成邊緣不規則的斷裂從而形成開闊水道。

表2 船體結構直行破冰模擬中的相關計算參數Tab.2 Main parameters in the simulation of ship icebreaking in straight line

圖4 船體結構破冰過程的離散元模擬Fig.4 DEM simulation of the icebreaking process of ship hull

圖5 船體結構破冰過程中的海冰破壞模式Fig.5 Failure mode of sea ice during icebreaking of ship hull

圖5 是船體結構破冰過程中的海冰破壞模式對比。從圖中可看出，離散元模擬中船體結構兩側出現了“鋸齒”狀連續的斷裂，這種破壞模式是典型的彎曲破壞。在實際船舶走航的海冰斷裂監測中，如雪龍號和芬蘭某型冰區船，也出現了類似的破壞模式。但是實船監測中海冰還呈現“環狀”的破壞，并不是規整的“鋸齒”狀。在離散元模擬中船體結構是剛體，只有x方向的剛性平動；而在實船破冰過程中，由于船體兩側非同時的海冰結構和破壞，船體還會產生橫搖，從而對船體兩側的海冰造成擠壓作用。因此，在實船破冰過程中海冰還會發生剪切破壞，其“環狀”的裂紋是剪切和彎曲破壞等復雜作用下的結果[22]。

船體結構破冰過程中的總體冰載荷如圖6 所示。從圖中可以看出，x 和z 方向的冰載荷具有類似的特點，即在較高的水平線上下波動，這種波動體現了海冰的斷裂過程，同時說明了船體結構運動過程中海冰的持續作用。y 方向的冰載荷在0 附近上下波動，且會出現較大的峰值載荷，說明了船舶在破冰過程中海冰會從兩側對結構產生擠壓，而這種擠壓作用會對結構安全性造成嚴重隱患。因此船舶兩側的冰載荷影響還需要對局部的冰壓力進行深入的分析。

圖6 船體結構上的總體冰載荷時程Fig.6 Time history of ice loads on ship hull

船體結構破冰過程中的局部冰壓力如圖7 所示。從圖中可以看出，船舶破冰過程中最大的冰壓力主要集中在船肩處，該部位應是船體結構安全的重點監測部位。由于破冰部位主要為船艏柱，該部位會與海冰發生較多的接觸，所以累計平均冰壓力集中于此，也說明船艏柱是冰激結構疲勞分析重點關注的部位。

圖7 船體結構上的冰壓力分布Fig.7 Distribution of ice pressure on ship hull

2.2 船體結構冰壓力的對比驗證

國際船級社組織（International Association of Classification Societies,IACS）對特定工況下船體結構的冰壓力給出了規范計算方法[23]。該計算針對船體結構與大塊浮冰的自由碰撞過程，通過能量守恒原理并結合一定的經驗處理，可總結該碰撞過程的理論模型。在理論模型中根據碰撞點處的船型參數給出了等效的壓力板大小，根據該壓力板上的冰載荷即可求出碰撞點上的冰壓力[24]。通過離散元計算了船體結構與大塊浮冰的碰撞過程，分析了碰撞點上的冰壓力并與IACS 的規范進行對比，從而驗證本文離散元方法的合理性。采用IACS 冰級PC5的相應參數作為離散元計算參數，主要參數列于表3中，其他相關參數參照表1。其中，船速為2.25 m/s，海冰厚度為3.0 m，海冰彎曲強度為1.0 MPa[25]。根據式(7)可由海冰厚度和彎曲強度確定粘結失效強度的具體數值。船體結構與大塊浮冰自由碰撞的離散元模型如圖8所示。設置船體結構與海冰的碰撞點為距船艏柱5 m、15 m、25 m和35 m位置，并分別根據規范推薦方法計算對應碰撞點上的冰壓力。

表3 船體結構冰壓力驗證分析中的相關計算參數Tab.3 Main parameters in the validation of ice pressure on ship hull

圖8 冰壓力規范驗證的離散元模型Fig.8 Sketch of DEM simulation for the validation with IACS standard

船體結構與大塊浮冰碰撞的離散元模擬過程如圖9 所示。從圖中可以看出，船艏部位與海冰發生碰撞造成海冰的擠壓破碎。隨著船體結構的運動，海冰不斷發生破碎且破碎后的海冰呈粉末狀浮在水面上。圖10 為四個算例中碰撞點處的冰壓力與規范對比情況。其中，x 為碰撞點距船艏柱的距離，alpha和beta是碰撞點處的船型參數，w和b為根據規范計算的壓力板寬度和高度，P為離散元模擬的冰壓力，Pstd為IACS規范值。

圖9 船體結構與大塊浮冰碰撞的離散元模擬Fig.9 DEM simulation of collision between ship hull and ice floe

圖10 碰撞點處船體結構的冰壓力Fig.10 Ice pressure on ship hull on the collision position

離散元模擬結果與規范進行對比的情況列于表4 中。離散元的計算結果與規范值相比存在上下誤差，壓力值保持在相同量級范圍內，誤差范圍在6.7%～18.1%之間。因此，本文的離散元方法對船體結構的冰壓力可進行準確可靠的分析，具有良好的工程適用性。

表4 離散元模擬結果與IACS標準對比Tab.4 Comparison between DEM simulation and IACS standard

3 結 語

本文基于平行粘結模型建立了基于球體單元的海冰離散元方法，并通過Mohr-Coulomb 準則實現了單元間的斷裂，從而可對海冰與結構相互作用的破碎過程進行模擬；采用海冰離散元方法模擬了科學考察船的破冰過程，研究了船體結構上的總體冰載荷時程；同時，結合船舶破冰的特點，對累計最大冰壓力和累計平均冰壓力進行了有效的分析；為驗證本文計算結果的可靠性，通過IACS 規范對本文離散元方法計算的船體結構冰壓力進行了校核。結果表明，本文離散元方法的冰壓力計算誤差范圍在6.7%～18.1%之間，充分說明了該方法的合理性。