計及流體影響的船舶回轉(zhuǎn)冰阻力數(shù)值模擬

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

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

0 引 言

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

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

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

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

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

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

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

表 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)格。對于流體的材料特性，采用軟件中的“空材料”描述，并通過定義流體狀態(tài)方程來描述壓力和密度之間的關(guān)系。以海水為例，狀態(tài)方程為Gruneisen狀態(tài)方程[9]，是一種可壓縮材料的狀態(tài)方程，其公式為

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

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

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

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

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

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

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

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

為了驗證網(wǎng)格的收斂性，在不同冰模型網(wǎng)格尺寸下進行了數(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)格尺寸的冰阻力曲線隨時間變化逐漸重合。此外，無論網(wǎng)格尺寸如何變化，冰阻力的劇烈振蕩總是存在，這表明振蕩來自平板和冰柱碰撞，而不是網(wǎng)格尺寸。

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

2.2 驗證船與冰碰撞

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

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

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

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

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

圖 6 不同速度和冰厚條件下數(shù)值模擬結(jié)果和經(jīng)驗公式值對比 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)驗公式計算值進行比較。所有參數(shù)與第2節(jié)相同。如式（2）所示，浸沒阻力包括潛在能量損失和摩擦阻力損失。在層冰破碎后，碎冰沿著船舶滑動會引起摩擦阻力。在本文的模擬中，當冰模型單元的應(yīng)力或應(yīng)變超過有限元模型中的設(shè)定值時，冰模型單元將被判定為無效并從模型中刪除，因而無法模擬碎冰沿船體滑動，也無法計算破冰后的摩擦阻力。因此，在本文的流體—結(jié)構(gòu)相互作用模型中，只包含潛在能量損失。由此，方程（2）簡化為

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

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

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

在數(shù)值模型驗證的基礎(chǔ)上，模擬船體在層冰條件下的水平回轉(zhuǎn)運動。對于回轉(zhuǎn)運動的模擬，假設(shè)船在以恒定角速度轉(zhuǎn)彎之前以恒定速度直線前進。這里需要指出的是，由于操縱性涉及螺旋槳和舵力等復(fù)雜的推進力，本文并沒有直接解決操縱性問題。但強制轉(zhuǎn)彎運動期間的模擬結(jié)果對于船舶操縱性仍然具有參考意義，因為它們具有類似的破冰模式。全局坐標系中的速度可表示為：

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

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

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

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

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

4 結(jié) 論

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

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

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

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