基于FEM/SPH耦合方法的極地運輸船舶冰阻力預報研究

王志鵬，郝寨柳，田于逵，

（中國船舶科學研究中心，江蘇無錫214082）

0 引 言

北極航道是聯(lián)系亞歐美三大洲的最短航線，與通常的蘇伊士運河和巴拿馬運河航線相比，能夠縮短數(shù)千公里航程，由此將大大降低商業(yè)運輸成本，影響和意義深遠。隨著氣候變暖，海冰融化加速，北極航道通航條件不斷改善，未來對極地運輸船舶的需求會明顯增大。極地運輸船舶不僅要滿足開敞水域性能和經(jīng)濟指標，而且需具備一定的破冰能力以滿足冰區(qū)航行需求。在極地運輸船舶開發(fā)過程中，冰阻力預報是船舶破冰能力和動力需求分析的重要輸入，也是總體設計的關鍵環(huán)節(jié)之一。在船舶冰阻力預報方面，主要有模型試驗、經(jīng)驗公式、數(shù)值模擬等方法，其中，模型試驗是最可靠的方式，基于實船測量和模型試驗數(shù)據(jù)總結(jié)的經(jīng)驗公式也被廣泛應用，目前常用的經(jīng)驗方法有Lindqvist 方法、Kei?nonen 方法、Riska 方法和Spencer 方法等[1-4]。相比于經(jīng)驗公式中模型過分簡化，數(shù)值計算能夠考慮海冰破壞模式和更多船型信息對冰阻力的影響。大量學者對海冰數(shù)值模型和破冰航行過程開展了研究，主要包括有限元、離散元、半經(jīng)驗數(shù)值模型等方法，使數(shù)值計算的精度和適用范圍得到了持續(xù)改進和提高。Arne Gürtner 等[5]基于有限元方法，利用LS-DYNA 軟件基于黏結(jié)單元建立海冰數(shù)值模型，對冰層與燈塔作用過程進行了數(shù)值模擬，獲得了作用區(qū)域載荷分布和變化過程，并與實測結(jié)果進行了對比，驗證了黏結(jié)單元法模擬海冰破壞的可行性；Wang 等[6]利用LS-DYNA 軟件，通過用戶自定義冰材料，將海冰破壞前表述為線彈性，采用多表面失效準則進行破壞判斷，并將數(shù)值計算結(jié)果和破冰船實測結(jié)果進行了對比，驗證了數(shù)值方法的可行性；任奕舟等[7]基于有限元理論，采用可破碎泡沫模型模擬冰材料，并通過冰錐試驗結(jié)果驗證其可行性，同時采用經(jīng)驗公式計算附連水質(zhì)量并疊加在船體質(zhì)量上，通過對破冰船不同航速下破冰過程進行數(shù)值模擬進而預報了船舶冰阻力；Su 等[8]采用離散元數(shù)值積分方法，耦合求解冰力和船體三自由度運動方程，建立船-冰作用時船體載荷數(shù)值模型，分析了船舶在破冰航行過程中冰阻力和船體局部載荷情況；季順迎等[9]研究了不同形態(tài)海冰離散單元模型特點，對極區(qū)海冰動力特性、海冰重疊堆積、海冰-海洋結(jié)構(gòu)物相互作用開展了數(shù)值模擬及分析；Tan 等[10]提出了基于半經(jīng)驗方法的數(shù)值模型，考慮壓力-面積關系計算船舶與冰層間接觸力，以此研究了船舶六自由度運動對破冰模式的影響，并通過實尺度破冰船數(shù)據(jù)進行了對比驗證；王鈺涵等[11]采用船冰接觸過程理想化、擠壓-彎曲過程理想化和海冰破損幾何形狀理想化等理想化破冰假設，對直航情況下連續(xù)破冰過程進行模擬，獲得了冰力時歷曲線和船舶運動響應，分析了相關破冰參數(shù)對連續(xù)破冰作用下破冰形狀和平均冰力的影響程度。

在上述冰阻力數(shù)值預報方法中，計算模型通常考慮船冰間相互作用，采用船冰水耦合求解的研究相對較少。開展船冰水耦合作用數(shù)值計算，能夠獲得破冰過程中冰層破壞模式、冰塊運動行為、船舶破冰排冰現(xiàn)象等更多信息，可為船舶性能預報和船型設計優(yōu)化提供更為直接有效的支撐，具有重要的理論意義和實際工程價值。本文基于有限元/光滑粒子流體動力學耦合算法，初步實現(xiàn)了破冰過程中的船/冰/水耦合過程模擬，并對船舶冰阻力進行了預報。

1 數(shù)值模型

海冰采用彈塑性模型，塑性屈服前海冰的力學行為描述為

式中：σ 為海冰應力；K = E/3( 1- 2υ )為海冰體積彈性模量；G =E/3( 1+ 2υ )為海冰剪切彈性模量；E 為海冰彈性模量；υ 為海冰泊松比；Ddel為膨脹應變；D′e= De- DeI/2 為彈性應變偏量；De為海冰彈性應變。塑性模型屈服準則采用修正的Drucker-Prager 帽蓋模型。船舶破冰航行過程中，冰層內(nèi)發(fā)生應力集中產(chǎn)生裂紋，隨著裂紋的擴張，冰層出現(xiàn)斷裂破壞現(xiàn)象。為模擬海冰裂紋的產(chǎn)生，引入粘聚單元損傷模型，即將冰層離散成孤立實體單元，在相鄰單元間插入粘聚單元作為潛在的斷裂面，如圖1所示，達到臨界強度時，粘聚單元失效，冰層內(nèi)出現(xiàn)破壞，其中，實體單元采用三棱柱單元，粘聚單元采用零厚度六面體單元。粘聚單元應力與相對應變在達到黏結(jié)強度之前表現(xiàn)為線彈性：

圖1 單元間粘聚方式示意圖Fig.1 Illustration of the ice mesh topology

式中，t表示應力，δ表示位移，下標n表示法向，下標s和t表示兩個切線方向。通過單元剛度矩陣K定義應力和位移間耦合關系，如果僅對Knn、Kss和Ktt參數(shù)定義，設定其他剛度系數(shù)值為0，即表示法向和切向為非耦合本構(gòu)關系。

船舶破冰過程中，海水對海冰力學行為的影響不可忽視。本文采用光滑粒子流體動力學方法模擬海水環(huán)境，實現(xiàn)破冰航行過程的船冰水動態(tài)耦合求解。光滑粒子法是一種無網(wǎng)格的純拉格朗日方法，采用一系列任意分布的粒子質(zhì)點來代表連續(xù)介質(zhì)流體，通過核函數(shù)W 來定義一定光滑長度h范圍內(nèi)其他鄰近粒子對目標粒子的影響程度，則流體域內(nèi)的速度、壓力及其梯度等變量均可通過一組無序質(zhì)點的核函數(shù)插值集合表示。

流體力學基本方程組離散為如下近似形式：

式中，dρi/dt、dUi/dt分別為粒子i相應的隨體導數(shù)。

有限元網(wǎng)格與水粒子之間通過罰函數(shù)接觸算法傳遞接觸狀態(tài)和作用力[12-13]，接觸面的法向接觸力表達式為

式中，ki是接觸剛度，fs是ki的比例因子，li是粒子i 相對網(wǎng)格單元的穿透深度，ni是網(wǎng)格單元的法向單位矢量。接觸面的切向接觸力為fti= μ ||fs，其中μ是粒子與網(wǎng)格單元的摩擦系數(shù)。

計算得到耦合作用力后，將作用力以外力的形式加入到流體動力方程和有限元動力方程中：

圖2 船/冰/水耦合求解過程Fig.2 The solution procedure of a ship-ice and water coupled problem

船冰水耦合求解方法如圖2 所示，船體為剛體模型，船冰間進行接觸判斷并計算相互作用力，同時根據(jù)水粒子的分布和位置計算水載荷，得到冰單元應力，當應力達到損傷起始值時采用損傷模型進行單元破壞判定，當達到破壞極限時，粘聚單元失效，更新碎冰單元接觸面使其能夠參與新的接觸判斷進而計算碎冰單元運動，同時更新冰單元邊界，如此迭代實現(xiàn)船冰水耦合求解。

2 海冰強度數(shù)值模擬

海冰壓縮強度和彎曲強度是影響船舶冰阻力預報的主要參數(shù)，采用本文建立的海冰數(shù)值模擬方法，參考ITTC推薦規(guī)程[14]中海冰物理強度測試方法，建立相應尺寸海冰試驗數(shù)值模型。采用相同加載方式，選取斯瓦爾巴群島海域當年冰為目標海冰，根據(jù)實測統(tǒng)計結(jié)果[15]，選取冰樣厚度為0.6 m，壓縮強度為620 kPa，彎曲強度為470 kPa。分別開展單軸壓縮強度和彎曲強度計算，驗證本文建立的海冰數(shù)值模擬方法，海冰數(shù)值模型主要參數(shù)見表1。

表1 海冰材料參數(shù)表Tab.1 Material properties of ice

單軸壓縮強度計算：參考物理模型試驗方式，冰樣尺寸為1h × 2h × 4h，h為冰樣厚度，采用上下壓盤對冰樣進行固定和載荷施加，模擬計算的海冰破壞對比及載荷時歷曲線如圖3 所示，可以得出：海冰破壞形式與試驗相符，采用σc= F/( Wh )計算得到海冰壓縮強度為637.5 kPa，與實測統(tǒng)計壓縮強度（620 kPa）相比，數(shù)值模擬結(jié)果偏大2.8%。

圖3冰樣單軸壓縮壓縮試驗數(shù)值模擬Fig 3 Numerical simulation of the uniaxial compression test of ice sample

彎曲強度計算：采用三點彎曲方式驗證彎曲強度，梁的尺寸為1h × 2h × 6h，采用三個圓柱體對冰樣進行固定和載荷施加，模擬計算的海冰破壞對比及載荷時歷曲線如圖4 所示，可以得出：海冰破壞形式與試驗相符，采用σf= M/W = 3Fl/( )2bh2計算得到海冰彎曲強度為491.3 kPa，與實測統(tǒng)計彎曲強度（470 kPa）相比，數(shù)值模擬結(jié)果偏大4.5%。

圖4 冰樣三點彎曲試驗數(shù)值模擬Fig.4 Numerical simulation of the three point bending test of ice sample

綜上，采用本文建立的海冰數(shù)值模擬方法，建立海冰數(shù)值模型，開展冰樣壓縮破壞和彎曲破壞數(shù)值模擬，海冰破壞形式與試驗相符，彎曲強度和壓縮強度與實測統(tǒng)計結(jié)果相差分別為2.8%和4.5%，海冰數(shù)值模型可作為開展船舶破冰航行過程模擬和冰阻力計算的輸入。

3 船舶冰阻力計算分析

本文以某極地船舶為研究對象，船型參數(shù)見表2，船體外形如圖5所示。

表2 船舶主尺度Tab.2 Principal dimensions of the carrier

圖5 船體外形圖Fig.5 Outside view of the ship

根據(jù)目標海域冰況特點，采用建立的海冰數(shù)值模型，選取1.3 m 和1.6 m 兩種厚度冰況，采用艏破冰方式，對船舶破冰過程進行計算，計算得到破冰過程和現(xiàn)象如圖6～7 所示。船舶艏部采用擠壓方式破冰前進，船舶前進方向，冰層出現(xiàn)徑向裂紋，隨著航行前進，冰層斷裂表現(xiàn)為環(huán)向裂紋，進而斷裂，破冰現(xiàn)象在下一次接觸破冰過程中循環(huán)出現(xiàn)；冰層破壞形成的冰塊被排向船體兩側(cè)，少量冰塊在浮力和摩擦力作用下沿船體滑移，部分冰塊在流場擾動作用下漂浮至航道中；破冰和排冰接觸作用構(gòu)成了冰阻力的不同組成部分。

圖7 冰層局部破壞現(xiàn)象Fig.7 Local damage of ice sheet

船舶破冰阻力時歷曲線如圖8 所示，初始破冰過程中，船舶靠近冰層利用艏部擠壓進行破冰，船冰接觸面積不斷增加，冰阻力逐漸增大；而后冰層受到持續(xù)破壞，粘聚單元持續(xù)失效造成接觸力卸載，冰阻力時歷曲線表現(xiàn)為高頻振蕩（曲線a），通過對阻力信號進行頻譜分析，發(fā)現(xiàn)其能量譜峰值主要集中在100 Hz 以下低頻范圍內(nèi)，高頻區(qū)域?qū)φw冰阻力的影響較為有限，采用低通濾波方式濾除高頻噪聲，獲得較為直觀的阻力（曲線b）。同時，采用Lindqvist公式[1,16]估算船舶冰阻力，假定70%的船體濕表面積被碎冰塊覆蓋，得到兩種冰厚下阻力曲線，如圖9所示。

通過分析船舶破冰現(xiàn)象和阻力變化曲線可知：（1）隨著航速和冰厚增加，船舶阻力明顯增大；數(shù)值計算結(jié)果中，阻力點與阻力-航速曲線偏差較大。數(shù)值計算考慮海冰的彎曲和剪切破壞，海冰破壞存在隨機性，造成了阻力點隨航速變化的離散性，使阻力點偏差增大。（2）對比經(jīng)驗公式與數(shù)值計算結(jié)果，1.3 m 厚度冰況阻力相差范圍為7.32%～13.91%，1.6 m 厚度冰況阻力相差范圍為5.13%～8.07%，數(shù)值計算得到的冰阻力均大于經(jīng)驗公式預報結(jié)果。Lindqvist 公式計算冰阻力時，航速與冰阻力是線性相關的，難以考慮冰水耦合作用及碎冰塊運動的影響，而數(shù)值計算對船冰水耦合作用進行模擬，冰阻力結(jié)果較大。（3）經(jīng)驗公式與數(shù)值計算的冰阻力差別隨著航速的增加而增加，隨著層冰厚度的增加而減小。船舶破冰過程中，經(jīng)驗公式與數(shù)值計算的冰阻力差別主要由冰水耦合及碎冰運動引起；隨著航速增加，冰層破壞后碎冰速度和碰撞程度均明顯增加，冰阻力差別隨航速增加明顯；隨著冰厚增加，冰塊破碎后尺寸增大，由于慣性原因碎冰塊運動受到海水影響減弱，引起冰阻力差別隨厚度的增加而降低。

圖8 船舶破冰阻力時歷曲線Fig.8 Time history of ice-breaking resistance

圖9 不同冰層厚度下破冰阻力和Lindqvist結(jié)果比較Fig.9 Comparison of the ice-breaking resistance of the polar vessel with Lindqvist’s results in ice with different thickness

4 結(jié) 語

本文基于有限元理論，采用粘聚單元法構(gòu)建海冰數(shù)值模型，開展了冰樣單軸壓縮強度和三點彎曲強度試驗過程數(shù)值模擬，驗證了該數(shù)值模型的合理性。在此基礎上，采用有限元/光滑粒子流體動力學耦合算法，建立船舶破冰過程中船/冰/水耦合作用數(shù)值計算方法，模擬破冰過程裂紋產(chǎn)生和擴展、冰塊翻轉(zhuǎn)與滑移等現(xiàn)象，針對某極地運輸船舶，開展了冰阻力數(shù)值計算和分析，得到不同航速、冰厚狀態(tài)下冰阻力變化情況。該方法可為船舶連續(xù)破冰模式下阻力預報提供一種技術(shù)手段，可為浮/碎冰工況數(shù)值計算方法研究提供參考。