冰區加強型集裝箱船碎冰航道航行阻力數值模擬

張遠雙，齊江輝，鄭亞雄，吳述慶

(1.武漢船舶職業技術學院，湖北 武漢 430050；2.武漢第二船舶設計研究所，湖北 武漢 430064)

0 引 言

冰區加強船舶不需要自主破冰，通常需要在破冰船已經開辟的航道中航行，即在碎冰區域航行。船舶在冰水混合的碎冰區航行時，航行阻力與開敞水域有明顯不同。關于船舶在冰水混合流動中的阻力性能預報，通常有理論計算、經驗公式估算、數值模擬和船模試驗等方法。理論計算及經驗公式估算受船型尺度等限制，通常只針對特定船型且精度較差。

在數值模擬方面，國內外開展的相關研究并不多。目前應用最多的是離散元法（DEM）、有限元法（FEM）以及光滑粒子流體動力學方法（SPH）等。Jungyong Wang 等[1]基于Ls-dyna 軟件模擬了碎冰條件下Terrry Fox 號破冰船的阻力性能，分析了碎冰密集度對阻力性能的影響。Moon-Chan Kim 等[2]對一艘冰區加強型散貨船的碎冰區阻力性能進行了數值模擬，通過與試驗結果的對比驗證了數值結果的準確性，其數值結果與試驗值誤差在20%左右。郭春雨等[3]同樣采用Ls-dyna 軟件計算各種工況下的船舶航行阻力特性，并與試驗值進行了對比。王超等[4]基于離散元模型結合歐拉多相流對碎冰區船舶的冰阻力進行了數值模擬，得到了碎冰阻力隨航速等的變化規律。在自編程計算方面，李紫麟[5]基于離散元理論，采用三維圓盤模擬海冰單元，對不同航速碎冰作用于船體的力進行了數值模擬。J.Tuhkuri 等[6]對DEM 方法應用于冰水船相互作用的研究進行了系統的評述。總的來說，碎冰區域航行船舶的阻力性能數值預報仍需要提高精度。

在船模試驗方面，冰水池船舶試驗通常分為凍結模型冰試驗和非凍結模型冰試驗。隨著冰水池技術的發展，國內外冰池船模試驗技術迅速發展，并取得了較多研究成果。首次破冰船阻力試驗是在1964 年，Corlett 和Snaith[7]以石蠟代替冰，對1 艘小型破冰船進行了阻力試驗。隨后德國、芬蘭、韓國、加拿大、俄羅斯等冰水池相繼在船舶冰區阻力、操縱性等方面開展了大量的模型試驗。在國內，郭春雨等[8–9]依托哈爾濱工程大學常溫拖曳水池采用非凍結模型冰進行了一系列的船-冰相互作用試驗。黃焱等[10–11]在天津大學冰力學實驗室，開展了碎冰和平整冰阻力試驗。

本文采用STARCCM+商用軟件，基于離散元方法，結合拉格朗日方法的顆粒接觸動力學模型，建立船-冰-水的CFD-DEM 耦合數值計算模型，對一艘冰區加強型集裝箱船在碎冰區航行的阻力特性進行研究。參照冰水池試驗結果，選取不同的顆粒形狀進行數值模擬，將數值計算現象與試驗現象進行對比，對不同航速、不同冰厚下的阻力特性進行分析。

1 理論基礎

1.1 數值計算方法

計算模型中，流體滿足不可壓縮及連續性方程，同時忽略各相之間的熱傳遞過程[12]。基于VOF 方法捕捉自由液面，湍流模型為SST-K-omege 模型。為防止出流邊界處的流體反射對船體周圍流場的影響，對出流邊界附近采用Choi 方法進行數值消波[13]。數值離散方法采用離散元方法。

離散元方法在20 世紀70 年代初由Cundalland Strack[14]首次提出。離散元方法將研究對象離散成一系列的單元，單元之間通過節點進行連接。其求解過程一般為：1）將研究對象離散為單元顆粒，選定合理的單元連接形式，給定初始條件；2）根據接觸條件判斷單元之間是否發生接觸，若發生接觸則根據單元間的本構模型，計算單元之間相互作用力，從而得到各單元之間的相對速度和位移，更新單元信息。通過對計算域內所有單元的物理量的求解，可以得到整體的研究對象的變形及運動情況。

1.2 離散項接觸力模型

兩單元之間接觸作用模型如圖1 所示，在離散元中接觸力公式實際是彈簧-阻尼器模型的一種變形。彈簧產生排斥力將顆粒分開，而阻尼器表示粘性阻尼并允許模擬除完全彈性外的碰撞類型。作用在接觸點處的接觸力可以看做是一對彈簧-阻尼器振子。

圖1 接觸力模型Fig.1 Model of contact force

Hertz-Mindlin 無滑移接觸模型是一種非線性彈簧-阻尼器接觸模型的變形[15]。作用在2 個單元A 和B 之間的力可表示為：

其中： Fn為法向力分量； Ft為切向力分量。

法向力分量可以寫成：

切向力分量可以寫成：

其中：Nndamp和 Ntdamp分 別為法向阻尼系數和切向阻尼系數，Cnrest和Ctrest分 別為法向和切向的恢復系數，若Cnrest=0 則 Nndamp=1，若 Ctrest=0 則 Nndamp=1； Req， Meq和Geq分別對應等效半徑、等效楊氏模量和等效剪切模量。

對于單元-壁面碰撞，上述公式保持不變，將壁面半徑和質量設定為 Rwall=∞ 和 Mwall=∞，因此等效半徑為 Req=Rparticle， 等效質量 Mwall=Mparticle。

2 數值計算流程

2.1 離散元模型

離散項冰粒子的建模參照冰水池試驗的碎冰圖像，碎冰粒子表現為一系列的圓球組成的復合粒子，碎冰粒子形狀選取為金字塔形和不規則棱柱形，如圖2所示。冰的密度為917 kg/m3。冰的彈性模量為E=9 GPa，根據冰屬性的縮尺比關系[16]，冰粒子模型的彈性模量為300 MPa，泊松比為0.3。

圖2 碎冰幾何形狀Fig.2 Geipetric shape of crushed ice

圖3 數值航道中的碎冰粒子Fig.3 Ice particles in numerical channels

根據北極地區碎冰尺寸統計數據，碎冰尺寸大致服從對數正態分布規律，其分布形式為：

本文數值模擬中，顆粒（2 種碎冰模型）、流體（水）和壁面（船體壁面、航道邊界壁面）多相之間會產生碰撞等相互作用，在計算中考慮7 種相互作用類型，如圖4 所示。其中，船-冰之間的摩擦系數為0.1，冰-冰之間的摩擦系數為0.1，動摩擦系數為0.3。冰粒子采用噴射器的形式進入流場，噴射器定義冰粒子進入的速度、位置及方式等。

2.2 計算模型設置

本文計算模型為1 艘典型冰區加強型集裝箱船，模型縮尺比為30，實船及模型主尺度參數如表1 所示。

圖4 多相相互作用示意圖Fig.4 Schematic diagram of mutiphase interaction

表1 船體主尺度參數Tab.1 The main parameters of ship

碎冰航道計算域以船體為中心，船體前后距計算域邊界距離均為8 L，船體兩側距計算域邊界距離均為1.5 Lpp，在垂直方向上計算域尺寸為?1.5 Lpp～2.0 Lpp。船體、平整冰及計算域的4 個壁面均為無滑移壁面邊界條件，速度入口及壓力出口邊界如圖5 所示。本文計算域分為空氣域和水域，計算過程中船體保持不動，通過改變速度入口處的速度調節船舶的航速。為準確捕捉自由液面位置、船體周圍流場等，在相應的區域進行網格加密處理，邊界層網格形式為棱柱形邊界層，邊界層采用y+壁面處理，y+值范圍為30～60。船體表面網格尺度為5 mm，船體表面第1 層網格厚度為2 mm，計算域中等網格總數約為212 萬，圖6 為網格劃分結果。

圖5 計算域及邊界條件Fig.5 Calculation domain and boundary conditions

圖6 網格劃分結果Fig.6 Result of meshing

3 數值計算結果及分析

3.1 船冰接觸驗證

圖7 為漢堡冰水池進行的某冰區加強型散貨船的碎冰阻力試驗圖像[17]。船舶在碎冰區域航行時，碎冰在與船舶發生碰撞時會發生翻轉現象，而隨著船舶的航行，在船體首部區域會發生碎冰堆積現象，同時一部分碎冰運動到船體底部或舷側隨船體滑移。在這個過程中，碎冰之間以及碎冰與船體之間會發生碰撞。從圖7 可以看出數值計算現象與試驗現象吻合度較好，證明本文采用的碎冰粒子模型及數值方法可以有效模擬碎冰區船舶與碎冰接觸時的現象。

圖7 船-冰碰撞現象模擬Fig.7 Simulation of ship-ice collision

3.2 碎冰運動狀態分析

船舶在碎冰區域航行時，碎冰與船體之間的碰撞、摩擦等作用產生的阻力是船體阻力的重要成分。分析船體周圍碎冰的運動狀態對阻力特性的預報十分重要。計算V=0.563 m/s 時船體周圍碎冰分布如圖8 所示。可以看出碎冰在船體首部發生碰撞后堆積，碎冰與船體間的相對速度驟降。而在船體尾部，由于船體尾流的作用使得碎冰分布較為分散，聚集在螺旋槳周圍的碎冰較少，這可以有效減少冰-槳接觸，有利于螺旋槳的推進效率。從圖9 也可以看出，隨著航速的增加，船體周圍的興波已經較為明顯，興波可以有效地將船體周圍的碎冰排開。同時可以看到，在船體尾部形成了一條浮冰很少的航道，該航道寬度大致相當于船體的寬度，這會有效降低船體-冰之間的接觸力。

3.3 航速對碎冰阻力的影響

圖8 船體周圍碎冰的運動狀態Fig.8 Movement position of crushed ice around hull

圖9 不同航速碎冰狀態對比Fig.9 Movement position of crushed ice in different speed

船舶在敞水區域航行時，航速對船體阻力的影響很大。計算航速分別為0.187 m/s，0.376 m/s，0.563 m/s，0.751 m/s 和0.939 m/s（分別對應實船航速2 kn，4 kn，6 kn，8 kn 和10 kn）時，船體受到的冰阻力時歷曲線如圖10 所示。時歷曲線中阻力分為3 部分，contact force 為船-冰接觸力，resistance 為水阻力，net force 為船體受到的總阻力。可以看出船-冰接觸力隨著航速的增大而增大，同時接觸力有強烈的震蕩特性。接觸力的3 個分量中，縱向力以某一基準值為下限震蕩；側向力則在0 值附近上下震蕩，這是由于船體左右兩舷對稱的原因；而垂向力基本為0，這表明船-冰碰撞基本為水平方向碰撞。低速航行時，縱向接觸力震蕩下限遠大于水阻力；航速較高時，縱向接觸力的震蕩下限小于水阻力。這說明隨著航速的增大，水阻力在船體總阻力中所占的比重越來越高。

對船體阻力的各個分量取穩定段時間平均即得到各航速下的船體阻力分量如圖11 所示。可以看出，船體總阻力及各分量均隨航速的增大而增大，這與敞水航行阻力的規律類似。在航速從0.187 m/s 增大到0.939 m/s 時，航速變為原來的5 倍，總阻力增大為原來的3.5 倍，縱向接觸力增大為原來的2.2 倍，水阻力增大為原來的20.4 倍，說明隨著航速的增大，水阻力的增大速度遠高于船-冰接觸力的增大速度。但在航速增大的過程中，縱向接觸力占總阻力的比例越來越小，航速為0.939 m/s 時比重約為59%；水阻力占總阻力的比例越來越大，航速為0.939 m/s 時比重約為41%。

圖10 低速和高速時船體阻力分量時歷曲線Fig.10 Time history curve of hull resistance at low speed and high speed

圖11 船體總阻力及分量隨航速變化曲線Fig.11 Comparision of resistance components

圖12 不同航速船體表面接觸力分布Fig.12 Contace force distribution of hull surface in different speed

為進一步解釋航速增大而船-冰接觸力增大不明顯的現象，圖12 給出2 種航速下船體表面與碎冰接觸力的分布。可以看出，在低航速時，碎冰在船體周圍聚集使得船體與碎冰的接觸面積較大；而高航速時，高航速及船體的興波等使得船體排開碎冰的效果更明顯，船體周圍聚集的碎冰數量明顯減少，因此船體與碎冰的接觸面積變小，因此航速增大后接觸力的增長不明顯。

4 結 語

本文以1 艘典型冰區加強型集裝箱船為研究目標，建立DEM 和CFD 相結合的數值模擬方法對碎冰區航行船舶的船-冰-水耦合作用進行研究。討論了船體周圍碎冰的運動狀態、航行阻力特性分析以及航速對船體航行阻力的影響等，結論如下：

1）基于CFD-DEM 結合的數值方法可以有效模擬船舶在碎冰區航行的過程，對碎冰在船體附近的堆積、翻轉等運動狀態可以有效模擬，采用該方法進行船-冰-水耦合作用分析是有效可行的。

2）碎冰在與船體接觸作用后會在首部形成一個明顯的減速區，使得碎冰在首部堆積。首部的碎冰一部分沿著船體側邊或底部滑行。在尾部區域，由于尾流的作用會在船體尾部形成一條碎冰數量較少的航道，而隨著航速的增大，該航道碎冰變得更少，同時船體左右兩側的興波使得船體排開碎冰的作用更加明顯。

3）航速對碎冰區航行船舶阻力有重要影響，碎冰區航行阻力與敞水航行阻力曲線特性一致。隨著航速的增大，水阻力的增大速度遠高于船-冰接觸力的增大速度，水阻力占總阻力的比重隨著航速的增大而增大明顯。