孔 帥，季順迎

(大連理工大學(xué) 工業(yè)裝備結(jié)構(gòu)分析國(guó)家重點(diǎn)實(shí)驗(yàn)室，遼寧 大連 116023)

海冰廣泛存在于地球的南北極及高寒地區(qū)，是影響北極航道開(kāi)發(fā)和通航的主要環(huán)境要素[1]。船舶在極地航行中會(huì)受到來(lái)自于海冰和海浪等多方面的動(dòng)態(tài)激勵(lì)，其中海冰對(duì)船舶振動(dòng)的影響最大[2]；另外船冰之間的碰撞，如冰山撞擊，會(huì)對(duì)船舶結(jié)構(gòu)安全產(chǎn)生極大的安全隱患[3]。

通過(guò)數(shù)值方法提前確定冰區(qū)航行的冰載荷可為船體重點(diǎn)區(qū)域安全維護(hù)提供技術(shù)支持，同時(shí)為船舶冰區(qū)航行提供預(yù)警標(biāo)準(zhǔn)以及相應(yīng)安全的航行操作指導(dǎo)。在模擬海冰與船體作用時(shí)，根據(jù)研究需要已開(kāi)發(fā)出眾多的仿真技術(shù)，如采用LS-Dyna動(dòng)力學(xué)軟件平整冰與船體之間的作用過(guò)程[4]、根據(jù)彈性梁理論推導(dǎo)出的海冰彎曲斷裂模式模型[5]、依據(jù)現(xiàn)場(chǎng)觀測(cè)數(shù)據(jù)建立的二維平面內(nèi)環(huán)形裂紋破壞模型[6]等不同數(shù)值方法。

離散元方法由Cundall等[7]提出，最早是應(yīng)用于巖石、谷物等散體材料的動(dòng)力特性分析方法。由于浮冰大量存在于海冰初始生成及融化階段，因此可采用離散元方法對(duì)浮冰進(jìn)行模擬分析，如李紫麟等[8]將蓮葉冰離散成圓盤單元，用于分析船體結(jié)構(gòu)冰與蓮葉冰之間的相互作用，Sun等[9]模擬了浮冰、波浪和海流共同對(duì)直立結(jié)構(gòu)導(dǎo)管架平臺(tái)的作用。另外，在與船體結(jié)構(gòu)之間的作用中，海冰的破碎過(guò)程屬于典型的脆性破壞[10]。脆性材料的破碎過(guò)程可采用具有黏結(jié)-失效效應(yīng)的離散元模型進(jìn)行模擬。如模擬海冰與斜坡結(jié)構(gòu)物的多面體單元Hybrid黏結(jié)失效準(zhǔn)則可表征出漸進(jìn)性過(guò)程[11]，但相比多面體單元，球體單元因其高效的接觸搜索判斷及較為準(zhǔn)確地接觸力計(jì)算而更便于應(yīng)用。Potyondy和Cundall[12]基于球形顆粒提出了平行黏結(jié)模型，該模型已被PFC3D軟件采用并廣泛應(yīng)用于巖石等脆性材料的裂紋生成研究中[13-14]。海冰在結(jié)構(gòu)物作用下的破壞過(guò)程也可采用平行黏結(jié)模型，如龍雪等[15]將其應(yīng)用到錐體海洋結(jié)構(gòu)的尺寸參數(shù)對(duì)平臺(tái)抗冰性能影響分析中。但現(xiàn)有模擬海冰破碎的平行黏結(jié)模型僅考慮了黏結(jié)-失效單元整體的力學(xué)行為，缺少對(duì)球體黏結(jié)單元失效漸進(jìn)性過(guò)程的關(guān)注。

船舶結(jié)構(gòu)冰區(qū)航行過(guò)程中的冰阻力特性分析有助于研究船體結(jié)構(gòu)整體運(yùn)動(dòng)和航行推進(jìn)器[16-17]。為更加有效地模擬船舶結(jié)構(gòu)與海冰的作用過(guò)程，采用三維廣義接觸模型對(duì)海冰的破碎過(guò)程進(jìn)行模擬，并依據(jù)海冰強(qiáng)度試驗(yàn)對(duì)離散元中關(guān)鍵細(xì)觀參數(shù)的影響進(jìn)行了研究。最后將其擴(kuò)展應(yīng)用到冰區(qū)船舶冰載荷分析中。

1 基于廣義接觸模型模擬海冰破碎的離散元方法

廣義接觸模型(3D Generalized contact model，簡(jiǎn)稱GCM-3D)是由GCM-2D發(fā)展而來(lái)[18-19]，能夠很好地應(yīng)用于巖石材料破碎分析。GCM-3D的黏結(jié)單元由局部黏結(jié)點(diǎn)組成，根據(jù)局部黏結(jié)點(diǎn)的失效情況確定整個(gè)黏結(jié)單元是否失效，進(jìn)而判斷整體破壞是否發(fā)生。

1.1 黏結(jié)單元的廣義接觸模型

GCM-3D局部黏結(jié)點(diǎn)分布在黏結(jié)單元的中間面上，其典型分布形式如圖1所示，顆粒單元間的圓點(diǎn)即為局部黏結(jié)點(diǎn)。圖1(a)為局部黏結(jié)點(diǎn)在面(t,n)的分布形式，圖1(b)為局部黏結(jié)點(diǎn)在面(t,s)的分布形式。

圖1 GCM-3D模型的局部黏結(jié)點(diǎn)分布

在顆粒的相互作用過(guò)程中，考慮單元間相對(duì)速度和彈性變形而引起的作用力，顆粒單元之間的接觸可采用彈簧-阻尼器-滑塊的唯象模型[20]，并且單元間的接觸力可解耦為法向分量和切向分量，解耦后法向與切向的接觸模型如圖2所示，單元間的接觸力由各個(gè)局部黏結(jié)點(diǎn)對(duì)之間作用力疊加而成。

圖2 黏結(jié)單元的廣義接觸模型

局部黏結(jié)點(diǎn)對(duì)之間的法向和切向作用力可由線彈性模型計(jì)算，即：

(1)

(2)

(3)

(4)

1.2 黏結(jié)單元局部黏結(jié)點(diǎn)的剛度和強(qiáng)度

GCM-3D模型顆粒單元間剛度由局部黏結(jié)點(diǎn)對(duì)之間剛度疊加，即：

(5)

(6)

局部黏結(jié)點(diǎn)未失效時(shí)法向和切向能承受的最大力由式(7)和式(8)定義，局部黏結(jié)點(diǎn)失效時(shí)的切向力由式(9)定義：

(7)

(8)

Fs=μFn

(9)

綜上可知，局部黏結(jié)點(diǎn)的接觸力本構(gòu)模型采用如圖3所示的本構(gòu)模型。其中，法向最大接觸力由法向黏結(jié)強(qiáng)度決定，剪切力服從Mohr-Coulomb準(zhǔn)則[21]。當(dāng)法向接觸力或切向接觸力達(dá)到式(7)或式(8)的數(shù)值時(shí)，該局部黏結(jié)點(diǎn)對(duì)便失去黏結(jié)效應(yīng)。當(dāng)局部接觸點(diǎn)發(fā)生失效時(shí)，局部接觸點(diǎn)只考慮接觸點(diǎn)之間互相擠壓下的法向接觸力和切向接觸力，其中切向接觸力由單元間的摩擦系數(shù)確定(式(9))。

圖3 局部黏結(jié)點(diǎn)對(duì)間的本構(gòu)模型

2 廣義接觸黏結(jié)模型細(xì)觀參數(shù)對(duì)海冰力學(xué)性質(zhì)的影響

圖4 海冰與船舶結(jié)構(gòu)作用時(shí)的破壞模式

GCM-3D模型中影響海冰宏觀力學(xué)性能的關(guān)鍵細(xì)觀參數(shù)有單元間法向黏結(jié)強(qiáng)度、切向黏結(jié)強(qiáng)度和接觸摩擦系數(shù)等。船舶結(jié)構(gòu)在海冰發(fā)生持續(xù)作用時(shí)，海冰的破壞形式復(fù)雜，由Lindqvist[22]提出船體結(jié)構(gòu)與海冰相互作用時(shí)海冰破壞主要以彎曲破壞為主，局部區(qū)域有擠壓破壞和剪切破壞。圖4為“雪龍?zhí)枴睒O地科學(xué)考察船艏部區(qū)域的海冰破壞圖像，分析圖像資料可知，彎曲斷裂和擠壓破壞為海冰的主要破壞模式。因此文中離散元參數(shù)影響研究依據(jù)海冰彎曲強(qiáng)度試驗(yàn)和單軸壓縮強(qiáng)度試驗(yàn)進(jìn)行。

2.1 海冰單軸壓縮和三點(diǎn)彎曲試驗(yàn)?zāi)M

為進(jìn)行海冰單軸壓縮強(qiáng)度和彎曲強(qiáng)度試驗(yàn)的模擬，采用Ji等[23]在渤海海域所做的海冰強(qiáng)度時(shí)的模型尺寸。海冰的單軸壓縮試驗(yàn)中，采用尺寸為70 mm × 70 mm × 175 mm的長(zhǎng)方體冰試樣，壓頭壓縮速率為2 mm/s，單元數(shù)目為7 800，粒徑大小為5.328×10-3m，并在豎直方向施加載荷；海冰彎曲強(qiáng)度的獲取采用三點(diǎn)彎曲試驗(yàn)方法，試驗(yàn)中海冰為75 mm × 75 mm × 700 mm的長(zhǎng)方體試樣，壓頭壓縮速率為3.5 mm/s，單元數(shù)目為23 940，粒徑大小為6.033×10-3m，并在試樣上表面施加豎直方向載荷。離散元模擬所需的參數(shù)如表1所示。

表1 離散元模擬的主要計(jì)算參數(shù)

海冰試樣在壓頭持續(xù)的作用下會(huì)發(fā)生局部破壞最終導(dǎo)致試樣整體破壞。圖5和圖6為采用離散元模擬的不同加載時(shí)刻海冰在壓頭持續(xù)作用下顆粒間作用力的分布。其中，單軸壓縮試驗(yàn)在未發(fā)生破壞時(shí)顆粒間作用力均勻分布，最終模型整體因發(fā)生剪切錯(cuò)動(dòng)而破壞，與脆性材料力學(xué)的單軸壓縮試驗(yàn)現(xiàn)象一致；三點(diǎn)彎曲試驗(yàn)中在未發(fā)生破壞時(shí)其顆粒間作用力較大區(qū)域集中于壓頭作用位置的上下表面，模型在壓頭持續(xù)作用下最終在壓頭位置處破壞。

圖5 海冰單軸壓縮試驗(yàn)的離散元數(shù)值模擬

圖6 海冰三點(diǎn)彎曲試驗(yàn)的離散元數(shù)值模擬

海冰單軸壓縮和三點(diǎn)彎曲數(shù)值模擬中壓頭施加在加載面上力的變化曲線分別如圖7(a)、(b)所示。將海冰試樣發(fā)生破壞時(shí)的最大壓應(yīng)力用于確定單軸壓縮強(qiáng)度和彎曲強(qiáng)度，其強(qiáng)度分別為4.38 MPa和1.40 MPa。

圖7 單軸壓縮和三點(diǎn)彎曲試驗(yàn)?zāi)M中壓力變化曲線

2.2 局部接觸黏結(jié)點(diǎn)對(duì)間的強(qiáng)度參數(shù)影響

圖和對(duì)海冰單軸壓縮強(qiáng)度的影響

圖和對(duì)海冰三點(diǎn)彎曲強(qiáng)度的影響

圖 和單軸壓縮強(qiáng)度與三點(diǎn)彎曲強(qiáng)度比值的影響

根據(jù)Timco等[24]對(duì)北極海冰強(qiáng)度的統(tǒng)計(jì)分析可知，當(dāng)年冰的單軸壓縮強(qiáng)度可以設(shè)定在0.5～5.0 MPa之間，彎曲強(qiáng)度在1.0 MPa左右。一般情況下，海冰單軸壓縮強(qiáng)度為彎曲強(qiáng)度的2～6倍，數(shù)值模擬中的海冰單軸壓縮強(qiáng)度與彎曲強(qiáng)度之間比值要與實(shí)際情況相吻合[23]。綜合本節(jié)分析可知，北極海冰的單軸壓縮強(qiáng)度、三點(diǎn)彎曲海冰強(qiáng)度可由局部接觸黏結(jié)對(duì)間的強(qiáng)度參數(shù)線性插值獲得；同時(shí)為保證線性插值的穩(wěn)定性，模擬中接觸摩擦系數(shù)確定為μc=0.2。

3 船舶結(jié)構(gòu)的冰阻力特性分析

海冰物理力學(xué)性質(zhì)的復(fù)雜性、船舶結(jié)構(gòu)的復(fù)雜性、海冰-船體之間作用的差異化導(dǎo)致了船舶結(jié)構(gòu)的冰阻力尚未有準(zhǔn)確的計(jì)算公式，目前的冰阻力計(jì)算公式絕大多數(shù)是針對(duì)平整冰區(qū)穩(wěn)定航行的半經(jīng)驗(yàn)公式。這些半經(jīng)驗(yàn)公式結(jié)合了船體破冰過(guò)程的理論分析和現(xiàn)場(chǎng)實(shí)測(cè)數(shù)據(jù)，并將大量相關(guān)的影響參數(shù)引入到公式中，可以較為準(zhǔn)確地反映出不同參數(shù)對(duì)冰阻力的影響，對(duì)破冰船的設(shè)計(jì)選型具有一定的指導(dǎo)意義。但經(jīng)驗(yàn)公式不可對(duì)破冰過(guò)程中的局部冰載荷分布(如線載荷分布、冰壓分布)、冰力時(shí)程特征、冰載荷極值分布等特征進(jìn)行描述。因此可通過(guò)經(jīng)驗(yàn)公式驗(yàn)證數(shù)值方法的有效性，進(jìn)而利用數(shù)值方法全方位描述冰載荷特征的性能。

3.1 Riska平整冰區(qū)船舶冰阻力半經(jīng)驗(yàn)公式

Riska公式是基于室內(nèi)模型試驗(yàn)及波羅的海的大量實(shí)船數(shù)據(jù)分析得出的[25]，并被芬蘭-瑞典冰級(jí)規(guī)范作為冰區(qū)船舶螺旋槳推進(jìn)力推算公式[26]，其公式如下：

F=C1+C2v

(10)

式中：v為船速。該公式考慮了船體形狀因素(如船體進(jìn)水角、船艏長(zhǎng)度和船舯長(zhǎng)度等)、船體航行速度及航行過(guò)程中海冰厚度的影響，C1和C2的取值如公式所示：

(11)

(12)

式中：α為水線處進(jìn)水角，T為冰區(qū)水線，B為型寬，L為水線處全船長(zhǎng)度，Lpar為船舯長(zhǎng)度，Lbow為艏部長(zhǎng)度，hi為冰厚，f1～f4、g1～g3等經(jīng)驗(yàn)系數(shù)由試驗(yàn)測(cè)定，其取值分別為0.23 kN/m3、4.58 kN/m3、1.47 kN/m3、0.29 kN/m3、18.9 kN/(m·s-1·m1.5)、0.67 kN/(m·s-1·m1.5)、1.55 kN/(m·s-1·m1.5)。

3.2 Lindqvist平整冰區(qū)船舶冰阻力半經(jīng)驗(yàn)公式

Lindqvist冰力公式也是目前使用較多的破冰船冰阻力計(jì)算經(jīng)驗(yàn)公式[22]，主要考慮了海冰在彎曲斷裂模式下的冰載荷。涉及到的參數(shù)包括主尺度、船型、冰厚、摩擦以及海冰的彎曲強(qiáng)度。該方法中冰阻力包括破冰阻力、海冰彎曲阻力和浸沒(méi)阻力三部分，其中，破冰阻力Rc為：

(13)

彎曲阻力Rb為：

(14)

式中：E為海冰彈性模量，υ為泊松比，B為船體型寬，ρw為海水密度，g為重力加速度。

浸沒(méi)阻力Rs為：

Rs=(ρw-ρi)ghtotBK

(15)

(16)

式中：ρi為海冰密度，htot為冰與雪的厚度和，T為冰區(qū)水線，L為船體長(zhǎng)度。

考慮速度影響下的總冰阻力Ri：

(17)

3.3 基于廣義接觸模型的平整冰區(qū)船舶冰阻力分析

球體離散元是一種簡(jiǎn)化的計(jì)算模型，無(wú)法完全體現(xiàn)海冰各向異性的特征，如柱狀晶體生長(zhǎng)結(jié)構(gòu)，但可通過(guò)對(duì)如厚度方向、海冰平面內(nèi)方向設(shè)置有差別的黏結(jié)強(qiáng)度模擬海冰強(qiáng)度各向異性的特征。現(xiàn)有的離散元模擬中的黏結(jié)強(qiáng)度雖未考慮各向異性特征，但仍可模擬出結(jié)構(gòu)物在海冰作用下的冰載荷特征，包括時(shí)程曲線、峰值特征及分布規(guī)律[15, 20]。同時(shí)，離散元模擬海冰與結(jié)構(gòu)物作用過(guò)程中，海冰的粒徑越小越能表征出海冰破碎的特征，但會(huì)導(dǎo)致計(jì)算資源大幅度增加，通常選擇在厚度方向布置兩層球體單元即可較為精確地描述出冰載荷特征[15]。

為分析平整冰區(qū)船體冰阻力特性，建立如圖11左側(cè)區(qū)域的平整冰區(qū)，區(qū)域大小為240 m×95 m；圖11右側(cè)為厚度方向上球形顆粒的分布形式，通過(guò)網(wǎng)格搜索后的相鄰顆粒之間初始都采用廣義接觸模型進(jìn)行黏結(jié)。

圖11 由球體顆粒黏結(jié)組成的平整冰區(qū)示意圖

表2為采用離散元模擬時(shí)破冰船在平整冰區(qū)穩(wěn)定航行時(shí)的參數(shù)，其中海冰顆粒之間的法向、切向黏結(jié)強(qiáng)度為0.5 MPa。經(jīng)2.2節(jié)中細(xì)觀參數(shù)-宏觀參數(shù)對(duì)應(yīng)關(guān)系可知，海冰彎曲強(qiáng)度和單軸壓縮強(qiáng)度分別為0.93 MPa和2.79 MPa；同時(shí)，海冰彈性模量為4.41 GPa。以“雪龍?zhí)枴睒O地科學(xué)考察船為分析對(duì)象，船體結(jié)構(gòu)由3 344個(gè)三角形單元組成，并以1.54 m/s的恒定船速在平整冰區(qū)航行。表2也列出了Riska公式[25]和Lindqvist公式[22]中的計(jì)算參數(shù)，其海冰與結(jié)構(gòu)間的摩擦系數(shù)采用了Lindqvist建議的數(shù)值區(qū)間[0.1,0.16]，文中取μ=0.15。

表2 離散元模擬中的主要計(jì)算參數(shù)

圖12展示了破冰船在0.8 m厚的平整冰區(qū)運(yùn)行時(shí)不同角度的破冰狀況。從圖12(a)和(b)可以看出海冰會(huì)在艏部的持續(xù)作用下產(chǎn)生彎曲破壞的斷裂模式，并有海冰在船肩位置處因擠壓的作用發(fā)生擠壓破壞并發(fā)生堆積。因此船舶重點(diǎn)防護(hù)和加固的區(qū)域應(yīng)設(shè)置在船艏及船肩位置處。

圖12 破冰船在平整冰區(qū)航行不同角度的展示

破冰船在破冰過(guò)程中的冰力主要由兩部分組成：一部分來(lái)自船艏部位沖擊破冰時(shí)冰的破碎和彎曲；另一部分來(lái)自于船肩、船舯等部位與海冰之間的摩擦作用，并且摩擦力與有效的接觸長(zhǎng)度有關(guān)，長(zhǎng)度越長(zhǎng)，摩擦作用越明顯。圖13為破冰船在冰區(qū)航行時(shí)的冰阻力變化曲線，可以看出船舶結(jié)構(gòu)的冰載荷因海冰不斷的發(fā)生彎曲斷裂破壞，是一個(gè)典型的具有周期性的載荷；另外隨著船體不斷的駛?cè)氡鶇^(qū)，其冰阻力也在逐步提升。當(dāng)全船進(jìn)入冰后(x=150 m)，其冰阻力也趨于穩(wěn)定。選取x=150～230 m的穩(wěn)定破冰階段冰載荷畫于圖13的右上角(其固定邊界僅在附近5 m范圍內(nèi)即x=235 m才會(huì)有影響，因而固定邊界對(duì)x=150～230 m的冰載荷沒(méi)有影響)，均值為1 367.8 kN，其均值可用于確定冰區(qū)穩(wěn)定航行時(shí)螺旋槳的推進(jìn)力，即：

Fp=Fice+Fwater

(18)

式中：Fp為螺旋槳冰區(qū)穩(wěn)定航行時(shí)所需推進(jìn)力，F(xiàn)ice為冰區(qū)穩(wěn)定航行時(shí)的冰阻力，F(xiàn)water為開(kāi)闊水域航行時(shí)的阻力。

圖13 冰厚為0.8 m時(shí)船舶冰阻力時(shí)程曲線

圖14(a)為冰厚在0.8～1.2 m時(shí)穩(wěn)定破冰階段的冰阻力與經(jīng)驗(yàn)公式確定的冰載荷對(duì)比，圖14(b)為三種船速(1.03 m/s、1.54 m/s和2.06 m/s)下1.2 m厚冰的冰阻力與經(jīng)驗(yàn)公式數(shù)值對(duì)比?？梢钥闯鰯?shù)值仿真的結(jié)果整體上與Lindqvist的結(jié)果更加接近，與Riska的結(jié)果相差較多。Riska公式數(shù)值較小主要是其參考波羅的海的破冰船的測(cè)量數(shù)據(jù)，但波羅的海的冰相對(duì)較??；且該公式對(duì)冰載荷與船速、冰厚之間進(jìn)行了線性簡(jiǎn)化，使其數(shù)值普遍小于Lindqvist冰阻力公式[27]。數(shù)值仿真的結(jié)果從整體上符合船舶結(jié)構(gòu)冰載荷的規(guī)律，為確定冰區(qū)航行螺旋槳推進(jìn)力和冰載荷分布規(guī)律提供了可行的計(jì)算方法。

圖14 數(shù)值仿真中冰載荷均值與經(jīng)驗(yàn)公式的對(duì)比

4 結(jié) 語(yǔ)

采用廣義接觸模型建立了船舶結(jié)構(gòu)冰載荷計(jì)算的離散元方法，同時(shí)也分析了不同冰厚和航速下冰載荷的變化規(guī)律。為分析廣義接觸模型中關(guān)鍵細(xì)觀參數(shù)對(duì)冰載荷模擬的影響，對(duì)海冰單軸壓縮試驗(yàn)和三點(diǎn)彎曲試驗(yàn)進(jìn)行了數(shù)值模擬，確定了強(qiáng)度參數(shù)對(duì)不同海冰強(qiáng)度的影響，并采用線性插值的方式獲得北極海冰的強(qiáng)度參數(shù)。該模型中的強(qiáng)度參數(shù)可以用于根據(jù)現(xiàn)場(chǎng)冰情、船舶航行特征確定相應(yīng)的離散元計(jì)算參數(shù)。通過(guò)對(duì)極地科學(xué)考察船“雪龍?zhí)枴痹谄秸鶇^(qū)航行時(shí)的冰載荷分析可知，該數(shù)值模型可以較好地對(duì)應(yīng)冰載荷經(jīng)驗(yàn)公式。該數(shù)值方法可用于冰區(qū)船舶螺旋槳功率配比及船舶在冰區(qū)航行冰阻力特性研究。為提升該數(shù)值方法的適用性，后續(xù)將會(huì)增加更多冰況、航行方式分析船舶結(jié)構(gòu)的整體冰阻力、局部冰載荷分布特性。