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

（1.大連理工大學(xué)工業(yè)裝備結(jié)構(gòu)分析國家重點實驗室，遼寧大連116024；2.中國船級社上海規(guī)范研究所，上海200135）

0 引 言

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

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

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

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

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

1 破冰船冰區(qū)航行的海冰離散元方法

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

1.1 海冰單元的粘結(jié)-破碎模型

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

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

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

式中，F(xiàn) 和M 分別代表力和力矩；下標(biāo)b 代表粘結(jié)，上標(biāo)n 和s 分別代表法向和切向分量。在粘結(jié)模型中考慮梁的拉伸、彎曲和扭轉(zhuǎn)作用，兩個球形單元之間的最大法向應(yīng)力和剪切應(yīng)力為

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

式中：τb即為切向粘結(jié)強度；C 為材料粘聚力強度；σ 是法向應(yīng)力；μ表示內(nèi)摩擦系數(shù)，μ = tanθ，θ 表示曲線傾角，如圖2 所示。粘結(jié)的失效準(zhǔn)則可以表示為

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

式中，σb為法向拉伸強度。圖2中實線外側(cè)即為斷裂區(qū)域，即上式表示的失效準(zhǔn)則。

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

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

1.2 船體結(jié)構(gòu)的冰壓力計算

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

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

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

從船體結(jié)構(gòu)上的所有節(jié)點來看，累計最大冰壓力體現(xiàn)了船體結(jié)構(gòu)上的高壓危險位置，該區(qū)域可能會產(chǎn)生結(jié)構(gòu)強度破壞；累計平均冰壓力體現(xiàn)了船體結(jié)構(gòu)上冰壓力較大且作用頻率較高的位置，會產(chǎn)生較為持續(xù)且強度較高的冰壓力作用，容易產(chǎn)生材料的疲勞強度破壞。

2 破冰船在平整冰區(qū)的離散元模擬

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

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

圖3 科學(xué)考察船的船體結(jié)構(gòu)三角單元劃分Fig.3 Triangular meshes of ship hull for the scientific research ship

表1 船體結(jié)構(gòu)破冰過程模擬的離散元參數(shù)Tab.1 DEM parameters in the simulation of icebreaking process of ship hull

2.1 船體結(jié)構(gòu)破冰過程的離散元分析

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

表2 船體結(jié)構(gòu)直行破冰模擬中的相關(guān)計算參數(shù)Tab.2 Main parameters in the simulation of ship icebreaking in straight line

圖4 船體結(jié)構(gòu)破冰過程的離散元模擬Fig.4 DEM simulation of the icebreaking process of ship hull

圖5 船體結(jié)構(gòu)破冰過程中的海冰破壞模式Fig.5 Failure mode of sea ice during icebreaking of ship hull

圖5 是船體結(jié)構(gòu)破冰過程中的海冰破壞模式對比。從圖中可看出，離散元模擬中船體結(jié)構(gòu)兩側(cè)出現(xiàn)了“鋸齒”狀連續(xù)的斷裂，這種破壞模式是典型的彎曲破壞。在實際船舶走航的海冰斷裂監(jiān)測中，如雪龍?zhí)柡头姨m某型冰區(qū)船，也出現(xiàn)了類似的破壞模式。但是實船監(jiān)測中海冰還呈現(xiàn)“環(huán)狀”的破壞，并不是規(guī)整的“鋸齒”狀。在離散元模擬中船體結(jié)構(gòu)是剛體，只有x方向的剛性平動；而在實船破冰過程中，由于船體兩側(cè)非同時的海冰結(jié)構(gòu)和破壞，船體還會產(chǎn)生橫搖，從而對船體兩側(cè)的海冰造成擠壓作用。因此，在實船破冰過程中海冰還會發(fā)生剪切破壞，其“環(huán)狀”的裂紋是剪切和彎曲破壞等復(fù)雜作用下的結(jié)果[22]。

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

圖6 船體結(jié)構(gòu)上的總體冰載荷時程Fig.6 Time history of ice loads on ship hull

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

圖7 船體結(jié)構(gòu)上的冰壓力分布Fig.7 Distribution of ice pressure on ship hull

2.2 船體結(jié)構(gòu)冰壓力的對比驗證

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

表3 船體結(jié)構(gòu)冰壓力驗證分析中的相關(guān)計算參數(shù)Tab.3 Main parameters in the validation of ice pressure on ship hull

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

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

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

圖10 碰撞點處船體結(jié)構(gòu)的冰壓力Fig.10 Ice pressure on ship hull on the collision position

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

表4 離散元模擬結(jié)果與IACS標(biāo)準(zhǔn)對比Tab.4 Comparison between DEM simulation and IACS standard

3 結(jié) 語

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