浮冰區(qū)船舶冰阻力數(shù)值計(jì)算研究

何賓峰 詹成勝 趙橋生 軒慎宇 國 威

(高性能船舶技術(shù)教育部重點(diǎn)實(shí)驗(yàn)室(武漢理工大學(xué))1) 武漢 430063)(武漢理工大學(xué)船海與能源動(dòng)力工程學(xué)院2) 武漢 430063) (中國船舶科學(xué)研究中心3) 無錫 214082)

0 引 言

全球氣候變暖造成極地冰層逐年融化，海冰覆蓋面積持續(xù)減小，影響北極航道通航[1-2].浮冰海況是極地船舶航行將遭遇的一種典型冰況，即便是在夏季通航期，北極航道多處航段都存在冰況不等的浮冰.船舶在浮冰海域航行過程中，易發(fā)生船速降低、舵效不佳、操縱困難、船舶的回轉(zhuǎn)圈大幅增加等現(xiàn)象.對(duì)于浮冰海域航行的船舶，船體冰阻力的大小不僅取決于船舶的船型和航速，還與浮冰的密集度、厚度、尺寸大小等參數(shù)有關(guān).

Cundall[3]提出了一種處理非連續(xù)介質(zhì)問題的數(shù)值模擬方法，即離散元方法(DEM).離散元法對(duì)浮冰冰塊進(jìn)行幾何及物理狀態(tài)建模，能形象地反應(yīng)應(yīng)力場(chǎng)、位移及速度的變化，非常適合于船舶在浮冰區(qū)域航行時(shí)的運(yùn)動(dòng)及冰載荷求解問題.Hopkins[4]將離散元方法引入船舶與破碎冰的相互作用問題，對(duì)浮冰區(qū)的船、冰作用力進(jìn)行分析計(jì)算.Konno等[5]運(yùn)用顆粒離散元方法，用二維圓盤和三維圓球建立浮冰模型，求解離散浮冰中航行船舶受到的冰載荷.Lau等[6]運(yùn)用由塊體離散元方法發(fā)展而來的商用軟件DECICE，計(jì)算了浮冰區(qū)船舶和浮冰的作用力，驗(yàn)證了離散元方法在處理浮冰問題上的可行性.國內(nèi)學(xué)者對(duì)浮冰區(qū)船舶冰阻力的研究起步較晚，季順迎等[7]采用三維圓盤離散元方法，對(duì)碎冰區(qū)航行船舶的冰阻力性能進(jìn)行了數(shù)值計(jì)算研究，得到了不同浮冰參數(shù)條件下的船體冰阻力.王超等[8-9]運(yùn)用離散元模型結(jié)合歐拉多相流方法，研究了船體和碎冰之間的相互作用，對(duì)冰船碰撞時(shí)的運(yùn)動(dòng)響應(yīng)進(jìn)行了探討.

離散元方法很好地解決了浮冰的不連續(xù)介質(zhì)特性，能更好地計(jì)算船體和浮冰之間的碰撞，是現(xiàn)階段國際上研究浮冰區(qū)船舶冰阻力問題所采用的主要方法.但國內(nèi)目前在對(duì)浮冰區(qū)船舶冰阻力的研究中對(duì)冰塊的建模停留在單個(gè)球形顆?；蛘邌蝹€(gè)三維圓盤的結(jié)構(gòu)，研究變量較少，且很少考慮流體作用對(duì)冰阻力的影響.在船冰發(fā)生相互作用時(shí)，冰塊的運(yùn)動(dòng)將對(duì)船體流場(chǎng)有較大的影響，同時(shí)流體的流動(dòng)也會(huì)對(duì)冰塊的運(yùn)動(dòng)產(chǎn)生浮托和阻曳作用，忽略流體作用將會(huì)對(duì)數(shù)值結(jié)果產(chǎn)生較大的誤差.船模試驗(yàn)是一種較為準(zhǔn)確的冰阻力預(yù)報(bào)方法，可以最大程度重現(xiàn)船舶在冰區(qū)航行時(shí)的狀態(tài)，但由于成本高、準(zhǔn)備時(shí)間長，并不適用于初步設(shè)計(jì)階段.文中運(yùn)用EDEM離散元軟件，用結(jié)構(gòu)相同的球形顆粒對(duì)浮冰塊進(jìn)行建模，考慮流體作用對(duì)冰阻力的影響，對(duì)船模試驗(yàn)中的不同航速、密集度、厚度和尺寸大小條件下的直航運(yùn)動(dòng)進(jìn)行模擬，分析不同航速和不同浮冰參數(shù)條件下的冰阻力特性.

1 數(shù)值方法

1.1 基本方程

船在較小尺寸的浮冰區(qū)航行時(shí)，船冰間、冰塊間的相互作用主要體現(xiàn)在碰撞、擠壓等行為，冰塊主要發(fā)生升沉、翻轉(zhuǎn)運(yùn)動(dòng)，冰塊的破碎現(xiàn)象較少，所以本文將船體和冰塊都視為剛體進(jìn)行數(shù)值模擬.水的流動(dòng)會(huì)對(duì)冰塊產(chǎn)生浮托和阻曳作用，忽略流體作用會(huì)對(duì)數(shù)值結(jié)果產(chǎn)生誤差，在離散元EDEM軟件中，添加水對(duì)冰塊的浮力和曳力作用.風(fēng)、浪對(duì)冰塊的影響較小，在數(shù)值計(jì)算中忽略自由液面的影響.

以一個(gè)浮冰冰塊為例，對(duì)離散元方法的基本方程進(jìn)行分析.根據(jù)牛頓第二定律，在每一個(gè)時(shí)間步內(nèi)，建立冰塊的運(yùn)動(dòng)方程：

(1)

(2)

式中：mi為冰塊質(zhì)量；ui為冰塊形心的運(yùn)動(dòng)速度；Fi為對(duì)冰塊的作用力；Ii為冰塊對(duì)形心的轉(zhuǎn)動(dòng)慣性；ωi為冰塊的角速度；Mi為冰塊作用力產(chǎn)生的力矩.

1.2 接觸力

1.2.1接觸變量

當(dāng)兩個(gè)顆粒或者顆粒-結(jié)構(gòu)體之間發(fā)生具有一定重疊量的接觸時(shí)，通過接觸變量表征重疊量的大小，見圖1.

圖1 顆粒-顆粒與顆粒-結(jié)構(gòu)體接觸時(shí)的接觸變量

接觸變量在二維時(shí)為長度量綱，三維時(shí)為面積量綱.船體大小遠(yuǎn)大于冰塊顆粒單元大小，當(dāng)船體與冰塊顆粒單元接觸時(shí)，船體視為半徑無窮大的結(jié)構(gòu)體.

1.2.2Hertz-Mindlin接觸模型

對(duì)于一般的干態(tài)、剛性、不發(fā)生破碎等變化的固體冰塊顆粒及船體，采用Hertz-Mindlin接觸模型描述其之間的接觸作用力[10].

法向彈性力：

(3)

法向阻尼力：

(4)

顆粒之間的切向力：

Ft=-Stδt

(5)

顆粒之間的切向阻尼力：

(6)

切向力合力的最大值受庫侖摩擦定律的限制，即切向力的合力不超過最大靜摩擦力μsFn，其中:μs為靜摩擦系數(shù).

當(dāng)接觸點(diǎn)的兩個(gè)單元在發(fā)生相對(duì)轉(zhuǎn)動(dòng)時(shí)，滾動(dòng)摩擦力將產(chǎn)生額外的力矩作用：

Ti=-μrFnRiωi

(7)

式中：μr為滾動(dòng)摩擦系數(shù)；ωi為接觸點(diǎn)的相對(duì)旋轉(zhuǎn)角速度.

1.3 曳力和浮力

冰塊所受的曳力計(jì)算式為

Fdrag=∑kRecAρf|vf-vp|(vf-vp)

(8)

式中：Re為流體雷諾數(shù)；A為離散元的網(wǎng)格面積；ρf為流體密度；vf和vp分別為流體與顆粒在離散元的網(wǎng)格處的速度；k為流體曳力系數(shù)，取0.627；c為雷諾數(shù)系數(shù)，取-0.012(k和c的值通過CFD直接數(shù)值模擬方法進(jìn)行標(biāo)定得到).

冰塊所受的曳力產(chǎn)生的力矩為

Tdrag=∑FdragLsinθ

(9)

式中：L為表面網(wǎng)格中心到顆粒質(zhì)心的矢量長度；θ為矢量F與矢量L的夾角.

冰塊所受的浮力計(jì)算式為

Fbuoyancy=∑ρgVi

(10)

式中：ρ為水的密度；Vi為冰塊浸入水中的體積.

冰塊所受的浮力產(chǎn)生的力矩為

Tbuoyancy=∑ρgViLsinθ

(11)

2 計(jì)算前處理

選取“雪龍?zhí)枴逼票Ｗ鳛橛?jì)算模型，船模與實(shí)船的縮尺比為30，船型參數(shù)見表1，船模模型見圖2.數(shù)值模擬過程中，船模進(jìn)行直航運(yùn)動(dòng)，約束船模的升沉及縱搖，忽略重力場(chǎng)的影響.

表1 船體模型參數(shù) 單位：m

圖2 船模模型圖

選取顆粒離散元方法對(duì)船模試驗(yàn)中的浮冰冰塊進(jìn)行建模，根據(jù)浮冰的材料特性采用球面擬合非球形冰塊顆粒外形，相接觸的浮冰周邊采用較小小球建模，不接觸的浮冰內(nèi)部采用較大小球進(jìn)行填充，以球面接觸代替冰塊和冰塊之間、冰塊和船體之間的復(fù)雜曲面接觸.在同一工況條件下的每個(gè)四棱柱冰塊的厚度H和邊長D相同，模型見圖3.數(shù)值模擬中浮冰物理屬性根據(jù)船模試驗(yàn)中的海冰參數(shù)，密度取919 kg/m3，剪切模量取8×108Pa，泊松比取0.25.

圖3 浮冰冰塊模型示意圖

文中計(jì)算域?yàn)殚L方體，見圖4.長16 m、寬4.5 m、高0.6 m，長度方向?yàn)閤方向，寬度方向?yàn)閥方向，高度方向?yàn)閦方向.本次數(shù)值模擬不包含水，但設(shè)置一個(gè)水平面，以便于觀察浮冰塊的升沉翻轉(zhuǎn)運(yùn)動(dòng).y方向采用周期性的邊界條件，保持計(jì)算域中的冰塊數(shù)量不變.計(jì)算網(wǎng)格采用離散元的計(jì)算網(wǎng)格，網(wǎng)格大小為0.03 m，時(shí)間步長設(shè)置5×10-6s，網(wǎng)格數(shù)量、冰塊數(shù)量及模擬總時(shí)間設(shè)置根據(jù)具體的工況為準(zhǔn).

3 計(jì)算結(jié)果分析

3.1 直航船冰數(shù)值模擬現(xiàn)象

圖5為0.563 m/s航速，浮冰密集度60%，冰厚0.026 7 m，邊長0.2 m參數(shù)下的3個(gè)時(shí)間段的航行運(yùn)動(dòng)情況，圖6為該工況下的船冰接觸現(xiàn)象.由圖6可見：在船舶航行的過程中，主要在船體和冰塊、冰塊和冰塊之間發(fā)生接觸碰撞，碰撞接觸的部位主要集中在船艏部和肩部，船體平行中體幾乎沒有冰塊的升沉和翻轉(zhuǎn)，船體后方因船體駛過形成開敞區(qū)域.在船冰接觸中，冰塊主要貼著船體表面在船艏部、肩部進(jìn)行升沉和翻轉(zhuǎn)運(yùn)動(dòng)，然后沿著船體表面向船的兩側(cè)滑動(dòng)，并對(duì)周圍的冰塊產(chǎn)生擠壓碰撞.

圖5 船體航行運(yùn)動(dòng)示意圖

圖6 船冰接觸現(xiàn)象圖

3.2 直航航速影響計(jì)算結(jié)果

對(duì)浮冰密集度60%，厚度0.026 7 m，邊長0.2 m條件下的不同航速0.188、0.282、0.376、0.470、0.563 m/s的船體冰阻力進(jìn)行了計(jì)算，得到的不同航速船體冰阻力曲線與文獻(xiàn)的對(duì)比情況，見圖7.由圖7可知：船體所受冰阻力大小隨著航速的增加而增大，隨著航速的進(jìn)一步增加，冰阻力增大的趨勢(shì)減緩.當(dāng)船速較低時(shí)，冰塊幾乎不發(fā)生升沉翻轉(zhuǎn)現(xiàn)象；隨著航速的增加，冰塊的升沉翻轉(zhuǎn)現(xiàn)象慢慢發(fā)生且程度越來越大；隨著航速的繼續(xù)增加，冰塊受到水的曳力影響明顯而使冰塊遠(yuǎn)離船體運(yùn)動(dòng)，冰塊與船體的碰撞頻率減小，導(dǎo)致冰阻力增大的趨勢(shì)減緩.

圖7 不同航速船體冰阻力

經(jīng)無因次化處理后，本次數(shù)值模擬情況，與郭春雨等[11]在某拖曳水池中進(jìn)行的冰區(qū)航行船舶非凍結(jié)冰模型試驗(yàn)的船體冰阻力在量級(jí)上和變化趨勢(shì)上一致，運(yùn)用Ls-dyna軟件的罰函數(shù)算法和流固耦合算法進(jìn)行的冰區(qū)船舶在碎冰中航行的數(shù)值模擬結(jié)果在量級(jí)上和變化趨勢(shì)一致.

3.3 浮冰密集度影響計(jì)算結(jié)果

對(duì)航速0.563 m/s，厚度0.026 7 m，邊長0.2 m條件下的不同密集度40%、50%、60%、70%下的船體冰阻力進(jìn)行了計(jì)算，得到的不同密集度條件下的船體冰阻力曲線圖見圖8.由圖8可知，冰阻力隨著密集度的增加而增大，從60%密集度增加到70%密集度時(shí)，冰阻力有明顯的增長趨勢(shì).

圖8 不同密集度下的船體冰阻力

在船舶航行過程中，船艏能夠迅速地將冰塊沖散到船兩側(cè).在低密集度情況下，冰塊有足夠的空間散開，冰塊的作用范圍很?。辉诟呙芗惹闆r下，船艏并不能夠迅速地將冰塊沖散到船兩側(cè)，而是頂著堆積的冰塊向前行駛，并在船艏形成堆積現(xiàn)象，見圖9，船體與浮冰的碰撞因?yàn)槊芏鹊脑黾佣觿。瑥亩鴮?dǎo)致冰阻力有明顯的增長趨勢(shì).

圖9 70%密集度下的浮冰堆積圖

3.4 浮冰厚度及尺寸影響計(jì)算結(jié)果

對(duì)航速0.563 m/s，密集度60%條件下的不同邊長0.133 3、0.2、0.266 67 m，和不同冰厚0.016 67、0.026 67、0.04 m下的9種不同參數(shù)的船體冰阻力進(jìn)行了計(jì)算.圖10為計(jì)算得到的船體冰阻力曲線圖，在相同的浮冰邊長工況下，冰阻力都隨著厚度的增加而增大，總體趨勢(shì)呈線性相關(guān).冰厚的增加使冰塊的質(zhì)量增加，在船與冰接觸碰撞后，冰塊對(duì)船體的沖擊動(dòng)量變大，從而引起冰阻力的增加.

圖10 船體冰阻力

在相同的厚度工況下，冰阻力隨尺寸的增加而增加，當(dāng)?shù)竭_(dá)一定尺寸后，冰阻力保持穩(wěn)定.當(dāng)冰塊尺寸較小時(shí)，冰塊的作用范圍較小，冰塊主要受質(zhì)量影響，尺寸越大，質(zhì)量越大，冰阻力越大.隨著冰塊尺寸的增加，冰塊的作用范圍變大，一方面尺寸的增加會(huì)導(dǎo)致冰塊的質(zhì)量的增加，使冰阻力增大；另一方面尺寸會(huì)使單位計(jì)算區(qū)域內(nèi)的浮冰數(shù)量減小，冰塊與船體的接觸頻率減小，使冰阻力減小，最終的結(jié)果使船體冰阻力保持穩(wěn)定.

4 結(jié) 論

1)船體冰阻力隨航速的增加而增大，隨著航速的增加，在航速大于0.470 m/s航速時(shí)，冰阻力的增加趨勢(shì)減緩.

2)隨浮冰密集度的增加而增大，在密集度為70%時(shí)，冰阻力增長明顯.

3)隨浮冰厚度的增加而線性增大.

4)隨冰塊邊長的增加而增大，當(dāng)邊長大于0.2 m時(shí)，船體冰阻力的無明顯變化.

本次計(jì)算主要考慮浮冰冰塊由多個(gè)球形顆粒構(gòu)成和流體對(duì)浮冰的影響兩個(gè)方面，數(shù)值模擬結(jié)果通過文獻(xiàn)對(duì)比，驗(yàn)證了該方法的可行性.但本次計(jì)算中浮冰塊仍為剛體，未考慮浮冰的破碎現(xiàn)象及浮冰的隨機(jī)排列，未來的數(shù)值研究應(yīng)在考慮浮冰的破碎現(xiàn)象及浮冰的隨機(jī)性的基礎(chǔ)上進(jìn)行優(yōu)化.