張遠雙,齊江輝,鄭亞雄,吳述慶
(1.武漢船舶職業(yè)技術(shù)學(xué)院,湖北 武漢 430050;2.武漢第二船舶設(shè)計研究所,湖北 武漢 430064)
冰區(qū)加強船舶不需要自主破冰,通常需要在破冰船已經(jīng)開辟的航道中航行,即在碎冰區(qū)域航行。船舶在冰水混合的碎冰區(qū)航行時,航行阻力與開敞水域有明顯不同。關(guān)于船舶在冰水混合流動中的阻力性能預(yù)報,通常有理論計算、經(jīng)驗公式估算、數(shù)值模擬和船模試驗等方法。理論計算及經(jīng)驗公式估算受船型尺度等限制,通常只針對特定船型且精度較差。
在數(shù)值模擬方面,國內(nèi)外開展的相關(guān)研究并不多。目前應(yīng)用最多的是離散元法(DEM)、有限元法(FEM)以及光滑粒子流體動力學(xué)方法(SPH)等。Jungyong Wang 等[1]基于Ls-dyna 軟件模擬了碎冰條件下Terrry Fox 號破冰船的阻力性能,分析了碎冰密集度對阻力性能的影響。Moon-Chan Kim 等[2]對一艘冰區(qū)加強型散貨船的碎冰區(qū)阻力性能進行了數(shù)值模擬,通過與試驗結(jié)果的對比驗證了數(shù)值結(jié)果的準確性,其數(shù)值結(jié)果與試驗值誤差在20%左右。郭春雨等[3]同樣采用Ls-dyna 軟件計算各種工況下的船舶航行阻力特性,并與試驗值進行了對比。王超等[4]基于離散元模型結(jié)合歐拉多相流對碎冰區(qū)船舶的冰阻力進行了數(shù)值模擬,得到了碎冰阻力隨航速等的變化規(guī)律。在自編程計算方面,李紫麟[5]基于離散元理論,采用三維圓盤模擬海冰單元,對不同航速碎冰作用于船體的力進行了數(shù)值模擬。J.Tuhkuri 等[6]對DEM 方法應(yīng)用于冰水船相互作用的研究進行了系統(tǒng)的評述??偟膩碚f,碎冰區(qū)域航行船舶的阻力性能數(shù)值預(yù)報仍需要提高精度。
在船模試驗方面,冰水池船舶試驗通常分為凍結(jié)模型冰試驗和非凍結(jié)模型冰試驗。隨著冰水池技術(shù)的發(fā)展,國內(nèi)外冰池船模試驗技術(shù)迅速發(fā)展,并取得了較多研究成果。首次破冰船阻力試驗是在1964 年,Corlett 和Snaith[7]以石蠟代替冰,對1 艘小型破冰船進行了阻力試驗。隨后德國、芬蘭、韓國、加拿大、俄羅斯等冰水池相繼在船舶冰區(qū)阻力、操縱性等方面開展了大量的模型試驗。在國內(nèi),郭春雨等[8–9]依托哈爾濱工程大學(xué)常溫拖曳水池采用非凍結(jié)模型冰進行了一系列的船-冰相互作用試驗。黃焱等[10–11]在天津大學(xué)冰力學(xué)實驗室,開展了碎冰和平整冰阻力試驗。
本文采用STARCCM+商用軟件,基于離散元方法,結(jié)合拉格朗日方法的顆粒接觸動力學(xué)模型,建立船-冰-水的CFD-DEM 耦合數(shù)值計算模型,對一艘冰區(qū)加強型集裝箱船在碎冰區(qū)航行的阻力特性進行研究。參照冰水池試驗結(jié)果,選取不同的顆粒形狀進行數(shù)值模擬,將數(shù)值計算現(xiàn)象與試驗現(xiàn)象進行對比,對不同航速、不同冰厚下的阻力特性進行分析。
計算模型中,流體滿足不可壓縮及連續(xù)性方程,同時忽略各相之間的熱傳遞過程[12]?;赩OF 方法捕捉自由液面,湍流模型為SST-K-omege 模型。為防止出流邊界處的流體反射對船體周圍流場的影響,對出流邊界附近采用Choi 方法進行數(shù)值消波[13]。數(shù)值離散方法采用離散元方法。
離散元方法在20 世紀70 年代初由Cundalland Strack[14]首次提出。離散元方法將研究對象離散成一系列的單元,單元之間通過節(jié)點進行連接。其求解過程一般為:1)將研究對象離散為單元顆粒,選定合理的單元連接形式,給定初始條件;2)根據(jù)接觸條件判斷單元之間是否發(fā)生接觸,若發(fā)生接觸則根據(jù)單元間的本構(gòu)模型,計算單元之間相互作用力,從而得到各單元之間的相對速度和位移,更新單元信息。通過對計算域內(nèi)所有單元的物理量的求解,可以得到整體的研究對象的變形及運動情況。
兩單元之間接觸作用模型如圖1 所示,在離散元中接觸力公式實際是彈簧-阻尼器模型的一種變形。彈簧產(chǎn)生排斥力將顆粒分開,而阻尼器表示粘性阻尼并允許模擬除完全彈性外的碰撞類型。作用在接觸點處的接觸力可以看做是一對彈簧-阻尼器振子。
圖1 接觸力模型Fig.1 Model of contact force
Hertz-Mindlin 無滑移接觸模型是一種非線性彈簧-阻尼器接觸模型的變形[15]。作用在2 個單元A 和B 之間的力可表示為:
其中: Fn為法向力分量; Ft為切向力分量。
法向力分量可以寫成:
切向力分量可以寫成:
其中:Nndamp和 Ntdamp分 別為法向阻尼系數(shù)和切向阻尼系數(shù),Cnrest和Ctrest分 別為法向和切向的恢復(fù)系數(shù),若Cnrest=0 則 Nndamp=1,若 Ctrest=0 則 Nndamp=1; Req, Meq和Geq分別對應(yīng)等效半徑、等效楊氏模量和等效剪切模量。
對于單元-壁面碰撞,上述公式保持不變,將壁面半徑和質(zhì)量設(shè)定為 Rwall=∞ 和 Mwall=∞,因此等效半徑為 Req=Rparticle, 等效質(zhì)量 Mwall=Mparticle。
離散項冰粒子的建模參照冰水池試驗的碎冰圖像,碎冰粒子表現(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)計數(shù)據(jù),碎冰尺寸大致服從對數(shù)正態(tài)分布規(guī)律,其分布形式為:
本文數(shù)值模擬中,顆粒(2 種碎冰模型)、流體(水)和壁面(船體壁面、航道邊界壁面)多相之間會產(chǎn)生碰撞等相互作用,在計算中考慮7 種相互作用類型,如圖4 所示。其中,船-冰之間的摩擦系數(shù)為0.1,冰-冰之間的摩擦系數(shù)為0.1,動摩擦系數(shù)為0.3。冰粒子采用噴射器的形式進入流場,噴射器定義冰粒子進入的速度、位置及方式等。
本文計算模型為1 艘典型冰區(qū)加強型集裝箱船,模型縮尺比為30,實船及模型主尺度參數(shù)如表1 所示。
圖4 多相相互作用示意圖Fig.4 Schematic diagram of mutiphase interaction
表1 船體主尺度參數(shù)Tab.1 The main parameters of ship
碎冰航道計算域以船體為中心,船體前后距計算域邊界距離均為8 L,船體兩側(cè)距計算域邊界距離均為1.5 Lpp,在垂直方向上計算域尺寸為?1.5 Lpp~2.0 Lpp。船體、平整冰及計算域的4 個壁面均為無滑移壁面邊界條件,速度入口及壓力出口邊界如圖5 所示。本文計算域分為空氣域和水域,計算過程中船體保持不動,通過改變速度入口處的速度調(diào)節(jié)船舶的航速。為準確捕捉自由液面位置、船體周圍流場等,在相應(yīng)的區(qū)域進行網(wǎng)格加密處理,邊界層網(wǎng)格形式為棱柱形邊界層,邊界層采用y+壁面處理,y+值范圍為30~60。船體表面網(wǎng)格尺度為5 mm,船體表面第1 層網(wǎng)格厚度為2 mm,計算域中等網(wǎng)格總數(shù)約為212 萬,圖6 為網(wǎng)格劃分結(jié)果。
圖5 計算域及邊界條件Fig.5 Calculation domain and boundary conditions
圖6 網(wǎng)格劃分結(jié)果Fig.6 Result of meshing
圖7 為漢堡冰水池進行的某冰區(qū)加強型散貨船的碎冰阻力試驗圖像[17]。船舶在碎冰區(qū)域航行時,碎冰在與船舶發(fā)生碰撞時會發(fā)生翻轉(zhuǎn)現(xiàn)象,而隨著船舶的航行,在船體首部區(qū)域會發(fā)生碎冰堆積現(xiàn)象,同時一部分碎冰運動到船體底部或舷側(cè)隨船體滑移。在這個過程中,碎冰之間以及碎冰與船體之間會發(fā)生碰撞。從圖7 可以看出數(shù)值計算現(xiàn)象與試驗現(xiàn)象吻合度較好,證明本文采用的碎冰粒子模型及數(shù)值方法可以有效模擬碎冰區(qū)船舶與碎冰接觸時的現(xiàn)象。
圖7 船-冰碰撞現(xiàn)象模擬Fig.7 Simulation of ship-ice collision
船舶在碎冰區(qū)域航行時,碎冰與船體之間的碰撞、摩擦等作用產(chǎn)生的阻力是船體阻力的重要成分。分析船體周圍碎冰的運動狀態(tài)對阻力特性的預(yù)報十分重要。計算V=0.563 m/s 時船體周圍碎冰分布如圖8 所示??梢钥闯鏊楸诖w首部發(fā)生碰撞后堆積,碎冰與船體間的相對速度驟降。而在船體尾部,由于船體尾流的作用使得碎冰分布較為分散,聚集在螺旋槳周圍的碎冰較少,這可以有效減少冰-槳接觸,有利于螺旋槳的推進效率。從圖9 也可以看出,隨著航速的增加,船體周圍的興波已經(jīng)較為明顯,興波可以有效地將船體周圍的碎冰排開。同時可以看到,在船體尾部形成了一條浮冰很少的航道,該航道寬度大致相當于船體的寬度,這會有效降低船體-冰之間的接觸力。
圖8 船體周圍碎冰的運動狀態(tài)Fig.8 Movement position of crushed ice around hull
圖9 不同航速碎冰狀態(tài)對比Fig.9 Movement position of crushed ice in different speed
船舶在敞水區(qū)域航行時,航速對船體阻力的影響很大。計算航速分別為0.187 m/s,0.376 m/s,0.563 m/s,0.751 m/s 和0.939 m/s(分別對應(yīng)實船航速2 kn,4 kn,6 kn,8 kn 和10 kn)時,船體受到的冰阻力時歷曲線如圖10 所示。時歷曲線中阻力分為3 部分,contact force 為船-冰接觸力,resistance 為水阻力,net force 為船體受到的總阻力??梢钥闯龃?冰接觸力隨著航速的增大而增大,同時接觸力有強烈的震蕩特性。接觸力的3 個分量中,縱向力以某一基準值為下限震蕩;側(cè)向力則在0 值附近上下震蕩,這是由于船體左右兩舷對稱的原因;而垂向力基本為0,這表明船-冰碰撞基本為水平方向碰撞。低速航行時,縱向接觸力震蕩下限遠大于水阻力;航速較高時,縱向接觸力的震蕩下限小于水阻力。這說明隨著航速的增大,水阻力在船體總阻力中所占的比重越來越高。
對船體阻力的各個分量取穩(wěn)定段時間平均即得到各航速下的船體阻力分量如圖11 所示??梢钥闯觯w總阻力及各分量均隨航速的增大而增大,這與敞水航行阻力的規(guī)律類似。在航速從0.187 m/s 增大到0.939 m/s 時,航速變?yōu)樵瓉淼? 倍,總阻力增大為原來的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
為進一步解釋航速增大而船-冰接觸力增大不明顯的現(xiàn)象,圖12 給出2 種航速下船體表面與碎冰接觸力的分布。可以看出,在低航速時,碎冰在船體周圍聚集使得船體與碎冰的接觸面積較大;而高航速時,高航速及船體的興波等使得船體排開碎冰的效果更明顯,船體周圍聚集的碎冰數(shù)量明顯減少,因此船體與碎冰的接觸面積變小,因此航速增大后接觸力的增長不明顯。
本文以1 艘典型冰區(qū)加強型集裝箱船為研究目標,建立DEM 和CFD 相結(jié)合的數(shù)值模擬方法對碎冰區(qū)航行船舶的船-冰-水耦合作用進行研究。討論了船體周圍碎冰的運動狀態(tài)、航行阻力特性分析以及航速對船體航行阻力的影響等,結(jié)論如下:
1)基于CFD-DEM 結(jié)合的數(shù)值方法可以有效模擬船舶在碎冰區(qū)航行的過程,對碎冰在船體附近的堆積、翻轉(zhuǎn)等運動狀態(tài)可以有效模擬,采用該方法進行船-冰-水耦合作用分析是有效可行的。
2)碎冰在與船體接觸作用后會在首部形成一個明顯的減速區(qū),使得碎冰在首部堆積。首部的碎冰一部分沿著船體側(cè)邊或底部滑行。在尾部區(qū)域,由于尾流的作用會在船體尾部形成一條碎冰數(shù)量較少的航道,而隨著航速的增大,該航道碎冰變得更少,同時船體左右兩側(cè)的興波使得船體排開碎冰的作用更加明顯。
3)航速對碎冰區(qū)航行船舶阻力有重要影響,碎冰區(qū)航行阻力與敞水航行阻力曲線特性一致。隨著航速的增大,水阻力的增大速度遠高于船-冰接觸力的增大速度,水阻力占總阻力的比重隨著航速的增大而增大明顯。