船舶碎冰區橫蕩運動的水動力導數研究

王 超，劉 正，朱江波，汪春輝，郭春雨

（1.哈爾濱工程大學船舶工程學院，哈爾濱 150001；2.中國艦船研究設計中心，上海 201108；3.啟東中遠海運海洋工程有限公司，江蘇啟東 226200）

0 引 言

船舶水動力導數表示船舶以某一速度運動，保證其他運動參數不變的情況下，只改變某一個運動參數所引起的船體所受水動力的改變與運動參數的比值。船舶水動力導數與合適的水動力模型相互配合能夠很好地預測船舶在航行時所受到的阻力以及力矩，避免了大量的實船以及船模實驗測量工作，從而節省了大量的人力以及物力。水動力導數研究具有實際性，通常針對某一類型船舶進行研究，不同類型的船舶測量得到的水動力導數相差較大，水動力導數的研究具有復雜性以及針對性。水動力導數的研究方法主要存在經驗公式法、約束船模試驗法和數值模擬方法等。其中，約束船模試驗是目前公認的求解水動力導數較為準確的方法。海軍工程大學吳興亞等[1]針對民用打撈型船采用CFD 軟件STAR-CCM 模擬了敞水中船舶純橫蕩運動，將其結果與理論公式進行對比，驗證了CFD 方法求解水動力導數的可行性，說明了使用該軟件進行計算存在一定的依據。上海交通大學樓鵬宇[2]采用FLUNT軟件針對KVLCC1船型進行了限制水域的純橫蕩運動數值模擬，忽略了自由液面的影響，進行了不同的水域寬度以及深度下船舶水動力導數影響研究，為限制航道船舶操縱性提供了安全指導。Miller 等[3]使用了CFDSHIP-IOWA 程序對船舶斜航、橫蕩、艏搖等數值模擬的RANS 方程進行了求解，得到了部分水動力導數。Turnock 等[4]對船模KVLCC2 進行了橫蕩運動的流場以及水動力數值模擬，將計算結果與實驗進行了對比，較小的誤差驗證了數值手段的合理性。Stern 等[5]在國際會議SIMMAN 2008 的基礎上總結了目前數值模擬船舶操縱性的CFD 手段以及提出采用DES 方法時更細網格可以給出更為精確的結果。Cura-Hochbaum 等[6]對KVLCC1 進行了基于RANS 方法的PMM 運動數值模擬并將求解的水動力導數代入操縱預報方程進行了運動軌跡預報對比，驗證了RANS 方法的可靠性。Kim 等[7]采用基于DSGS 模型的求解器SHIP_Motion，研究了KCS 船模的純橫蕩以及純艏搖試驗，數值計算水動力導數結果與實驗結果吻合較好，說明CFD 手段可以替代實驗研究水動力導數。Yoon等[8-9]基于水面艦艇DTMB 5415，采用自主開發程序CFDShip-Iwoa 進行了PMM 試驗數值模擬，將數值預報得到的相關的水動力導數與試驗結果進行對比分析，大部分線性水動力導數數值預報精度較高，非線性水動力導數預報誤差較大。由此可以看出，非線性水動力導數預報仍然存在較大的難度，由于計算時非線性水動力導數引起的阻力部分不容易被捕捉，進行數值計算時要著重考慮該部分的計算手段。Kim等[10]以及Vroegrigk[11]分別應用有限元軟件LS-DYNA 以及STAR-CCM+軟件，模擬了船舶在碎冰航道中航行時相關性能，并針對冰阻力進行了研究，但是上述兩位學者在研究冰船相互接觸時，只進行了單項耦合設置并未涉及雙向耦合，也沒有對水動力進行相關研究。Lau等[12]使用離散元手段針對海冰與船舶作用下破碎模型以及船舶阻力進行了研究，證明了離散元方法在研究海冰與船舶等結構物的碰撞是可行的，也為后續采用DEM 方法模擬海冰提供了一定的思路以及方法。楊勇[13]等通過求解非定常RANS方程對KVLCC1裸船體在深水以及淺水中進行了純橫蕩運動數值模擬，計算了相關的力和力矩，忽略了自由液面影響，得出了船舶在淺水中所受到的水動力以及線性水動力導數遠大于深水中相應數值的結論，證實了船舶的淺水效應。劉晨飛等[14]運用CFD 以及重疊網格技術模擬了斜航和純橫蕩運動，分別計算了月池封閉以及打開時的船舶水動力導數，得到了水動力導數會因為月池存在而有所增加的結論，證明了使用重疊網格在計算船舶水動力導數時具有一定的準確性，可用于控制船舶運動進行水動力導數求解。哈爾濱工程大學郭春雨等[15]使用LS-DYNA 軟件模擬了船舶碎冰阻力以及船-冰作用過程，在較低速度條件下得出數值模擬結果在定性上與試驗結果比較符合,數值計算方法研究船舶與碎冰相互作用具有一定的可行性。季順迎等[16]運用改進型DEM模型模擬冰的運動變化過程，并將此方法應用于船-冰作用及海洋結構物-冰作用上，針對密集度、流速、冰塊尺寸以及厚度對冰阻力的影響進行了研究，得出了海冰參數的增加均會導致海冰對船體的沖擊動量以及頻率增加的結論，盡管該研究在浮冰形態、船體性能以及流體動力方面進行了諸多假設，但是計算方法以及思路具有較大的參考意義。涂勛程[17]使用軟件LS-DYNA 在浮冰環境和平整環境下對物探船的冰阻力預報及相關參數敏感性開展了數值模擬研究。王超等[18]運用離散元模型結合歐拉多相流，對船體與碎冰之間的相互作用關系進行了探索，得出以下結論：船體所受冰阻力主要是由碎冰與船體表面的摩擦和碰撞產生，并隨航速的增大而增大，但當航速增大到一定值后，碎冰阻力不再增加，甚至還有減小的趨勢。

碎冰區冰水耦合下船舶水動力導數的研究較少，目前研究人員考慮冰水耦合主要針對于碎冰存在水浮力時碎冰阻力以及獨立的水阻力研究，并未考慮耦合后對流體阻力本身的影響。本文擬采用STAR-CCM+軟件中DEM模塊下的雙向耦合模式對碎冰區中某型船舶進行不同運動頻率下的純橫蕩運動模擬，求取船舶所受的側向力以及轉艏力矩，進而經無因次化求解出冰水耦合下的部分水動力導數，為今后考慮耦合下流體增阻后的水動力導數以及操縱性能預報方法研究提供理論基礎和技術支撐。

1 基礎理論和計算模型

本研究采用STAR-CCM+軟件。湍流模型選擇標準k-ε模型，依靠STAR-CCM+中DEM 模塊下的雙向耦合功能達到冰水耦合作用，該耦合模式考慮了離散相對連續相的影響，主要包括位移、相間動量、質量以及熱傳遞的影響。由拉格朗日離散相方程對網格單元進行積分可得出每個顆粒的動量、質量以及能量的變化。通過對穿透體積的所有顆粒進行求和，從而得到網格單元內與連續相交換的總的動量、質量和能量，將其代入連續相方程從而得到冰塊與水的相互耦合作用。使用STAR-CCM+軟件中重疊網格使船模具有橫向運動y=asin（ωt），從而模擬船模橫蕩運動，監測耦合后船舶所受的水動力。

1.1 DEM基礎理論

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

離散相DEM 冰粒子平動運動動量守恒方程為

離散相DEM冰粒子旋轉運動角動量守恒方程為

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

1.2 重疊網格

在涉及到物體自由運動的求解計算中，重疊網格具有較大的優勢。重疊網格是由Steger在上世紀50 年代提出的，最大的特點是不同的區域網格交界能夠發生重疊、嵌套和覆蓋等拓撲關系，能夠明顯減小網格生成的難度，且生成的網格具有較高的質量，對復雜形狀的適應性很強。在重疊網格中，不同區域的相對位置關系可以隨時間變化。在模擬的過程中，重疊網格之間節點的數據是根據插值計算來進行的。重疊網格各個子區域之間的相互作用關系往往是通過三個步驟來實現的：定義區域網格、搜索貢獻單元以及網格之間的數據插值。

1.3 數學模型的建立

單一來說，冰粒子主要屬于固體力學研究范疇，流體力部分屬于流體力學研究范疇，運用STARCCM+中DEM 模塊模擬碎冰粒子，DEM 參數的設置具有開放性，可進行變參數研究。橫蕩運動存在兩部分流體域，其中運動域如圖1所示，該運動域用于控制船體做y=asin（ωt）運動。計算域如圖2所示，為模擬無限水域減少池壁效應，計算域速度入口距離船首部3L，壓力出口距離船尾部3.5L，左右兩側距離船中3.3L，頂部以及底部分別距離船舶2L、3.5L。

圖2 計算域Fig.2 Computational domain

1.4 船體及碎冰模型的建立

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

圖3 碎冰粒子圖Fig.3 Brash ice particle

圖4 探針布置圖Fig.4 Probe layout

表1 某冰區船型參數及碎冰物理參數Tab.1 Ship type parameters of ice region and physical parameters of brash ice

1.5 帶有重疊網格的計算域網格劃分

網格劃分模型選用切割體網格、棱柱層網格以及表面重構。總體網格基礎尺寸為2.3%L，船體表面網格大小選取為基礎尺寸的12.5%，船舶尾部進行網格加密，尺寸為基礎尺寸的6.25%，棱柱層總厚度為0.025 m，棱柱層數為6，棱柱層延伸率為1.2。船體網格如圖5（a）所示，尾部網格加密如圖5（b）所示。由于DEM 粒子的求解是否收斂很大程度上取決于水線面x-y-z向網格密度，分別對自由液面進行x-y-z向加密，z向采用兩層加密，網格尺度分別為基礎尺寸的25%、12.5%。x-y方向網格尺度均為基礎尺寸的80%。計算域整體網格圖如圖5（c）所示，z向加密如圖5（d）所示。考慮重疊網格自身的屬性，將運動域外表面網格尺寸定為基礎尺寸的40%，同時計算域采用兩層過渡，計算域最后一層網格尺寸也為基礎尺寸的40%，達到運動域表面的重疊網格要盡量與背景網格貼合。運動域與背景網格貼合如圖5（e）、（f）所示。重疊網格迭代方式為線性迭代。

圖5 網格劃分圖Fig.5 Distribution of meshes

2 計算方法的驗證

數值計算方法的驗證依托哈爾濱工程大學船模拖曳水池中進行的76K 冰區加強散貨船[19]的碎冰試驗。76K 船模網格以及DEM 相關設置與本文破冰船數值模擬設置一致，通過驗證76K 數值模擬的網格以及DEM 設置的合理性從而驗證本次數值模擬方法的準確性。本次驗證選擇一種碎冰密集度60%、一個航速點0.613 m/s，進行敞水船模阻力、碎冰阻力以及碎冰姿態驗證。實驗采用的碎冰尺寸具有隨機性，存在多種大小不同的碎冰，數值模擬時如果采用多種碎冰容易造成冰水耦合求解發散，極大地消耗計算資源，并且不同的碎冰尺寸所需要的網格設置不同，采用多種碎冰需要在一次計算時設置不同的計算網格。文獻[16]曾說明海冰密集度以及厚度保持不變的情況下，冰塊尺寸的變化主要改變冰塊與船體的碰撞頻率，平均阻力變化不明顯。基于上述原因，本次數值模擬中碎冰厚度采用與實驗中碎冰相同厚度，數值模擬碎冰長寬尺寸選擇與實驗中碎冰尺寸比較小的一類相一致。圖6（a）中靜水阻力數值模擬值取后5 s 均值為5.32 N，靜水阻力實驗值為5.19 N，誤差為3%；圖6（b）中碎冰阻力圖中數值模擬碎冰阻力值均值為4.89 N，實驗值為4.33 N，誤差為13%。使用重疊網格時，背景網格與重疊網格之間進行網格信息的傳輸時不可避免地會發生些許信息的丟失，所以阻力部分的誤差將會偏大，但誤差仍處于可接受范圍之內。從圖7 中可以看出，76K 冰區加強散貨船以及某型破冰船的數值仿真結果與76 K 試驗結果在船舶尾部均出現了碎冰沿著船兩側移動最后在船尾形成的尾跡。從圖8中可看出，76K冰區加強散貨船以及某型破冰船的數值仿真結果與76K試驗結果在船舶首部均出現了碎冰堆積現象，并且沿著船首向船兩側移動。綜上所述，本次數值模擬具有一定的準確性，可以較為真實地模擬船舶在碎冰區航行時相關的阻力性能以及冰塊運動姿態。

圖6 76K冰區加強散貨船碎冰區阻力Fig.6 Resistance of ice strengthened 76k bulk carriers in floe ice region

圖7 尾部區域碎冰分布圖Fig.7 Floe ice distribution at the stern area

圖8 船艏碎冰運動姿態圖Fig.8 Floe ice motion posture at the stem area

3 計算結果及分析

3.1 數值求解方法

求解方法如下：公式（3）為船舶水動力計算公式，其中Yv、Nv為線速度導數即位置導數，Yv?、Nv?為線加速度導數。對公式（3）進行無量綱化得公式（4）。

將船舶純橫蕩運動條件代入公式（4），得

式中，A1、A2、B1、B2分別為

只需要對各橫向力和力矩曲線進行非線性曲線擬合，得到y=asinωt+bcosωt的形式，即可得出系數A1、A2、B1、B2，從而得到船舶的水動力導數。

3.2 敞水工況橫蕩運動水動力導數數值計算結果

純橫蕩運動選擇頻率f分別為0.06 Hz、0.08 Hz、0.1 Hz、0.12 Hz和0.14 Hz，考慮冰區船舶在航行時運動幅度變化不大，按照40縮尺后船模的運動振幅為a=0.075 m，采用重疊網格方法建立運動域，將運動y=asin(wt)賦予船體，同時設置來流速度V=0.8 m/s。分別監測各個頻率下船舶所受的側向力Fy以及轉艏力矩Nz，如圖9和圖10所示。通過ORIGIN軟件進行各個頻率下的橫向力和力矩曲線擬合并將結果進行無量綱化，擬合系數列于表2 中。根據各非線性系數計算結果即可求得各頻率下船模的橫蕩水動力導數，列于表3中。

圖9 敞水工況橫蕩0.06～0.14 Hz船舶所受的側向力Fig.9 Lateral forces on ships in open water conditions at 0.06-0.14 Hz

圖10 敞水工況橫蕩0.06～0.14 Hz船舶所受的轉艏力矩Fig.10 Turning moments on ships in open water conditions at 0.06-0.14 Hz

表2 各頻率下非線性擬合系數A1、B1、A2、B2值Tab.2 Non-linear fitting coefficients of A1,B1,A2 and B2 at different frequencies

表3 V=0.8 m/s船模做純橫蕩運動時水動力導數計算結果Tab.3 Hydrodynamic derivative calculation results of ship model performing pure sway motion at V=0.8 m/s

從圖9、圖10 中可以看出，敞水工況船舶的轉艏力矩以及側向力都隨著頻率的增加而增大，符合理論推導。井上模型[20-21]統計公式求得水動力導數Y'v=0.016 5，N'v=0.005 8。由表3 可知，敞水工況下數值模擬的船舶水動力導數與統計公式相差不大，其中水動力導數Y'v平均誤差較統計公式偏差小，水動力導數N'v平均誤差大。這是由于當船模處于橫蕩運動模式時，改變頻率對船舶所受側向力的變化較轉艏力矩的變化更為敏感，從而導致水動力導數Y'v預報更接近統計值。并且統計公式采用的是井上模型，井上模型是根據10 艘各種類型的船模（油輪3 艘，貨船3 艘，集裝箱船、液化天然氣船、滾裝船、汽車輪渡各1 艘）進行相關的PMM 運動試驗整理出的水動力導數計算公式，本次計算模型采用的是冰區某型破冰船，其型線以及艙室布置、尾部形狀都進行了改進，與傳統船型存在一定的差別，所以一定意義上導致了水動力導數N'v平均誤差大。

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

碎冰工況下船舶運動參數設置與敞水工況一致，即頻率f分別為0.06 Hz、0.08 Hz、0.1 Hz、0.12 Hz和0.14 Hz，運動振幅為a=0.075 m。考慮計算耗時問題，本次模擬碎冰密集度較低為60%，經調研某型破冰船在浮冰工況下（浮冰尺寸大于50～60 cm），連續破冰航速是連續自破冰下航速不超過10 kn，0.2～0.3 m的浮冰基本不受影響，考慮本次模擬情況為低密集度小型碎冰并且不存在破冰船引航，選取船舶速度為9.8 kn，依據傅汝德數相似換算來流速度為0.8 m/s。同時將船體表面設置為壁面，物理模型選擇拉格朗日多相模型中的多相耦合模型。多相相互作用設置冰-冰接觸、冰-船接觸、冰-水接觸、冰-空氣接觸、水-空氣接觸。通過冰-船耦合作用改變冰塊運動姿態進而冰塊對船體周圍流場產生干擾從而達到船-冰-水相互耦合。滾動阻力模型選用Hertz Mindlin，阻力系數模型使用Schiller-Naumann模型。各個頻率下船舶所受的側向力Fy以及轉艏力矩Nz如圖11～12所示。對于各頻率下的橫向力和力矩擬合后得到的系數列于表4中，求得的各頻率下船模的純橫蕩水動力導數值列于表5中。

圖11 碎冰工況橫蕩0.06～0.14 Hz船舶所受的側向力Fig.11 Lateral force on ships in floe ice area conditions at 0.06-0.14 Hz

圖12 碎冰工況橫蕩0.06～0.14 Hz船舶所受的轉首力矩Fig.12 Turning moment on ships in floe ice area conditions at 0.06-0.14 Hz

表4 各頻率下非線性擬合系數A1、B1、A2、B2值Tab.4 Non-linear fitting coefficients of A1,B1,A2 and B2 at different frequencies

表5 V=0.8 m/s船模做純橫蕩運動時水動力導數計算結果Tab.5 Hydrodynamic derivative calculation results of ship model performing pure sway motion at V=0.8 m/s

從圖11～12 中可以看出，冰水耦合下船舶的側向力以及轉艏力矩都呈現上下波動趨勢。這是由于碎冰與船舶碰撞時，運動姿態發生改變，部分會堆積在船舶周圍，部分會沿著船體表面滑行、翻轉，還有極少數的冰塊會貼在船體隨船舶一起運動。這些冰塊不確定的運動擾動了船體周圍流場，從而對船舶的水動力產生上下波動的情況。表5中船舶水動力導數Y'v、N'v比敞水工況值偏大，由此可以判斷冰水耦合下由Y'v、N'v所引起的船舶流體阻力將會增加。耦合后Y'v、N'v水動力導數值較敞水工況值偏大，主要原因在于船舶在碎冰中航行時，冰塊與船體相互作用具有隨機性導致碎冰運動姿態同樣具有隨機性，隨機運動的碎冰對船體周圍的流場擾動不規則，擾動了船體周圍的邊界層，增大了船舶航行時所受的阻力。表5中低頻率時船舶的水動力導數Y'v?、N'v?明顯大于敞水工況下水動力導數值，當頻率大于0.1 Hz時冰水耦合后水動力導數值增加不明顯，主要原因在于隨著頻率的增加，船舶周圍的流場變化劇烈，壓力差明顯，冰塊更容易貼在船體表面，貼在船體表面的碎冰使得水流與船體接觸變少，導致了冰水耦合后在較高頻率時水動力波動減弱，形成了冰水耦合時低頻率水動力導數值要明顯大于敞水工況值的現象。但是Y'v?比敞水工況值偏小，這是由于船舶在航行過程中冰塊對船體的碰撞力是一個瞬時阻礙船舶運動的力，所以冰碎耦合后船舶的加速度較敞水工況值偏小，加速度部分引起船舶的阻力占總的水阻力分量將會減小，從而導致了Y'v?比敞水工況值偏小。水動力導數N'v?與敞水工況值總體上相差不大，這是由于船舶橫向加速度在船舶轉艏力矩中靈敏度偏低，船舶的轉艏力矩主要由N'v引起。隨著運動頻率的增加，部分碎冰會貼在船體兩側，將船體與周圍的流場隔離開，水動力波動減弱，更不易捕捉船舶橫向加速度在船舶轉艏力矩中的變化，所以此項水動力導數值與敞水值偏差不大。

4 結 論

本文使用STAR-CCM+軟件進行了碎冰區船舶橫蕩運動數值模擬，求解出敞水以及冰水耦合下船舶橫蕩運動時的水動力導數，所得結論如下：

（1）通過靜水阻力、碎冰阻力以及碎冰與船舶艏尾運動姿態與實驗結果對比，驗證了本文數值方法的有效性。

（2）由于冰塊與船發生碰撞作用時冰塊的運動姿態和軌跡發生變化，造成船舶的側向力以及轉艏力矩在冰水耦合的情況下都呈現了上下波動的情況。

（3）碎冰工況下考慮冰水耦合船舶水動力導數Y'v、N'v、比敞水工況值偏大，而Y'v?比敞水工況值偏小。

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