浮冰區船舶冰阻力數值計算研究

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

(高性能船舶技術教育部重點實驗室(武漢理工大學)1) 武漢 430063)(武漢理工大學船海與能源動力工程學院2) 武漢 430063) (中國船舶科學研究中心3) 無錫 214082)

0 引 言

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

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

離散元方法很好地解決了浮冰的不連續介質特性，能更好地計算船體和浮冰之間的碰撞，是現階段國際上研究浮冰區船舶冰阻力問題所采用的主要方法.但國內目前在對浮冰區船舶冰阻力的研究中對冰塊的建模停留在單個球形顆粒或者單個三維圓盤的結構，研究變量較少，且很少考慮流體作用對冰阻力的影響.在船冰發生相互作用時，冰塊的運動將對船體流場有較大的影響，同時流體的流動也會對冰塊的運動產生浮托和阻曳作用，忽略流體作用將會對數值結果產生較大的誤差.船模試驗是一種較為準確的冰阻力預報方法，可以最大程度重現船舶在冰區航行時的狀態，但由于成本高、準備時間長，并不適用于初步設計階段.文中運用EDEM離散元軟件，用結構相同的球形顆粒對浮冰塊進行建模，考慮流體作用對冰阻力的影響，對船模試驗中的不同航速、密集度、厚度和尺寸大小條件下的直航運動進行模擬，分析不同航速和不同浮冰參數條件下的冰阻力特性.

1 數值方法

1.1 基本方程

船在較小尺寸的浮冰區航行時，船冰間、冰塊間的相互作用主要體現在碰撞、擠壓等行為，冰塊主要發生升沉、翻轉運動，冰塊的破碎現象較少，所以本文將船體和冰塊都視為剛體進行數值模擬.水的流動會對冰塊產生浮托和阻曳作用，忽略流體作用會對數值結果產生誤差，在離散元EDEM軟件中，添加水對冰塊的浮力和曳力作用.風、浪對冰塊的影響較小，在數值計算中忽略自由液面的影響.

以一個浮冰冰塊為例，對離散元方法的基本方程進行分析.根據牛頓第二定律，在每一個時間步內，建立冰塊的運動方程：

(1)

(2)

式中：mi為冰塊質量；ui為冰塊形心的運動速度；Fi為對冰塊的作用力；Ii為冰塊對形心的轉動慣性；ωi為冰塊的角速度；Mi為冰塊作用力產生的力矩.

1.2 接觸力

1.2.1接觸變量

當兩個顆粒或者顆粒-結構體之間發生具有一定重疊量的接觸時，通過接觸變量表征重疊量的大小，見圖1.

圖1 顆粒-顆粒與顆粒-結構體接觸時的接觸變量

接觸變量在二維時為長度量綱，三維時為面積量綱.船體大小遠大于冰塊顆粒單元大小，當船體與冰塊顆粒單元接觸時，船體視為半徑無窮大的結構體.

1.2.2Hertz-Mindlin接觸模型

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

法向彈性力：

(3)

法向阻尼力：

(4)

顆粒之間的切向力：

Ft=-Stδt

(5)

顆粒之間的切向阻尼力：

(6)

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

當接觸點的兩個單元在發生相對轉動時，滾動摩擦力將產生額外的力矩作用：

Ti=-μrFnRiωi

(7)

式中：μr為滾動摩擦系數；ωi為接觸點的相對旋轉角速度.

1.3 曳力和浮力

冰塊所受的曳力計算式為

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

(8)

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

冰塊所受的曳力產生的力矩為

Tdrag=∑FdragLsinθ

(9)

式中：L為表面網格中心到顆粒質心的矢量長度；θ為矢量F與矢量L的夾角.

冰塊所受的浮力計算式為

Fbuoyancy=∑ρgVi

(10)

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

冰塊所受的浮力產生的力矩為

Tbuoyancy=∑ρgViLsinθ

(11)

2 計算前處理

選取“雪龍號”破冰船船模作為計算模型，船模與實船的縮尺比為30，船型參數見表1，船模模型見圖2.數值模擬過程中，船模進行直航運動，約束船模的升沉及縱搖，忽略重力場的影響.

表1 船體模型參數 單位：m

圖2 船模模型圖

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

圖3 浮冰冰塊模型示意圖

文中計算域為長方體，見圖4.長16 m、寬4.5 m、高0.6 m，長度方向為x方向，寬度方向為y方向，高度方向為z方向.本次數值模擬不包含水，但設置一個水平面，以便于觀察浮冰塊的升沉翻轉運動.y方向采用周期性的邊界條件，保持計算域中的冰塊數量不變.計算網格采用離散元的計算網格，網格大小為0.03 m，時間步長設置5×10-6s，網格數量、冰塊數量及模擬總時間設置根據具體的工況為準.

3 計算結果分析

3.1 直航船冰數值模擬現象

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

圖5 船體航行運動示意圖

圖6 船冰接觸現象圖

3.2 直航航速影響計算結果

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

圖7 不同航速船體冰阻力

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

3.3 浮冰密集度影響計算結果

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

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

在船舶航行過程中，船艏能夠迅速地將冰塊沖散到船兩側.在低密集度情況下，冰塊有足夠的空間散開，冰塊的作用范圍很小；在高密集度情況下，船艏并不能夠迅速地將冰塊沖散到船兩側，而是頂著堆積的冰塊向前行駛，并在船艏形成堆積現象，見圖9，船體與浮冰的碰撞因為密度的增加而加劇，從而導致冰阻力有明顯的增長趨勢.

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

3.4 浮冰厚度及尺寸影響計算結果

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

圖10 船體冰阻力

在相同的厚度工況下，冰阻力隨尺寸的增加而增加，當到達一定尺寸后，冰阻力保持穩定.當冰塊尺寸較小時，冰塊的作用范圍較小，冰塊主要受質量影響，尺寸越大，質量越大，冰阻力越大.隨著冰塊尺寸的增加，冰塊的作用范圍變大，一方面尺寸的增加會導致冰塊的質量的增加，使冰阻力增大；另一方面尺寸會使單位計算區域內的浮冰數量減小，冰塊與船體的接觸頻率減小，使冰阻力減小，最終的結果使船體冰阻力保持穩定.

4 結 論

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

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

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

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

本次計算主要考慮浮冰冰塊由多個球形顆粒構成和流體對浮冰的影響兩個方面，數值模擬結果通過文獻對比，驗證了該方法的可行性.但本次計算中浮冰塊仍為剛體，未考慮浮冰的破碎現象及浮冰的隨機排列，未來的數值研究應在考慮浮冰的破碎現象及浮冰的隨機性的基礎上進行優化.