劉 璐，胡 冰 ，季順迎

(1. 大連理工大學(xué) 工業(yè)裝備結(jié)構(gòu)分析國家重點(diǎn)實(shí)驗(yàn)室，遼寧 大連 116024；2. 上海交通大學(xué) 海洋工程國家重點(diǎn)實(shí)驗(yàn)室，上海 200240；3. 中國遠(yuǎn)洋海運(yùn)集團(tuán)有限公司，上海 200027)

北極航道大大減少了我國與歐洲的海運(yùn)航程，同時(shí)避免了航行經(jīng)過不穩(wěn)定地區(qū)帶來的安全性問題。采用極地船舶獨(dú)立進(jìn)行冰區(qū)航行會(huì)大大提高航運(yùn)成本，而采用破冰船引航的方式則可有效降低成本[1]。船舶在冰區(qū)的運(yùn)行面臨破冰過程中的結(jié)構(gòu)安全問題，對(duì)船舶設(shè)計(jì)和冰區(qū)航行保障的要求較高，需采用有效方法對(duì)冰荷載進(jìn)行分析并總結(jié)不同航行條件下的冰荷載特點(diǎn)[2]。

目前，冰荷載的研究方法主要包括現(xiàn)場監(jiān)測、模型試驗(yàn)、經(jīng)驗(yàn)公式和數(shù)值模擬。國外開展船舶結(jié)構(gòu)冰荷載的現(xiàn)場測量工作較早，積累了大量復(fù)雜冰況的冰荷載數(shù)據(jù)[3]。對(duì)船冰相互作用的物理過程進(jìn)行系統(tǒng)總結(jié)和歸納并結(jié)合大量的實(shí)測數(shù)據(jù)結(jié)果，可建立經(jīng)驗(yàn)化的直接計(jì)算公式[4]。近些年來，國內(nèi)逐步開展了船舶結(jié)構(gòu)冰荷載分析工作[5-6]，在冰荷載現(xiàn)場原型測試、數(shù)值模擬等領(lǐng)域開展了大量工作[7-8]。數(shù)值模擬方法具有較好的參數(shù)可控性、結(jié)果可重復(fù)性，經(jīng)濟(jì)成本也較低，還可對(duì)海冰破碎、冰塊碰撞重疊等過程進(jìn)行細(xì)化分析，其中離散元方法使用較為廣泛[9]。離散元方法可對(duì)船冰相互作用過程中的海冰斷裂機(jī)理進(jìn)行深入分析，通過引入流固耦合模型可模擬船冰作用過程中流體對(duì)海冰和船舶結(jié)構(gòu)的影響[10]。

基于閔可夫斯基和(Minkowski sum)的擴(kuò)展多面體單元方法可準(zhǔn)確描述海冰的非規(guī)則形態(tài)，對(duì)海冰與結(jié)構(gòu)的碰撞過程中的冰載荷進(jìn)行有效分析[11]。在海冰離散元方法中，海冰單元的運(yùn)動(dòng)和斷裂完全決定于單元間的作用力傳遞和材料的自身參數(shù)，體現(xiàn)出更好的自適應(yīng)性和健壯性。本文采用基于閔可夫斯基和原理構(gòu)造擴(kuò)展多面體單元，建立擴(kuò)展多面體離散元方法。采用該方法計(jì)算海冰與船舶結(jié)構(gòu)的相互作用過程，分析船舶結(jié)構(gòu)的冰載荷并與經(jīng)典的Lindqvist公式和Riska公式對(duì)比驗(yàn)證結(jié)果的可靠性。采用離散元方法模擬破冰船引航條件下的船冰作用過程，分析不同航速和船寬比條件下的冰載荷特點(diǎn)。

1 海冰的擴(kuò)展多面體離散元方法

根據(jù)閔可夫斯基和原理構(gòu)造擴(kuò)展多面體單元，采用Voronoi切割算法構(gòu)造擴(kuò)展多面體形態(tài)的海冰單元。同時(shí)，發(fā)展了擴(kuò)展多面體單元間的粘結(jié)失效模型，并采用非線性的接觸力模型計(jì)算粘結(jié)失效后單元間接觸力，從而建立海冰的擴(kuò)展多面體離散元方法。

1.1 海冰的擴(kuò)展多面體單元構(gòu)造

對(duì)于空間中的兩個(gè)幾何體A和B，二者的閔可夫斯基和可以表示為：

采用2D-Voronoi切割算法劃分平面單元，在此基礎(chǔ)上生成具有一定厚度的海冰單元。Voronoi單元的幾何形態(tài)可通過規(guī)則度參數(shù)進(jìn)行控制[12]。如圖2所示，首先將2D-Voronoi單元按照擴(kuò)展半徑進(jìn)行“腐蝕”形成新的多邊形，通過給定厚度生成具有一定厚度的三維多面體單元，最后通過閔可夫斯基和運(yùn)算生成擴(kuò)展多面體單元。

圖 1 基于閔可夫斯基和的擴(kuò)展多面體Fig. 1 Dilated polyhedron based on Minkowski sum

圖 2 2D-Voronoi切割算法生成海冰單元Fig. 2 Sea ice element based on 2D-Voronoi tessellation algorithm

1.2 擴(kuò)展多面體的粘結(jié)失效模型

考慮粘結(jié)作用在單元間的分布，本文采用在單元交界面上設(shè)置粘結(jié)點(diǎn)的方式粘結(jié)相鄰且共面的單元。如圖3所示，在粘結(jié)面上設(shè)置粘結(jié)點(diǎn)，粘結(jié)點(diǎn)之間的法向應(yīng)變?chǔ)舗可寫作：

式中：d為兩粘結(jié)節(jié)點(diǎn)間的距離向量；n為交界面的法向，這里取為兩單元粘結(jié)面的外法向之差，即n =unit(ni- nj)；Cij為兩個(gè)粘結(jié)單元之間的特征長度，與剛體有限元中的定義相同[13]，即Cij= hi+ hj。考慮切向方向t = unit (d - (d·n)n)，則切向應(yīng)變可寫作：

圖 3 相鄰單元在交界面上的粘結(jié)作用Fig. 3 Bonding on interface between adjacent elements

根據(jù)三維條件下的彈性矩陣，交界面上兩單元間的彈性應(yīng)力σe可寫作：

式中：E為材料彈性模量；ν為材料泊松比；ε為考慮粘結(jié)節(jié)點(diǎn)之間的法向和切向應(yīng)變，即。

式中：A為交界面面積；n為一個(gè)粘結(jié)面上的粘結(jié)節(jié)點(diǎn)個(gè)數(shù)。將粘結(jié)力轉(zhuǎn)換到單元質(zhì)心上，結(jié)合接觸力等其他受力，即可計(jì)算每個(gè)單元質(zhì)心受到的合力。

針對(duì)每個(gè)粘結(jié)點(diǎn)受到的作用力，可采用不同的失效模型控制粘結(jié)失效。這里采用混合的斷裂能和判斷粘結(jié)失效，即同時(shí)考慮法向拉伸和剪切計(jì)算臨界強(qiáng)度σn和σs，在粘結(jié)力超過臨界應(yīng)力強(qiáng)度后考慮剛度的軟化過程。結(jié)合拉伸和剪切破壞并采用混合的斷裂能模型計(jì)算斷裂能，從而根據(jù)斷裂能計(jì)算粘結(jié)點(diǎn)間的臨界變形[14]。若粘結(jié)點(diǎn)間的變形超過該臨界變形，則該粘結(jié)失效。本文發(fā)展的粘結(jié)失效準(zhǔn)則更全面地考慮了脆性材料的破碎特性，可對(duì)海冰的破壞過程進(jìn)行更精確模擬。

1.3 擴(kuò)展多面體的接觸模型

擴(kuò)展多面體間的接觸可采用基于近似包絡(luò)函數(shù)的接觸搜索算法計(jì)算接觸點(diǎn)和接觸重疊量[11]。這里采用非線性的接觸力模型簡化計(jì)算擴(kuò)展多面體之間的接觸力[15]。其中，法向接觸力Fn包括法向彈性力和法向黏滯力，即：

式中：kn為法向接觸剛度，且；E*和R*分別為等效彈性模量等效單元半徑[16]；為兩單元的法向接觸重疊量；為法向相對(duì)速度；cn為法向阻尼系數(shù)，且有，mAB為等效質(zhì)量[16]，ζn為無量綱阻尼系數(shù)，e為回彈系數(shù)[17]。

切向接觸力與切向重疊量和摩擦力有關(guān)，一般切向力滿足庫倫摩擦定律。這里采用簡化的方式避免與摩擦力的頻繁比較和判斷。其中，切向彈性力寫作：

式中：ct為切向阻尼系數(shù)，這里取[18]；為切向相對(duì)速度。該切向力計(jì)算方法可避免在每個(gè)時(shí)間步中與摩擦力進(jìn)行對(duì)比，提高了計(jì)算效率。

2 破冰船輔助航行中冰荷載的離散元模擬及驗(yàn)證

采用擴(kuò)展多面體離散元方法模擬船舶結(jié)構(gòu)在平整冰區(qū)的航行過程，并通過船舶結(jié)構(gòu)冰荷載的Lindqvist公式和Riska公式驗(yàn)證模擬冰荷載的準(zhǔn)確性。通過模擬破冰船輔助航行時(shí)貨船的冰荷載特點(diǎn)，分析不同參數(shù)條件下冰荷載的變化趨勢。

2.1 離散元模擬參數(shù)

圖4(a)為“雪龍?zhí)枴比壤?，船體結(jié)構(gòu)的長度為167.0 m，型寬22.6 m，吃水9.0 m。將船體結(jié)構(gòu)表面劃分為三角片元用于海冰單元與船體結(jié)構(gòu)之間的接觸搜索，模擬中船體結(jié)構(gòu)視為不變形的剛體結(jié)構(gòu)。平整冰區(qū)的3個(gè)方向采用固定邊界約束，船舶以速度vs在平整冰區(qū)航行。

圖 4 雪龍?zhí)柶票c平整冰作用的離散元模擬Fig. 4 Simulation sketch of the interaction between icebreaker Xuelong and level ice

擴(kuò)展多面體離散元方法可合理分析海冰與結(jié)構(gòu)相互作用中的彎曲破壞模式，在與ISO199006的冰力規(guī)范公式對(duì)比中也體現(xiàn)出了良好的準(zhǔn)確度。模擬中主要的計(jì)算參數(shù)列于表1[19]，考慮海冰受到水的拖曳作用，拖曳力和拖曳力矩可寫作[20]：

表 1 船舶結(jié)構(gòu)破冰過程模擬的主要計(jì)算參數(shù)Tab. 1 Main parameters in the simulation of icebreaking process of ship hull

2.2 船舶結(jié)構(gòu)冰荷載模擬結(jié)果及驗(yàn)證

采用擴(kuò)展多面體離散元方法模擬單個(gè)船體在冰區(qū)航行過程，將穩(wěn)定段冰載荷時(shí)程的均值作為冰阻力，并通過經(jīng)驗(yàn)公式驗(yàn)證分析本文數(shù)值方法的準(zhǔn)確性。在模擬中，平整冰區(qū)的尺寸L×W為600 m × 150 m。

圖5是船舶結(jié)構(gòu)上3個(gè)方向的冰荷載時(shí)程曲線?？梢?，船體剛進(jìn)入冰區(qū)時(shí)x和z方向上的冰荷載經(jīng)過一段時(shí)間的上升，之后冰荷載趨于穩(wěn)定且在一定水平上持續(xù)波動(dòng)。由于船體結(jié)構(gòu)破冰過程中海冰會(huì)經(jīng)歷由船艏至船舯型線不斷收縮變小的過程，海冰會(huì)在船肩處劇烈地由兩側(cè)擠壓船體，因此，y方向上船體結(jié)構(gòu)也會(huì)受到較大的冰荷載作用。這里主要考慮船體是否具有向前的破冰能力，即重點(diǎn)考慮x方向上的冰荷載。在后面的分析中，冰荷載和冰阻力均指x方向的結(jié)構(gòu)受力。

圖 5 3個(gè)方向的船體冰荷載時(shí)程曲線Fig. 5 Time history of ice loads on ship hull in three directions

采用穩(wěn)定階段的冰載荷均值作為冰阻力，將冰阻力與船舶結(jié)構(gòu)冰阻力的Lindqvist公式和Riska公式計(jì)算結(jié)果對(duì)比。Lindqvist公式將船舶結(jié)構(gòu)冰阻力分為船艏處海冰擠壓冰力、海冰彎曲破壞冰力和破碎海冰的浸沒阻力三部分，在冰區(qū)船舶設(shè)計(jì)中使用較為廣泛[4]。Lindqvist中每部分冰力均采用理想化的模型總結(jié)計(jì)算公式，因此其計(jì)算的冰阻力較為保守。Riska公式是在Lindqvist基礎(chǔ)上總結(jié)了波羅的海實(shí)船測量數(shù)據(jù)而得出的經(jīng)驗(yàn)公式，在冰區(qū)船舶的冰荷載計(jì)算中具有代表性[21]。圖5(a)中的虛線為對(duì)應(yīng)工況下Lindqvist公式計(jì)算的冰阻力。圖6是采用擴(kuò)展多面體離散元方法計(jì)算的冰阻力隨冰厚的變化，及其與Lindqvist公式和Riska公式對(duì)比結(jié)果。從圖6可見，三者計(jì)算的冰阻力均隨冰厚增大而增大，且具有明顯的線性增加趨勢；在具體數(shù)值上，離散元結(jié)果基本介于Lindqvist公式和Riska公式之間。因此，通過Lindqvist公式和Riska公式可以充分說明本文擴(kuò)展多面體離散元方法計(jì)算船舶結(jié)構(gòu)冰荷載具有良好的準(zhǔn)確性和可靠性。

圖 6 船體冰荷載的離散元計(jì)算結(jié)果與經(jīng)驗(yàn)公式對(duì)比Fig. 6 Comparison between DEM result and empirical formulas of ice loads on ship hull

圖 7 貨船結(jié)構(gòu)Fig. 7 Hull structure of cargo ship

2.3 破冰船引航下船舶結(jié)構(gòu)冰荷載的離散元分析

采用擴(kuò)展多面體離散元方法模擬破冰船輔助航行中破冰船和商船的冰荷載，即破冰船在前、貨船緊隨其后。平整冰區(qū)的尺寸為1 000 m × 120 m。貨輪長度為200 m，型寬27.8 m，吃水12.3 m（圖7）。

圖8(a)是無引航條件下破冰船和貨船的冰載荷時(shí)程曲線，虛線是穩(wěn)定階段冰載荷的均值。從圖8(a)可以看出，貨船的冰荷載則與破冰船差別較大，進(jìn)入冰區(qū)后冰荷載會(huì)出現(xiàn)類似脈沖形態(tài)的波動(dòng)，但是其整體水平較破冰船要小。由于貨船較破冰船要長，其冰荷載依然比破冰船小很多。圖8(b)是無引航條件下貨船的冰荷載時(shí)程曲線，虛線是穩(wěn)定階段的冰阻力。無引航條件下貨船冰阻力為5.9 MN，比引航條件下的貨船冰阻力要大。但由于貨船不具備破冰能力，其艏柱傾角較小且水線處存在較為尖銳的艏部結(jié)構(gòu)，容易使海冰發(fā)生擠壓破碎，這與破冰船艏部大傾角造成海冰彎曲破壞不同。因此，無引航條件下貨船冰荷載存在近100 MN的峰值載荷，對(duì)貨船的結(jié)構(gòu)安全構(gòu)成嚴(yán)重威脅。因此，在有破冰船引航條件下，貨船的冰阻力則可明顯降低，且不會(huì)存在較大的峰值荷載。

圖 8 有無破冰船引航條件下船舶結(jié)構(gòu)冰載荷時(shí)程曲線Fig. 8 Time history of ice loads on ship hulls with and without icebreaker pilotage

采用擴(kuò)展多面體離散元方法分析不同船速和船寬比條件下兩船的冰阻力特點(diǎn)。船寬比是貨船與破冰船寬度的比，即w2/w1。圖9(a)是不同船速條件下兩船的冰阻力對(duì)比，模擬中兩船的寬度比為1.0。從圖9(a)可以看出，若采用專業(yè)破冰船在破冰能力較弱的船舶之前進(jìn)入冰區(qū)進(jìn)行破冰作業(yè)，即使貨船的結(jié)構(gòu)尺寸比破冰船大，其受到的冰荷載阻力也會(huì)小很多，能夠大大緩解貨船在冰區(qū)航行的阻力。圖9(b)是不同的船寬比條件下兩船的冰阻力，模擬中船速為1.5 m/s。從圖9(b)可以看出，破冰船的破冰阻力并沒有受到貨船的影響，基本保持在相同的水平。在貨船寬度明顯小于破冰船時(shí)，即w2/w1= 0.5時(shí)，貨船冰阻力比破冰船小很多。隨著貨船寬度增大，特別是在w2/w1> 1時(shí)，貨船的冰阻力明顯增大。

圖 9 不同工況下破冰船和貨船的冰阻力Fig. 9 Ice resistance on icebreaker and cargo ship under different conditions

3 結(jié) 語

介紹了基于閔可夫斯基和的擴(kuò)展多面體單元方法，并將該方法應(yīng)用于船舶結(jié)構(gòu)與海冰的相互作用模擬，分析船舶結(jié)構(gòu)上的冰荷載。通過離散元計(jì)算的船舶冰荷載與Lindqvist公式和Riska公式對(duì)比驗(yàn)證了離散元方法，得出本文的計(jì)算結(jié)果在數(shù)值和趨勢上具有良好的準(zhǔn)確性。采用本文方法模擬破冰船引航條件下貨船的冰荷載，分析了不同船寬和船速條件下破冰船和貨船的冰阻力特點(diǎn)。計(jì)算結(jié)果表明，破冰船引航可以有效降低貨船的冰荷載，但在不同船寬比情況下，貨船冰荷載有明顯的差異。本文方法在海冰的離散元模擬中引入了非規(guī)則的計(jì)算單元，并在單元間的粘結(jié)失效中考慮了混合斷裂能和剛度軟化過程，從而使海冰的破壞更符合真實(shí)的物理過程。