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

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

0 引 言

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

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

1 數(shù)值模型

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

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

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

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

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

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

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

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

式中，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ìn)行接觸判斷并計算相互作用力，同時根據(jù)水粒子的分布和位置計算水載荷，得到冰單元應(yīng)力，當(dāng)應(yīng)力達(dá)到損傷起始值時采用損傷模型進(jìn)行單元破壞判定，當(dāng)達(dá)到破壞極限時，粘聚單元失效，更新碎冰單元接觸面使其能夠參與新的接觸判斷進(jìn)而計算碎冰單元運動，同時更新冰單元邊界，如此迭代實現(xiàn)船冰水耦合求解。

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

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

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

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

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

彎曲強(qiáng)度計算：采用三點彎曲方式驗證彎曲強(qiáng)度，梁的尺寸為1h × 2h × 6h，采用三個圓柱體對冰樣進(jìn)行固定和載荷施加，模擬計算的海冰破壞對比及載荷時歷曲線如圖4 所示，可以得出：海冰破壞形式與試驗相符，采用σf= M/W = 3Fl/( )2bh2計算得到海冰彎曲強(qiáng)度為491.3 kPa，與實測統(tǒng)計彎曲強(qiá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ù)值模擬，海冰破壞形式與試驗相符，彎曲強(qiáng)度和壓縮強(qiáng)度與實測統(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ù)目標(biāo)海域冰況特點，采用建立的海冰數(shù)值模型，選取1.3 m 和1.6 m 兩種厚度冰況，采用艏破冰方式，對船舶破冰過程進(jìn)行計算，計算得到破冰過程和現(xiàn)象如圖6～7 所示。船舶艏部采用擠壓方式破冰前進(jìn)，船舶前進(jìn)方向，冰層出現(xiàn)徑向裂紋，隨著航行前進(jìn)，冰層斷裂表現(xiàn)為環(huán)向裂紋，進(jìn)而斷裂，破冰現(xiàn)象在下一次接觸破冰過程中循環(huán)出現(xiàn)；冰層破壞形成的冰塊被排向船體兩側(cè)，少量冰塊在浮力和摩擦力作用下沿船體滑移，部分冰塊在流場擾動作用下漂浮至航道中；破冰和排冰接觸作用構(gòu)成了冰阻力的不同組成部分。

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

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

通過分析船舶破冰現(xiàn)象和阻力變化曲線可知：（1）隨著航速和冰厚增加，船舶阻力明顯增大；數(shù)值計算結(jié)果中，阻力點與阻力-航速曲線偏差較大。數(shù)值計算考慮海冰的彎曲和剪切破壞，海冰破壞存在隨機(jī)性，造成了阻力點隨航速變化的離散性，使阻力點偏差增大。（2）對比經(jīng)驗公式與數(shù)值計算結(jié)果，1.3 m 厚度冰況阻力相差范圍為7.32%～13.91%，1.6 m 厚度冰況阻力相差范圍為5.13%～8.07%，數(shù)值計算得到的冰阻力均大于經(jīng)驗公式預(yù)報結(jié)果。Lindqvist 公式計算冰阻力時，航速與冰阻力是線性相關(guān)的，難以考慮冰水耦合作用及碎冰塊運動的影響，而數(shù)值計算對船冰水耦合作用進(jìn)行模擬，冰阻力結(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ù)值模型，開展了冰樣單軸壓縮強(qiáng)度和三點彎曲強(qiáng)度試驗過程數(shù)值模擬，驗證了該數(shù)值模型的合理性。在此基礎(chǔ)上，采用有限元/光滑粒子流體動力學(xué)耦合算法，建立船舶破冰過程中船/冰/水耦合作用數(shù)值計算方法，模擬破冰過程裂紋產(chǎn)生和擴(kuò)展、冰塊翻轉(zhuǎn)與滑移等現(xiàn)象，針對某極地運輸船舶，開展了冰阻力數(shù)值計算和分析，得到不同航速、冰厚狀態(tài)下冰阻力變化情況。該方法可為船舶連續(xù)破冰模式下阻力預(yù)報提供一種技術(shù)手段，可為浮/碎冰工況數(shù)值計算方法研究提供參考。