倪寶玉，黃其，陳綰綬，薛彥卓

哈爾濱工程大學(xué) 船舶工程學(xué)院，黑龍江 哈爾濱 150001

0 引 言

近年來(lái)，北極地區(qū)的航線開(kāi)發(fā)和資源勘探問(wèn)題受到了學(xué)者們的關(guān)注。隨著極地冰的融化，極地航行難度降低，越來(lái)越多的船舶嘗試在極地航行。在此背景下，研究極地航行船舶的操縱性就具有很強(qiáng)的現(xiàn)實(shí)意義。極地航行船舶除了需要承受風(fēng)、波浪、水流等常規(guī)阻力外，還要承受冰阻力，在這些外部阻力中，冰阻力一般最為突出。當(dāng)一艘船航行于層冰冰面時(shí)，船與冰的強(qiáng)烈非線性碰撞會(huì)對(duì)船舶的航行產(chǎn)生很大影響。因此，從保證船舶航行安全來(lái)說(shuō)，能夠準(zhǔn)確預(yù)測(cè)船舶冰阻力至關(guān)重要。

自19世紀(jì)60年代起，針對(duì)于冰區(qū)航行船舶冰阻力及冰阻力影響下船舶運(yùn)動(dòng)的研究從未間斷。Kim等[1]模擬了碎冰條件下貨船的阻力性能，并將數(shù)值模擬結(jié)果與韓國(guó)水池的非冷凍模型冰試驗(yàn)結(jié)果以及加拿大冰池中的預(yù)鋸冰試驗(yàn)結(jié)果進(jìn)行了比較，取得了良好的一致性。Kajaste-Rudnitski和Kujala[2]采用商業(yè)軟件ABAQUS對(duì)層冰條件下的船舶三自由度運(yùn)動(dòng)進(jìn)行了模擬，結(jié)果顯示破冰是一個(gè)隨機(jī)的過(guò)程，冰阻力隨時(shí)間的變化呈現(xiàn)出一系列持續(xù)時(shí)間非常短的尖峰波動(dòng)。Valanto[3]通過(guò)理論計(jì)算和三維有限元模型預(yù)測(cè)了層冰條件下船舶受到的冰阻力，并進(jìn)一步考慮了船體型線及水的影響。Su[4]開(kāi)發(fā)了一種“點(diǎn)—面”接觸算法用來(lái)計(jì)算船舶表面的冰阻力，并在此基礎(chǔ)上通過(guò)自編程分析了船舶的阻力和機(jī)動(dòng)性。Lau和Sim?es Ré[5]利用離散元法軟件DECICE模擬了層冰和碎冰中破冰船的操縱性，并將其與模型試驗(yàn)結(jié)果進(jìn)行比較驗(yàn)證了離散元模型的準(zhǔn)確性。狄少丞等[6]基于離散元法模擬了破冰船在層冰和浮冰中的回轉(zhuǎn)運(yùn)動(dòng)，冰與船之間的接觸過(guò)程用球—三角接觸模型表示，分析了冰厚度和冰密集度對(duì)冰阻力和回轉(zhuǎn)直徑的影響。

上述的數(shù)值模擬研究主要集中于船—冰強(qiáng)碰撞過(guò)程中的船舶冰阻力，沒(méi)有考慮流體（海水）的作用，從而缺少流體對(duì)船舶冰阻力的影響研究。針對(duì)這一點(diǎn)，本文基于有限元分析中的流固耦合模型，模擬冰—水共同作用下船舶的回轉(zhuǎn)運(yùn)動(dòng)。

1 數(shù)值模擬參數(shù)

本文采用ANSYS/LS-DYNA顯式動(dòng)力學(xué)分析軟件，使用有限元法（FEM）模擬層冰中破冰船的回轉(zhuǎn)運(yùn)動(dòng)。建立了船、層冰、空氣域和海水域的三維模型，圖1所示為隱藏了空氣域后的數(shù)值模型示意圖。

圖 1 冰—水—船數(shù)值模型Fig. 1 Numerical model of ice, sea water and ship

模擬過(guò)程中，將船體視為剛體，船型的主要參數(shù)如表1所示。層冰的主要尺寸為600 m ×350 m × 1 m，冰網(wǎng)格采用1 m ×1 m ×1 m的六面體網(wǎng)格。已知冰層寬度（350 m）遠(yuǎn)大于船寬（22.6 m），邊界對(duì)碰撞結(jié)果的影響已經(jīng)很小，因此將冰的邊界條件設(shè)定為非反射邊界，從而模擬無(wú)限邊界的情況。由于冰材料的復(fù)雜性和隨機(jī)性，對(duì)于冰材料參數(shù)和模型的研究各有不同，目前較常用的一種做法是假設(shè)冰為各向同性彈塑性材料，如楊亮和馬駿[7]以及Liu[8]的研究。為此，本文也選擇此模型作為冰模型，相關(guān)參數(shù)的選取如表2所示，參數(shù)的選取及其合理性參見(jiàn)文獻(xiàn)[7]。當(dāng)單元的應(yīng)力或應(yīng)變超過(guò)有限元模型中的設(shè)定值時(shí)，該單元無(wú)效，將從模型中刪除，不影響后續(xù)計(jì)算。

表 1 船舶參數(shù)Table 1 Ship parameters

表 2 冰模型參數(shù)Table 2 Ice model parameters

空氣域尺寸為800 m × 500 m ×10 m，海水域尺寸為800 m ×500 m ×15 m。流體域的網(wǎng)格在接觸面連接，為保證流固耦合的精確度，采用等同于冰材料的1 m ×1 m ×1 m的六面體網(wǎng)格。對(duì)于流體的材料特性，采用軟件中的“空材料”描述，并通過(guò)定義流體狀態(tài)方程來(lái)描述壓力和密度之間的關(guān)系。以海水為例，狀態(tài)方程為Gruneisen狀態(tài)方程[9]，是一種可壓縮材料的狀態(tài)方程，其公式為

需指出的是，由于有限元軟件中沒(méi)有采用Naiver-Stokes方程求解流體運(yùn)動(dòng)，所以不能直接獲得流體動(dòng)力，不過(guò)通過(guò)該模型可以考慮冰在海水中的浮力。通過(guò)設(shè)定流體和固體的粘性與速度的耦合關(guān)系，將固體與流體耦合在一起。在考慮重力和動(dòng)態(tài)阻尼的同時(shí)，設(shè)定海水域和空氣域邊界為無(wú)限邊界，通過(guò)耦合固體和流域接觸點(diǎn)的位移和形變，來(lái)防止流固分離。模型中的船體和層冰模型均與流體進(jìn)行了相關(guān)的耦合。流固耦合方面的設(shè)置采用的是LS-DYNA軟件中的Arbitrary Lagrange-Euler（ALE）算法，通過(guò)定義關(guān)鍵字的方法將流體和固體耦合在一起，流固耦合關(guān)鍵字選擇“CONSTRAINED_LAGRANGE_IN_SOLID”。此后，設(shè)置流體與固體的耦合方式，LS-DYNA軟件提供了多種耦合方式[10]，本文選用“罰函數(shù)耦合法”（penalty coupling），耦合的方向設(shè)定為接觸節(jié)點(diǎn)的法向。此外，在流固耦合分析中還有很多必要的處理，如多物質(zhì)單元的耦合方式等，鑒于篇幅限制，這里不一一介紹，具體可參見(jiàn)文獻(xiàn)[10-11]。

在模型中，設(shè)定船舶的三自由度分別為縱蕩、橫搖和艏搖。初始時(shí)刻，隨船坐標(biāo)與全局坐標(biāo)一致，其中隨船坐標(biāo)系的坐標(biāo)原點(diǎn)為船舶重心，如圖2(a)所示，隨船坐標(biāo)系的x軸指向船艏，y軸指向左舷，z軸指向上方。圖2(b)所示為船舶艏部有限元模型的近視圖。數(shù)值方法和建模的更多細(xì)節(jié)可以參考文獻(xiàn)[12]。

圖 2 船舶坐標(biāo)和有限元模型示意圖Fig. 2 Sketch of coordinates and finite element model of the ship

2 網(wǎng)格收斂性分析及驗(yàn)證

2.1 網(wǎng)格收斂性分析

網(wǎng)格收斂性分析對(duì)于數(shù)值模擬具有重要意義。然而，由于船—冰—海水模型中存在大量網(wǎng)格的原因，很難減小網(wǎng)格尺寸。為了簡(jiǎn)化計(jì)算，選擇了一個(gè)類似的冰—結(jié)構(gòu)碰撞模型來(lái)檢查網(wǎng)格的收斂性。

用這個(gè)簡(jiǎn)單的模型模擬了鋼板擠壓帶有半球頭冰柱的過(guò)程。在模型中，鋼板設(shè)置為剛體，鋼板和冰的材料參數(shù)及破壞模式與第1節(jié)相同。鋼板的速度為1 m/s，冰柱底部固定。具體模型如圖3所示。

圖 3 平板—冰柱模型Fig. 3 Model of plate-ice cylinder

為了驗(yàn)證網(wǎng)格的收斂性，在不同冰模型網(wǎng)格尺寸下進(jìn)行了數(shù)次模擬。設(shè)置了3種網(wǎng)格：100 mm×100 mm×100 mm，90 mm×90 mm×90 mm和80 mm×80 mm×80 mm。3種網(wǎng)格尺寸下的鋼板冰阻力曲線如圖4所示。由圖可以看出，不同網(wǎng)格尺寸的冰阻力曲線隨時(shí)間變化逐漸重合。此外，無(wú)論網(wǎng)格尺寸如何變化，冰阻力的劇烈振蕩總是存在，這表明振蕩來(lái)自平板和冰柱碰撞，而不是網(wǎng)格尺寸。

圖 4 不同網(wǎng)格密度下的冰阻力Fig. 4 Ice resistance simulated by different mesh sizes

2.2 驗(yàn)證船與冰碰撞

在考慮流域影響之前，為了驗(yàn)證船和冰碰撞的數(shù)值模型，將數(shù)值模擬得到的船舶冰阻力與Lindqvist[13]經(jīng)驗(yàn)公式計(jì)算出的冰阻力進(jìn)行比較。保留與Lindqvist相同的部分輸入，例如船速、冰厚、船型主要參數(shù)等。除此外，一些冰參數(shù)，如塑性硬化模量、體積模量和塑性破壞應(yīng)變等在Lindqvist的方程中未使用，這些數(shù)據(jù)的選取需要根據(jù)先前的研究和經(jīng)驗(yàn)值進(jìn)行判斷。在船和冰碰撞驗(yàn)證中，取船速 V=5 m/s，其他參數(shù)同上一節(jié)。Lindqvist的公式為

此時(shí)，冰的邊界條件仍為無(wú)反射邊界條件。數(shù)值模擬結(jié)果與經(jīng)驗(yàn)公式的比較如圖5所示。

由圖5可以看出，數(shù)值模擬的平均值非常接近由經(jīng)驗(yàn)公式計(jì)算的值。經(jīng)驗(yàn)公式的結(jié)果為1.056 MN，數(shù)值模擬的平均值為1.055 MN，相對(duì)誤差為0.1%。因此，可以認(rèn)為數(shù)值模型的準(zhǔn)確度在一定程度上得到了驗(yàn)證。

為進(jìn)一步驗(yàn)證數(shù)值模擬的合理性，分別對(duì)1 m冰厚、不同航速，以及5 m/s航速、不同冰厚的情況進(jìn)行了模擬，并將模擬結(jié)果與簡(jiǎn)化后的Lindqvist經(jīng)驗(yàn)公式結(jié)果進(jìn)行了對(duì)比，如圖6所示，結(jié)果顯示模擬結(jié)果基本合理。

2.3 驗(yàn)證流體與結(jié)構(gòu)的相互作用

圖 6 不同速度和冰厚條件下數(shù)值模擬結(jié)果和經(jīng)驗(yàn)公式值對(duì)比 Fig. 6 Comparison of simulation results and empirical formulae under different speed and ice thickness conditions

在2.2節(jié)的基礎(chǔ)上研究流體與結(jié)構(gòu)的相互作用，在考慮流體存在的情況下，模擬層冰上的船舶冰阻力，并與經(jīng)驗(yàn)公式計(jì)算值進(jìn)行比較。所有參數(shù)與第2節(jié)相同。如式（2）所示，浸沒(méi)阻力包括潛在能量損失和摩擦阻力損失。在層冰破碎后，碎冰沿著船舶滑動(dòng)會(huì)引起摩擦阻力。在本文的模擬中，當(dāng)冰模型單元的應(yīng)力或應(yīng)變超過(guò)有限元模型中的設(shè)定值時(shí)，冰模型單元將被判定為無(wú)效并從模型中刪除，因而無(wú)法模擬碎冰沿船體滑動(dòng)，也無(wú)法計(jì)算破冰后的摩擦阻力。因此，在本文的流體—結(jié)構(gòu)相互作用模型中，只包含潛在能量損失。由此，方程（2）簡(jiǎn)化為

數(shù)值模擬結(jié)果與經(jīng)驗(yàn)公式的比較如圖7所示。經(jīng)驗(yàn)公式的結(jié)果為1.481 MN，數(shù)值模擬的平均值為1.466 MN，相對(duì)誤差為1%。因此，可以認(rèn)為流體—結(jié)構(gòu)相互作用方法在一定程度上得到了驗(yàn)證。比較圖5和圖7可以清楚地看到，冰阻增加了近28%，表明考慮海水對(duì)冰的支持作用后，增加了船舶的冰阻力。這說(shuō)明海水的影響不可輕易忽視，應(yīng)該在數(shù)值模擬中予以考慮。此外，對(duì)比圖5和圖7可知，當(dāng)采用流固耦合方法進(jìn)行分析時(shí)，冰阻力曲線的振蕩程度降低了，表明海水緩沖了船與冰之間的碰撞，降低了冰阻力的隨機(jī)性。

3 船與冰碰撞的數(shù)值模擬

3.1 不考慮流體域時(shí)的數(shù)值模擬

在數(shù)值模型驗(yàn)證的基礎(chǔ)上，模擬船體在層冰條件下的水平回轉(zhuǎn)運(yùn)動(dòng)。對(duì)于回轉(zhuǎn)運(yùn)動(dòng)的模擬，假設(shè)船在以恒定角速度轉(zhuǎn)彎之前以恒定速度直線前進(jìn)。這里需要指出的是，由于操縱性涉及螺旋槳和舵力等復(fù)雜的推進(jìn)力，本文并沒(méi)有直接解決操縱性問(wèn)題。但強(qiáng)制轉(zhuǎn)彎運(yùn)動(dòng)期間的模擬結(jié)果對(duì)于船舶操縱性仍然具有參考意義，因?yàn)樗鼈兙哂蓄愃频钠票Ｊ健Ｈ肿鴺?biāo)系中的速度可表示為：

圖9所示為全局坐標(biāo)系中，流固耦合作用下冰阻力隨時(shí)間的變化。從圖中可以看出，冰阻力呈現(xiàn)出一系列持續(xù)時(shí)間非常短的高峰，Kajaste-Rudnitski等[2]在冰船碰撞模擬中也觀察到了這一點(diǎn)，且在回轉(zhuǎn)運(yùn)動(dòng)中，冰阻力振蕩的振幅比在直線運(yùn)動(dòng)中更大。

所示。右舷的冰阻力是由船左轉(zhuǎn)時(shí)船體和冰航道內(nèi)部之間擠壓所產(chǎn)生。類似的現(xiàn)象也曾在其他學(xué)者基于離散元法[14]的研究中出現(xiàn)過(guò)。

3.2 考慮流體域時(shí)的數(shù)值模擬

為了研究海水域?qū)Ρ枇Φ挠绊?，比較無(wú)流固耦合船—冰碰撞模擬和流固耦合船—冰碰撞模擬的冰阻力情況。

圖 12 設(shè)計(jì)水線平面上冰阻力示意圖Fig. 12 Sketch of ice resistance on the designed waterline plane

4 結(jié) 論

本文基于非線性有限元法和流固耦合模型，模擬層冰中船舶的回轉(zhuǎn)運(yùn)動(dòng)，以研究水對(duì)冰阻力的影響，主要得出以下結(jié)論：

1） 海水的存在增加了船體在直線運(yùn)動(dòng)和回轉(zhuǎn)運(yùn)動(dòng)中各個(gè)方向上的冰阻力，但每個(gè)方向的增量都不同。

2） 回轉(zhuǎn)運(yùn)動(dòng)階段，在海水的作用下，受海水影響最大的是縱向冰阻力，其次是橫向和垂向冰阻力。

3） 當(dāng)采用流固耦合方法進(jìn)行分析時(shí)，冰阻力曲線的振蕩程度降低，表明海水緩沖了船與冰之間的碰撞，降低了冰阻力的隨機(jī)性。