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

（1.哈爾濱工程大學(xué)船舶工程學(xué)院，哈爾濱 150001；2.中國(guó)艦船研究設(shè)計(jì)中心，上海 201108；3.啟東中遠(yuǎn)海運(yùn)海洋工程有限公司，江蘇啟東 226200）

0 引 言

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

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

1 基礎(chǔ)理論和計(jì)算模型

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

1.1 DEM基礎(chǔ)理論

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

離散相DEM 冰粒子平動(dòng)運(yùn)動(dòng)動(dòng)量守恒方程為

離散相DEM冰粒子旋轉(zhuǎn)運(yùn)動(dòng)角動(dòng)量守恒方程為

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

1.2 重疊網(wǎng)格

在涉及到物體自由運(yùn)動(dòng)的求解計(jì)算中，重疊網(wǎng)格具有較大的優(yōu)勢(shì)。重疊網(wǎng)格是由Steger在上世紀(jì)50 年代提出的，最大的特點(diǎn)是不同的區(qū)域網(wǎng)格交界能夠發(fā)生重疊、嵌套和覆蓋等拓?fù)潢P(guān)系，能夠明顯減小網(wǎng)格生成的難度，且生成的網(wǎng)格具有較高的質(zhì)量，對(duì)復(fù)雜形狀的適應(yīng)性很強(qiáng)。在重疊網(wǎng)格中，不同區(qū)域的相對(duì)位置關(guān)系可以隨時(shí)間變化。在模擬的過(guò)程中，重疊網(wǎng)格之間節(jié)點(diǎn)的數(shù)據(jù)是根據(jù)插值計(jì)算來(lái)進(jìn)行的。重疊網(wǎng)格各個(gè)子區(qū)域之間的相互作用關(guān)系往往是通過(guò)三個(gè)步驟來(lái)實(shí)現(xiàn)的：定義區(qū)域網(wǎng)格、搜索貢獻(xiàn)單元以及網(wǎng)格之間的數(shù)據(jù)插值。

1.3 數(shù)學(xué)模型的建立

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

圖2 計(jì)算域Fig.2 Computational domain

1.4 船體及碎冰模型的建立

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

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

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

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

1.5 帶有重疊網(wǎng)格的計(jì)算域網(wǎng)格劃分

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

圖5 網(wǎng)格劃分圖Fig.5 Distribution of meshes

2 計(jì)算方法的驗(yàn)證

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

圖6 76K冰區(qū)加強(qiáng)散貨船碎冰區(qū)阻力Fig.6 Resistance of ice strengthened 76k bulk carriers in floe ice region

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

圖8 船艏碎冰運(yùn)動(dòng)姿態(tài)圖Fig.8 Floe ice motion posture at the stem area

3 計(jì)算結(jié)果及分析

3.1 數(shù)值求解方法

求解方法如下：公式（3）為船舶水動(dòng)力計(jì)算公式，其中Yv、Nv為線速度導(dǎo)數(shù)即位置導(dǎo)數(shù)，Yv?、Nv?為線加速度導(dǎo)數(shù)。對(duì)公式（3）進(jìn)行無(wú)量綱化得公式（4）。

將船舶純橫蕩運(yùn)動(dòng)條件代入公式（4），得

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

只需要對(duì)各橫向力和力矩曲線進(jìn)行非線性曲線擬合，得到y(tǒng)=asinωt+bcosωt的形式，即可得出系數(shù)A1、A2、B1、B2，從而得到船舶的水動(dòng)力導(dǎo)數(shù)。

3.2 敞水工況橫蕩運(yùn)動(dòng)水動(dòng)力導(dǎo)數(shù)數(shù)值計(jì)算結(jié)果

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

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

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

表2 各頻率下非線性擬合系數(shù)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船模做純橫蕩運(yùn)動(dòng)時(shí)水動(dòng)力導(dǎo)數(shù)計(jì)算結(jié)果Tab.3 Hydrodynamic derivative calculation results of ship model performing pure sway motion at V=0.8 m/s

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

3.3 碎冰工況冰水耦合后水動(dòng)力導(dǎo)數(shù)數(shù)值計(jì)算結(jié)果

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

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

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

表4 各頻率下非線性擬合系數(shù)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船模做純橫蕩運(yùn)動(dòng)時(shí)水動(dòng)力導(dǎo)數(shù)計(jì)算結(jié)果Tab.5 Hydrodynamic derivative calculation results of ship model performing pure sway motion at V=0.8 m/s

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

4 結(jié) 論

本文使用STAR-CCM+軟件進(jìn)行了碎冰區(qū)船舶橫蕩運(yùn)動(dòng)數(shù)值模擬，求解出敞水以及冰水耦合下船舶橫蕩運(yùn)動(dòng)時(shí)的水動(dòng)力導(dǎo)數(shù)，所得結(jié)論如下：

（1）通過(guò)靜水阻力、碎冰阻力以及碎冰與船舶艏尾運(yùn)動(dòng)姿態(tài)與實(shí)驗(yàn)結(jié)果對(duì)比，驗(yàn)證了本文數(shù)值方法的有效性。

（2）由于冰塊與船發(fā)生碰撞作用時(shí)冰塊的運(yùn)動(dòng)姿態(tài)和軌跡發(fā)生變化，造成船舶的側(cè)向力以及轉(zhuǎn)艏力矩在冰水耦合的情況下都呈現(xiàn)了上下波動(dòng)的情況。

（3）碎冰工況下考慮冰水耦合船舶水動(dòng)力導(dǎo)數(shù)Y'v、N'v、比敞水工況值偏大，而Y'v?比敞水工況值偏小。

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