船舶碎冰區斜航運動的水動力導數預報

王 超，劉 正，傅江妍，汪春輝,郭春雨

(1.哈爾濱工程大學 船舶工程學院，哈爾濱 150001; 2.中國艦船研究設計中心，上海 201108)

破冰船在冰區航行時，船體周圍流場受到冰面的干擾，其流場特性必然發生變化，故在求解船舶所受水動力時需額外考慮冰面干擾。水動力導數求解相對于船舶其他方面的研究較少，主要存在經驗公式法、約束船模實驗法、數值模擬方法等，其中約束船模試驗是目前公認的求解水動力導數較為準確的方法。針對敞水水動力導數研究主要有：吳興亞等[1]基于民用打撈型船采用CFD軟件STAR-CCM數值模擬敞水中船舶斜航運動，驗證了CFD方法求解水動力導數的可行性。Ohmori等[2]等采用有限體積法分別對Esso Osake船模作斜航運動的黏性流場和水動力進行了計算。Alessandrini等[3]對系列60船模作定常斜航以及回轉運動時的黏性流場和所受側向力和縱向分布進行了計算。Carrica等[4]采用RANS方程數值模擬船舶較快速度下周圍流場分布并計算了部分水動力，后期又針對回轉以及Z型操縱運動進行了相關的數值模擬。Miller[5]使用了CFDSHIP-IOWA程序對斜航、橫蕩、艏搖等RANS方程進行了求解從而進行了相關數值模擬。Ismail等[6]對裸船體模型在0°以及12°漂角進行了斜航數值模擬，得到了船舶水動力導數的預報采用2階TVD格式更精結論。Cura-Hochbaum[7]對KVLCC1的郵輪進行了PMM運動數值模擬，預報相應的水動力導數，并將其代入操縱模型中與自航試驗對比，驗證了RANS方法的可靠性。Simonsen等[8]對KCS船舶做艏搖運動進行了數值模擬，對所受的力以及力矩的不確定進行了分析。Simonsen等[9]基于求解器CFDShip-Iowa，數值模擬PMM運動試驗，求解部分水動力導數值，數值模擬了KCS船舶的Z型操縱運動以及全回轉運動。

Kim等[10]應用有限元軟件LS-DYNA模擬了碎冰航道阻力，并與冰水池實驗結果進行了對比，驗證了方法的準確性。Vroegrijk[11]利用STAR-CCM+軟件數值模擬了船舶在碎冰航道中的航行時相關性能。Lau等[12]使用離散元手段針對海冰與船舶作用下破碎模型以及船舶阻力進行了研究。劉晨飛等[13]運用CFD以及重疊網格技術模擬了斜航、純橫蕩運動分別計算了月池封閉以及打開時船舶水動力導數，得到了水動力導數會因為月池存在有所增加。孟慶杰[14]采用湍流模型中兩方程K-W SST模型模擬了淺水工況下船體的吸底效應，不同漂角下水動力性能以及船體周圍流場的變化。王建華[15]采用自主開發求解器naoe-FoAM-SJTU數值模擬了船舶斜航運動，得到了不同漂角下的轉艏力矩以及側向以及相應的水動力導數。針對冰區相關的研究主要有：郭春雨等[16]使用LS-DYNA軟件模擬了船舶碎冰阻力以及船-冰作用過程，得出數值模擬結果在較低速度條件下在定性上與試驗結果比較符合。任奕舟等[17]采用有限元方法建模對冰區船舶阻力進行了數值模擬，得到的冰阻力較為準確，適用于預報破冰船冰阻力。齊江輝等[18]使用STAR-CCM+軟件數值模擬了碎冰區的航行過程，分析了航速對航行阻力的影響以及碎冰在船體周圍的運動狀態，得出船舶航行阻力隨著航速的增加有所提高，碎冰阻力由于存在波浪排開作用隨著航速待提高緩慢增加的結論。涂勛程[19]使用軟件LS-Dyna在碎冰環境和平整環境下物探船的冰阻力預報及相關參數敏感性開展了數值模擬研究并與經驗公式對比誤差較小，驗證了耦合方法的準確性。王超等[20]運用離散元模型結合歐拉多相流，對船體與碎冰之間的相互作用關系進行探索，計算不同航速、不同碎冰密集度下船體的受力情況，從直觀上解釋碎冰阻力的變化原因。

碎冰區冰水耦合下船舶水動力導數預報的研究缺失，可查到的文獻考慮冰水耦合也僅僅用于碎冰阻力方面未考慮耦合后對流體阻力本身的影響，大部分只針對于船舶直航時碎冰阻力以及水阻力獨立研究。本文擬采用STAR-CCM+軟件中DEM模塊下雙向耦合模式對碎冰區中某型船舶進行不同漂角下的斜航運動模擬，求取船舶所受的側向力以及轉艏力矩，進而求解出冰水耦合下的部分水動力導數，為今后考慮耦合下流體增阻后的水動力導數以及操縱性能預報方法研究提供理論基礎和技術支撐。

1 基礎理論和計算模型

本文主要采用CFD軟件STAR-CCM+手段獲得，需要滿足連續性方程以及動量守恒方程。選用目前廣泛使用的標準的k-ε模型處理標準的兩方程中湍流問題。依靠STAR-CCM+中DEM模塊下的雙向耦合模式，進行冰塊與水的動量、能量交換達到耦合作用，監測耦合后的船舶所受的水動力。由于雙向耦合基礎理論內容較多請詳見STAR-CCM+軟件幫助文檔，此處不進行詳述。

1.1 DEM基礎理論

碎冰工況航行時存在船-冰-水耦合作用的過程，本文采用CFD軟件STAR-CCM+中DEM模塊模擬碎冰粒子，碎冰粒子由若干DEM基本顆粒組成。拉格朗日構架下，通過顯式計算進行計算域中所有粒子的軌跡跟蹤。CFD-DEM耦合計算包括冰粒子之間的相互作用以及離散相冰與流體相水之間的相互作的計算，離散相DEM冰粒子的運動通常包括平行運動和轉動運動兩種方式，需滿足平動動量守恒以及旋轉角動量守恒方程具體形式為：

(1)

(2)

式中:mi、vi和wi分別為離散冰項粒子編號i的質量、速度和角速度;Ii為編號i的粒子的轉動慣量；Ri為粒子i的半徑；Fg=mig為粒子i的重力；Ffluid為流體對粒子i的作用力；Fij為粒子i與粒子j或壁面之間的碰撞力以及其他作用在粒子上的非接觸力；Tij為接觸力矩，表示除粒子重力以外的接觸力在粒子上產生的扭矩，Tij=r×Fij-μr|Fij|ωi/|ωi|,其中r為粒子重心到粒子接觸點的矢量。

1.2 數學模型的建立

單一來說冰粒子主要屬于固體力學研究范疇，流體力部分屬于流體力學研究范疇，軟件STAR-CCM+中雙向耦合模塊可實現二者耦合。DEM參數的設置具有開放性，可進行變參數研究。圍繞船體建立數值水池，為模擬無限水域，速度入口距離船首部3.0 L，壓力出口距離船尾部3.5 L。水池左右兩側距離船中3.3 L，水池頂部以及底部分別距離船舶2.0、3.5 L，如圖1所示。

圖1 計算域劃分圖

1.3 幾何模型建立

碎冰區冰塊厚度范圍在0.3～3.0 m[19]，小型碎冰跨度在2.0 m左右，大型碎冰尺寸跨度為0.5～10.0 km[19]。考慮模擬采用的碎冰粒子是用DEM復合顆粒填充形成的碎冰粒子，尺寸過大噴射碎冰粒子耗時多，綜合考慮本次碎冰模型采用小型碎冰，長度、寬度、厚度分別為3.0、3.0、1.5 m，按縮尺比40進行建模，碎冰如圖2所示。碎冰的DEM粒子填充數為120。采用STAR-CCM+軟件中探針進行噴射碎冰，其中探針的長度為6.0 m,分辨率為38,如圖3所示，噴射時顆粒流率為每秒噴射9個。某型冰區船模主尺度以及碎冰物理性質見表1，其中表1中碎冰的物理參數指的是依據實尺度冰屬性與縮尺的關系換算得到的模型冰參數。船模幾何圖形如圖4所示。

圖2 碎冰粒子圖

圖3 探針布置圖

表1 某冰區船型參數

圖4 船體模型幾何圖

1.4 網格劃分

斜航工況漂角選擇0°,2°,4°,6°,8°網格劃分方法為各工況網格尺度不變，需要將船體分別旋轉2°,4°,6°,8°將船體與計算域進行布爾減運算，然后劃分網格。模型網格劃分選用切割體網格、棱柱層網格以及表面重構相結合的方式。總體網格基礎尺寸為2.3%L，船體表面網格大小選取為基礎尺寸的12.5%，船舶尾部進行網格加密，尺寸為基礎尺寸的6.25%，為使船體周圍的流場過渡緩和，設置棱柱層數6層，棱柱層延伸率1.2。船體網格圖如圖5(a)所示，尾部加密圖如圖5(b)所示。使用DEM模擬冰粒子時粒子求解是否收斂很大程度上取決于水線面處網格密度，考慮防止自由液面發散，分別對自由面進行x-y-z向加密，z向采用兩層加密，網格尺度分別為基礎尺寸的25.0%、12.5%。x-y方向網格尺度均為基礎尺寸的80%。以斜航0°為例。計算域整體網格圖如圖5(c)所示，z向加密圖如圖5(d)所示。

圖5 網格劃分圖

2 計算方法的驗證

本次數值計算方法的驗證依托哈爾濱工程大學船模拖曳水池中進行的76K冰區加強散貨船[21]的碎冰實驗。76 K船模網格設置以及DEM相關設置與本文某型破冰船數值模擬設置一致，通過驗證76 K數值模擬的網格以及DEM設置的合理性從而驗證本次數值模擬方法的準確性。本次驗證選擇一種碎冰密集度40%、一個航速點0.613 m/s，進行敞水船模阻力、碎冰阻力以及碎冰姿態驗證。進而說明本次計算網格設置的合理性。阻力圖如圖6所示，圖7為冰塊與船艏相互作用的數值模擬與實驗測量的對比圖。圖8為船體航行后尾部區域碎冰分布圖。靜水阻力取圖6(a)后5 s均值為-5.157 N，靜水阻力實驗值為-5.160 N，可看出數值模擬結果與試驗吻合度非常高。碎冰阻力值取圖6(b)25 s后穩定段均值為-4.950 N，實驗值為-4.330 N，誤差控制在15%以內，滿足目前行業內冰阻力預報誤差不超20%的要求。圖7中可看出76 K冰區加強散貨船以及某型破冰船的數值仿真結果與76 K試驗結果在船舶首部均出現了碎冰堆積現象。圖8中可以看出76 K冰區加強散貨船以及某型破冰船的數值仿真結果與76 K試驗結果在船舶尾部均出現了碎冰向船兩側移動形成的尾跡圖。綜上所述本次數值模擬具有一定的準確性，可以較為真實地模擬船舶在碎冰區航行時相關的阻力性能以及冰塊運動姿態。

圖6 76 K冰區加強散貨船碎冰區阻力圖

圖7 船艏碎冰運動姿態圖

圖8 尾部區域碎冰分布圖

3 計算結果及分析

3.1 敞水工況水動力導數數值計算結果

斜航運動選擇漂角為0°,2°,4°,6°,8°，數值模擬采用相對運動模擬船舶斜航，計算來流速度為0.8 m/s，分別計算船舶的側向力以及轉艏力矩如圖9所示，其中船舶的側向力是指沿著船寬方向的力。用最小二乘法分別擬合無因次后的側向力Fy、轉艏力矩Nz，所得曲線在原點處斜率即為水動力導數。側向力Fy、轉艏力矩Nz的計算結果見表2。

圖9 船舶所受側向力及轉艏力矩圖

表2 不同漂角下船舶所受側向力及轉艏力矩

將船舶在不同漂角所受的側向力以及轉艏力矩依據式(3)進行無因次化后擬合得到的曲線如圖10所示。

(3)

圖10 敞水側向力及轉艏力矩擬合圖

斜航敞水工況水動力導數計算結果以及統計公式計算結果(采用井上模型)見表3，可以看出Y′v、N′v相對偏差為20%、5.2%，Y′v誤差偏大的原因是由于處于斜航運動模式時變化漂角時船舶力矩的變化較側向力更為敏感求解結果更為準確，并且統計公式采用的是井上模型，井上模型是根據10艘各種類型的船模(油輪3艘，貨船3艘，集裝箱船、液化天然氣船、滾裝船、汽車輪渡各1艘)進行相關試驗整理出的水動力導數計算公式，本次計算模型采用的是冰區某型破冰船其型線以及艙室布置、尾部形狀都進行了改進與傳統船型存在一定的差別，綜上考慮文章數值計算方法具有一定的有效性。

表3 斜航敞水工況水動力導數計算結果

3.2 碎冰工況冰水耦合后水動力導數數值計算結果

碎冰工況的斜航運動參數設置與敞水工況一致，漂角為0°,2°,4°,6°,8°。船體表面設置為壁面，物理模型中選擇拉格朗日混合模型中多相耦合模塊達到水與冰塊流固耦合作用。多相相互作用設置冰-冰接觸、冰-船接觸、冰-水接觸、冰-空氣接觸和水-空氣接觸。通過冰-船耦合作用改變冰塊運動姿態進而冰塊對船體周圍流場產生干擾從而達到船冰水相互耦合。阻力模型選用Hertz Mindlin，阻力系數模型使用Schiller-Naumann模型，考慮計算耗時問題，本次模擬為碎冰低密集度40%下水動力相關模擬，經調研某型破冰船在浮冰工況下(浮冰尺寸大于50～60 cm)，連續破冰航速是：連續自破冰下，航速不超過18.52 km/h，0.2～0.3 m的浮冰基本不受影響，考慮本次模擬情況為低密集度小型碎冰并且不存在破冰船引航，選取船舶速度為18.149 6 km/h，依據弗汝德數相似換算得來流速度為0.8 m/s。分別計算冰水耦合后船舶的側向力以及轉艏力矩，如圖11～13所示。

圖11 碎冰工況斜航0°～8°船舶所受的側向力

圖12 碎冰工況最小側向力及最大側向力擬合圖

圖13 碎冰工況斜航0°-8°船舶所受的轉艏力矩

3.2.1 冰水耦合下斜航側向力

從圖11中可以看出在每一個漂角工況下，船舶的側向力都呈現了上、下波動的情況，這是由于冰塊與船發生碰撞作用時冰塊的運動姿態發生變化，部分會翻轉以及沿著船表面滑行還有一部分會被波浪排開，這些都對流場造成一定規律性擾動。本次數據處理考慮冰塊的干擾具有隨機性所以采用干擾的最大值以及最小值分別求解對應的水動力導數，從而形成一個水動力波動區間，更為真實地預報冰水耦合下的船舶水動力導數。其中最大以及最小值求取方法采用將15～33 s時間分為12個時間段，第i個時間段水動力波動的峰值為Ai、Bi，則整個時間段水動力波動峰值為

側向力具體值見表4。不同漂角所受的側向力進行無因次化擬合后得到的曲線如圖12所示。碎冰工況冰水耦合水動力導數計算結果見表5。

表4 碎冰工況下不同漂角時船舶所受側向力

表5 碎冰工況下斜航水動力導數計算結果

3.2.2 冰水耦合下斜航轉艏力矩

從圖13中可以看出在每個漂角工況下，船舶轉艏力矩同樣會由于冰塊的運動姿態發生改變(主要包括部分會翻轉、沿著船表面滑行以及一部分會被波浪排開)對船體周圍流場造成干擾從而導致力矩值上下波動的情況。力矩數據處理方法與側向力方法一致。轉艏力矩具體值見表6。將浮冰工況下船舶在不同漂角所受的轉艏力矩進行無因次化后擬合得到的曲線如圖14所示。

表6 碎冰工況下不同漂角時船舶所受轉艏力矩

圖14 碎冰工況最小轉艏力矩以及最大轉艏力矩擬合圖

由表3、表5和表7中水動力導數計算結果可以看出碎冰工況冰水耦合下最大側向力得到的線性水動力導數值F′ymax小于敞水計算結果，最小側向力得到的線性水動力導數值F′ymin大于敞水計算結果。最大轉艏力矩得到的線性水動力導數值N′zmax小于敞水計算結果，最小轉艏力矩得到的線性水動力導數值N′zmin略大于敞水計算結果，線性水動力導數區間端點的絕對值存在一端數值大于敞水計算結果一端數值小于敞水計算結果的現象，并且非線性水動力導數計算結果呈現正負值。這是由于船舶在碎冰中航行時冰塊與船體作用使得冰塊運動姿態發生變化，冰塊會在船體周圍翻轉、在船艏的堆積以及由于壓力差緊貼船體表面的滑行，使得冰面對船體周圍的流場造成了隨機性的干擾，非線性水動力導數計算結果呈現正負值可看出對流場的干擾更具有隨機性，并且非線性水動力導數值在冰水耦合下不容易捕捉求解，波動靈敏性太大。

表7 碎冰工況下斜航水動力導數計算結果

4 結 論

1)本文通過與本次數值模擬具有相同網格以及DEM相關參數設置的76 K冰區散貨船碎冰區的數值模擬，從靜水阻力、碎冰阻力以及碎冰與船舶艏尾運動姿態與實驗進行了對比驗證，并將敞水下船舶進行斜航模擬得到的水動力導數值與井上模型統計公式數值對比，數值較為接近，進一步說明了數值方法的可靠性。

2)船舶的側向力以及轉艏力矩在冰水耦合的情況下都呈現了上、下波動的情況，這是由于冰塊與船發生碰撞作用時冰塊的運動姿態發生變化，部分會翻轉以及沿著船表面滑行還有一部分會被波浪排開，這些都對流場造成一定規律性擾動。

3)水動力導數計算結果可以看出最大側向力以及最大轉艏力矩得到的線性水動力導數值F′ymax、N′zmax小于敞水計算結果，最小側向力以及最小轉艏力矩得到的線性水動力導數值F′ymin、N′zmin大于敞水計算結果，這是由于船舶在碎冰中航行時冰塊與船體作用使得冰塊運動姿態發生變化，冰面對船體周圍的流場造成了隨機性的干擾。隨機性取其區間值能體現冰塊對流場干擾的隨機性。

4)本文使用CFD-DEM方法雖然一定程度上解決了冰水耦合的問題，但還有許多設置需要改進，同時受到時間和計算能力等的限制，本文只計算了一種密集度下冰水耦合的情況，后期將針對于不同密集度以及碎冰的不同尺寸部分進行研究。