張遠(yuǎn)雙，齊江輝，鄭亞雄，吳述慶

(1.武漢船舶職業(yè)技術(shù)學(xué)院，湖北 武漢 430050；2.武漢第二船舶設(shè)計(jì)研究所，湖北 武漢 430064)

0 引 言

冰區(qū)加強(qiáng)船舶不需要自主破冰，通常需要在破冰船已經(jīng)開辟的航道中航行，即在碎冰區(qū)域航行。船舶在冰水混合的碎冰區(qū)航行時(shí)，航行阻力與開敞水域有明顯不同。關(guān)于船舶在冰水混合流動(dòng)中的阻力性能預(yù)報(bào)，通常有理論計(jì)算、經(jīng)驗(yàn)公式估算、數(shù)值模擬和船模試驗(yàn)等方法。理論計(jì)算及經(jīng)驗(yàn)公式估算受船型尺度等限制，通常只針對(duì)特定船型且精度較差。

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

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

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

1 理論基礎(chǔ)

1.1 數(shù)值計(jì)算方法

計(jì)算模型中，流體滿足不可壓縮及連續(xù)性方程，同時(shí)忽略各相之間的熱傳遞過程[12]?；赩OF 方法捕捉自由液面，湍流模型為SST-K-omege 模型。為防止出流邊界處的流體反射對(duì)船體周圍流場(chǎng)的影響，對(duì)出流邊界附近采用Choi 方法進(jìn)行數(shù)值消波[13]。數(shù)值離散方法采用離散元方法。

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

1.2 離散項(xiàng)接觸力模型

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

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

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

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

法向力分量可以寫成：

切向力分量可以寫成：

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

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

2 數(shù)值計(jì)算流程

2.1 離散元模型

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

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

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

根據(jù)北極地區(qū)碎冰尺寸統(tǒng)計(jì)數(shù)據(jù)，碎冰尺寸大致服從對(duì)數(shù)正態(tài)分布規(guī)律，其分布形式為：

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

2.2 計(jì)算模型設(shè)置

本文計(jì)算模型為1 艘典型冰區(qū)加強(qiáng)型集裝箱船，模型縮尺比為30，實(shí)船及模型主尺度參數(shù)如表1 所示。

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

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

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

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

圖6 網(wǎng)格劃分結(jié)果Fig.6 Result of meshing

3 數(shù)值計(jì)算結(jié)果及分析

3.1 船冰接觸驗(yàn)證

圖7 為漢堡冰水池進(jìn)行的某冰區(qū)加強(qiáng)型散貨船的碎冰阻力試驗(yàn)圖像[17]。船舶在碎冰區(qū)域航行時(shí)，碎冰在與船舶發(fā)生碰撞時(shí)會(huì)發(fā)生翻轉(zhuǎn)現(xiàn)象，而隨著船舶的航行，在船體首部區(qū)域會(huì)發(fā)生碎冰堆積現(xiàn)象，同時(shí)一部分碎冰運(yùn)動(dòng)到船體底部或舷側(cè)隨船體滑移。在這個(gè)過程中，碎冰之間以及碎冰與船體之間會(huì)發(fā)生碰撞。從圖7 可以看出數(shù)值計(jì)算現(xiàn)象與試驗(yàn)現(xiàn)象吻合度較好，證明本文采用的碎冰粒子模型及數(shù)值方法可以有效模擬碎冰區(qū)船舶與碎冰接觸時(shí)的現(xiàn)象。

圖7 船-冰碰撞現(xiàn)象模擬Fig.7 Simulation of ship-ice collision

3.2 碎冰運(yùn)動(dòng)狀態(tài)分析

船舶在碎冰區(qū)域航行時(shí)，碎冰與船體之間的碰撞、摩擦等作用產(chǎn)生的阻力是船體阻力的重要成分。分析船體周圍碎冰的運(yùn)動(dòng)狀態(tài)對(duì)阻力特性的預(yù)報(bào)十分重要。計(jì)算V=0.563 m/s 時(shí)船體周圍碎冰分布如圖8 所示?？梢钥闯鏊楸诖w首部發(fā)生碰撞后堆積，碎冰與船體間的相對(duì)速度驟降。而在船體尾部，由于船體尾流的作用使得碎冰分布較為分散，聚集在螺旋槳周圍的碎冰較少，這可以有效減少冰-槳接觸，有利于螺旋槳的推進(jìn)效率。從圖9 也可以看出，隨著航速的增加，船體周圍的興波已經(jīng)較為明顯，興波可以有效地將船體周圍的碎冰排開。同時(shí)可以看到，在船體尾部形成了一條浮冰很少的航道，該航道寬度大致相當(dāng)于船體的寬度，這會(huì)有效降低船體-冰之間的接觸力。

3.3 航速對(duì)碎冰阻力的影響

圖8 船體周圍碎冰的運(yùn)動(dòng)狀態(tài)Fig.8 Movement position of crushed ice around hull

圖9 不同航速碎冰狀態(tài)對(duì)比Fig.9 Movement position of crushed ice in different speed

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

對(duì)船體阻力的各個(gè)分量取穩(wěn)定段時(shí)間平均即得到各航速下的船體阻力分量如圖11 所示。可以看出，船體總阻力及各分量均隨航速的增大而增大，這與敞水航行阻力的規(guī)律類似。在航速?gòu)?.187 m/s 增大到0.939 m/s 時(shí)，航速變?yōu)樵瓉淼? 倍，總阻力增大為原來的3.5 倍，縱向接觸力增大為原來的2.2 倍，水阻力增大為原來的20.4 倍，說明隨著航速的增大，水阻力的增大速度遠(yuǎn)高于船-冰接觸力的增大速度。但在航速增大的過程中，縱向接觸力占總阻力的比例越來越小，航速為0.939 m/s 時(shí)比重約為59%；水阻力占總阻力的比例越來越大，航速為0.939 m/s 時(shí)比重約為41%。

圖10 低速和高速時(shí)船體阻力分量時(shí)歷曲線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

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

4 結(jié) 語(yǔ)

本文以1 艘典型冰區(qū)加強(qiáng)型集裝箱船為研究目標(biāo)，建立DEM 和CFD 相結(jié)合的數(shù)值模擬方法對(duì)碎冰區(qū)航行船舶的船-冰-水耦合作用進(jìn)行研究。討論了船體周圍碎冰的運(yùn)動(dòng)狀態(tài)、航行阻力特性分析以及航速對(duì)船體航行阻力的影響等，結(jié)論如下：

1）基于CFD-DEM 結(jié)合的數(shù)值方法可以有效模擬船舶在碎冰區(qū)航行的過程，對(duì)碎冰在船體附近的堆積、翻轉(zhuǎn)等運(yùn)動(dòng)狀態(tài)可以有效模擬，采用該方法進(jìn)行船-冰-水耦合作用分析是有效可行的。

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

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