陳銳，黃武剛，陳曉璐，康美澤

中國(guó)艦船研究設(shè)計(jì)中心，湖北 武漢 430064

0 引 言

南、北極地理位置特殊，且擁有豐富的資源。開(kāi)發(fā)極地航道和極地資源，開(kāi)展極地科學(xué)考察具有重要的戰(zhàn)略意義。隨著全球變暖，北冰洋海冰融化，北極地區(qū)蘊(yùn)藏的豐富資源和極具價(jià)值的航道，也日益成為世界各國(guó)關(guān)注的焦點(diǎn)。極地科學(xué)考察破冰船是探索兩極的必要裝備，如何準(zhǔn)確評(píng)估極地科學(xué)考察破冰船連續(xù)破冰的冰阻力和破冰能力，是破冰船總體設(shè)計(jì)的關(guān)鍵技術(shù)之一，具有重要的理論和現(xiàn)實(shí)意義[1-2]。

在破冰船的冰阻力研究過(guò)程中，學(xué)者們采用了多種研究手段。國(guó)外早期基于模型試驗(yàn)、理論分析和實(shí)船觀測(cè)等方法，初步建立了一些破冰船冰阻力的經(jīng)驗(yàn)計(jì)算公式，如Lindqvist 公式[3]、Vance 公式[4]、Lewis 公式[5]、改進(jìn)的Edwards 公式[6]和Riska 公式[7]等，具有較高的實(shí)用意義，其中Lindqvist 公式使用最為廣泛。近年來(lái)，部分學(xué)者開(kāi)始使用數(shù)值模擬方法研究船與冰的相互作用，離散元方法（DEM）與有限元方法（FEM）為應(yīng)用最為廣泛的2 種數(shù)值方法。DEM 注重不同離散單元間的連接方式，Zhan 等[8]采用了離散元軟件DECICE 模擬了船舶在浮冰區(qū)的操縱運(yùn)動(dòng)，季順迎等[9]研究了海洋結(jié)構(gòu)物所受冰載荷的離散元算法并開(kāi)發(fā)了Ice DEM 軟件等。國(guó)內(nèi)在FEM 方面的研究也逐漸興起，王峰等[10-11]使用有限元軟件LS-DYNA 對(duì)連續(xù)破冰的過(guò)程進(jìn)行了數(shù)值模擬，何菲菲[12]采用商業(yè)軟件DYTRAN 對(duì)破冰船的破冰能力進(jìn)行了研究。冰池模型阻力試驗(yàn)是驗(yàn)證破冰阻力的重要方法，但投資成本大，試驗(yàn)周期長(zhǎng)。目前，國(guó)外具備冰池試驗(yàn)條件的有芬蘭阿克北極、德國(guó)漢堡船池，國(guó)內(nèi)具備條件的有天津大學(xué)。

在極地破冰考察船的設(shè)計(jì)過(guò)程中，需要對(duì)設(shè)計(jì)船型方案的破冰阻力和破冰能力進(jìn)行評(píng)價(jià)。本文將采用經(jīng)驗(yàn)公式法和數(shù)值模擬法進(jìn)行計(jì)算和結(jié)果比對(duì)，在天津大學(xué)冰水池開(kāi)展模型試驗(yàn)以驗(yàn)證以上結(jié)果，并對(duì)經(jīng)驗(yàn)公式和數(shù)值模擬計(jì)算方法的適用性和結(jié)果特點(diǎn)進(jìn)行研究，以便為快速評(píng)價(jià)破冰船的冰阻力提供可用的方法，為進(jìn)一步優(yōu)化改進(jìn)破冰船船型提供有效的設(shè)計(jì)手段。

1 破冰船破冰能力估算方法

1.1 常用破冰阻力估算方法

冰阻力是破冰船舶在冰區(qū)行進(jìn)過(guò)程中所受的主要載荷，對(duì)船舶的破冰阻力有著直接影響。研究人員基于實(shí)船測(cè)試以及理論分析，從不同角度推出了包括Vance 公式、Lindqvist 公式、Lewis 公式、Su 公式和Riska 公式在內(nèi)的一系列破冰阻力計(jì)算公式。

本文將采用Lindqvist，Vance 和Lewis 這3 種經(jīng)驗(yàn)公式進(jìn)行計(jì)算，探究不同經(jīng)驗(yàn)公式對(duì)冰阻力預(yù)報(bào)存在的差異，并對(duì)經(jīng)驗(yàn)公式進(jìn)行參數(shù)敏感性分析。經(jīng)驗(yàn)公式預(yù)報(bào)結(jié)果也用于后文與數(shù)值模擬結(jié)果以及船模試驗(yàn)結(jié)果進(jìn)行比對(duì)分析。

1.2 平整冰破冰阻力試算結(jié)果及分析

經(jīng)驗(yàn)公式所需參數(shù)如表1 所示。針對(duì)不同的航速與海冰厚度，采用經(jīng)驗(yàn)公式對(duì)冰阻力進(jìn)行了計(jì)算，具體工況與計(jì)算結(jié)果如表2 所示。對(duì)于3 種冰厚（1.5，2 及2.5 m），選取1～4 kn 不同航速，用3 種經(jīng)驗(yàn)公式進(jìn)行冰阻力計(jì)算，如圖1 所示，對(duì)于不同的工況，各經(jīng)驗(yàn)公式所得的計(jì)算結(jié)果都表明，冰阻力隨著船舶航速與冰厚的增加而增大，但具體數(shù)值存在差異。

圖1 不同工況下經(jīng)驗(yàn)公式預(yù)報(bào)冰阻力的比較Fig. 1 Comparison of the results obtained with empirical formula under different conditions

冰阻力的大小與破冰過(guò)程中船舶、冰層的參數(shù)大小密切相關(guān)。參數(shù)敏感性分析是將每個(gè)自變量分別在一定范圍內(nèi)變化，根據(jù)其相應(yīng)的變化趨勢(shì)以及變化程度，評(píng)估其變化量對(duì)結(jié)果的影響大小，參數(shù)敏感性分析即為對(duì)影響程度大小的衡量[11]。敏感性系數(shù)越大，則自變量對(duì)結(jié)果的影響也越大。根據(jù)表1 所示的經(jīng)驗(yàn)公式計(jì)算所需參數(shù)，選取了航速、冰厚以及彎曲強(qiáng)度3 個(gè)變量，針對(duì)航速為2 kn，冰厚1.5 m，船寬26 m，彎曲強(qiáng)度900 kPa 的初始條件，進(jìn)行參數(shù)敏感性分析，結(jié)果如圖2 所示。

圖2 參數(shù)敏感性分析Fig. 2 Parameter sensitivity analysis

表1 經(jīng)驗(yàn)公式計(jì)算所需參數(shù)Table 1 Main parameters required by calculation of ice resistance formulas

從圖中可以看出，若船舶航速增大25%，Lindqvist，Vance 和Lewis 公式對(duì)冰阻力的預(yù)報(bào)結(jié)果分別增大了10.15%，2.60%以及8.78%。航速的變化對(duì)Lindqvist 預(yù)報(bào)的冰阻力變化率影響最大，對(duì)Vance 公式的影響最小。冰厚的變化對(duì)3 種經(jīng)驗(yàn)公式的預(yù)報(bào)值均有較大影響，若冰厚增大25%，Lindqvist，Vance 和Lewis 公式對(duì)冰阻力的預(yù)報(bào)結(jié)果分別增大了39.83%，54.63%以及47.60%。冰的彎曲強(qiáng)度對(duì)Lindqvist 預(yù)報(bào)的冰阻力變化率影響最大，若彎曲強(qiáng)度增大20%，Lindqvist 預(yù)報(bào)結(jié)果增大13.96%，Vance 和Lewis 公式增長(zhǎng)不大，增長(zhǎng)率均在3%以內(nèi)。

表2 經(jīng)驗(yàn)公式計(jì)算結(jié)果Table 2 Results of ice resistance formulas calculation

2 破冰船冰阻力的數(shù)值模擬

本節(jié)在ANSYS 經(jīng)典模塊中進(jìn)行船模與冰層建模與網(wǎng)格劃分。采用非線性有限元軟件DYNA，基于傳統(tǒng)有限元方法與黏聚單元法（ICEM）構(gòu)造冰層數(shù)值模擬模型，對(duì)破冰船在平整冰區(qū)航行時(shí)連續(xù)破冰過(guò)程中的破冰阻力進(jìn)行了數(shù)值模擬。

數(shù)值模擬工況參考委托天津大學(xué)建筑工程學(xué)院對(duì)某船舶開(kāi)展的阻力模型試驗(yàn)，采用縮尺比λ為30 的船模，對(duì)不同冰厚與航速下船模的破冰阻力進(jìn)行數(shù)值模擬，以便與船模試驗(yàn)結(jié)果進(jìn)行比對(duì)分析。船模的主要尺度如表3 所示，其中La為總長(zhǎng)，Bw為水線寬， ?為排水量，模擬工況如表4 所示。

表3 船模的主要尺度Table 3 Main dimensions of the ship model

表4 數(shù)值模擬工況Table 4 Different simulation conditions

2.1 基 于 傳 統(tǒng) 有 限 元 法 的 破 冰 阻 力 數(shù) 值模擬

有限元數(shù)值模擬技術(shù)發(fā)展迅速，顯式動(dòng)力學(xué)有限元軟件DYNA 被廣泛應(yīng)用于解決碰撞接觸問(wèn)題，本文在應(yīng)用LS-DYNA 時(shí)，連續(xù)體的處理方法采用Lagrange 法[13]。

數(shù)值模型如圖3 所示，平整冰長(zhǎng)9 m，約為2 倍船長(zhǎng)，寬度為4.5 m，約為5 倍船寬，在平整冰豎直方向上僅設(shè)置一層網(wǎng)格來(lái)模擬平整冰的彎曲斷裂。在不與船殼接觸的三邊界上采用三邊固支對(duì)冰層進(jìn)行約束。在船冰碰撞的數(shù)值模型中，船體定義為殼單元，冰單元定義為實(shí)體單元。船網(wǎng)格采用四邊形面網(wǎng)格劃分，在碰撞過(guò)程中，不考慮船殼變形，故采用剛性材料MAT020-RIGID。由于冰材料的結(jié)構(gòu)特性非常復(fù)雜，被視為一種材料屬性和本構(gòu)關(guān)系都十分復(fù)雜的材料，國(guó)內(nèi)外學(xué)者在數(shù)值模擬研究中都沒(méi)有得出較理想的研究成果[14]。本文在參考國(guó)內(nèi)外多位學(xué)者的研究后，采用各向同性彈性斷裂失效本構(gòu)模型模擬冰材料[14-15]，模型參數(shù)的選取如表5 所示。其中， ρ為密度，SM為剪切模量，SIGY為屈服應(yīng)力，ETAN為塑性硬化模量，BULK為體積模量，EPF為塑性失效應(yīng)變，PRF為截?cái)鄩毫Α?/p>

圖3 船模和平整冰的網(wǎng)格劃分Fig. 3 Mesh patterns of the model and level ice

速度為0.188 m/s，冰厚為0.05 m 時(shí)，破冰船在連續(xù)破冰過(guò)程中所受冰阻力的時(shí)歷曲線如圖4 所示。從圖中可以看出，在0～3 s 階段，船體不斷擠壓冰單元，由于船冰接觸面積的不斷增大，冰阻力曲線隨時(shí)間的增長(zhǎng)振蕩上升。大約在3 s 之后，由于船艏完全進(jìn)入冰層，單位時(shí)間內(nèi)破冰船擠壓侵蝕的冰單元趨于穩(wěn)定，故冰阻力的均值亦趨于穩(wěn)定狀態(tài)。而在碰撞過(guò)程中，由于冰單元的失效，會(huì)出現(xiàn)接觸力卸載現(xiàn)象，故冰阻力會(huì)呈現(xiàn)振蕩狀態(tài)。

表5 冰體材料參數(shù)Table 5 Material parameters of the ice

圖4 冰阻力時(shí)歷曲線Fig. 4 Time history of ice resistance

2.2 基于黏聚單元法的破冰阻力數(shù)值模擬

黏聚單元模型以斷裂力學(xué)為理論基礎(chǔ)，廣泛應(yīng)用于固體材料裂紋的形成、增長(zhǎng)與穩(wěn)定分析。將黏聚單元應(yīng)用于冰層的數(shù)值模型的建立，可模擬平整冰與破冰船相互作用時(shí)發(fā)生的彎曲斷裂以及大裂紋的隨機(jī)擴(kuò)展過(guò)程[16-17]。

不同于傳統(tǒng)有限元模型中用標(biāo)準(zhǔn)六面體單元直接相連建立冰層，黏聚單元法將實(shí)體冰單元進(jìn)行離散，離散后的各個(gè)實(shí)體冰單元的內(nèi)邊界上插入共節(jié)點(diǎn)的零厚度黏聚單元，通過(guò)黏聚單元將相鄰的實(shí)體冰單元進(jìn)行連接。2 種單元相互獨(dú)立且具有不同的材料屬性，又相互連接傳遞受力。其中，實(shí)體冰單元根據(jù)冰材料的本構(gòu)關(guān)系發(fā)生彈性和塑性變形。黏聚單元是潛在的斷裂平面，當(dāng)黏聚單元承受的強(qiáng)度超出臨界斷裂強(qiáng)度時(shí)，黏聚單元失效，實(shí)體單元間出現(xiàn)縫隙而脫落破碎[11]。黏聚單元法所建立的冰層模型如圖5 所示。圖5（a）為平整冰模型，白色縱橫交替的為黏聚單元，藍(lán)色為實(shí)體冰單元；圖5（b）為平整冰模型的局部放大圖，藍(lán)色為實(shí)體冰單元。

圖5 黏聚單元法建立的冰層模型Fig. 5 Level ice model constructed by cohesive element method

船殼采用剛性材料MAT020-RIGID。對(duì)于實(shí)體冰單元可以采用許多簡(jiǎn)化的不同材料進(jìn)行模擬，例如：各向同性彈性模型、簡(jiǎn)化的破損模型、彈塑性模型等。本文主要考慮冰材料的彈性與塑性，實(shí)體冰單元采用MAT012-ISOTROPIC-ELASTICPLASTIC。黏聚單元采用MAT186-COHESIVEGENERAL 材料，具體各個(gè)部分材料參數(shù)如表6所示[10,16]，表中G為楊氏模量，PR為泊松比，S為屈服應(yīng)力，ETAN為塑性模量，TS為拉伸強(qiáng)度，SS為剪切強(qiáng)度，En為 法向斷裂能量釋放率，Es為剪切斷裂能量釋放率

表6 材料參數(shù)Table 6 Material parameters of the ship shell and ice

以冰厚為0.067 m，3 個(gè)不同航速下的工況的數(shù)值模擬結(jié)果為例。如圖7 所示，取冰阻力穩(wěn)定后的2 s 左右時(shí)間段內(nèi)的平均值作為數(shù)值模擬的船模所受冰阻力大小。對(duì)以上9 個(gè)工況的冰阻力進(jìn)行數(shù)值模擬計(jì)算。從圖中可以看到，在連續(xù)破冰的過(guò)程中，經(jīng)過(guò)一段時(shí)間后，模擬所得的冰阻力始終在一定范圍內(nèi)波動(dòng)，取這一穩(wěn)定區(qū)間的平均值為該工況下的冰阻力大小?？梢钥吹?，在同一冰厚下，數(shù)值模擬所得冰阻力的大小隨航速增大而增加。

圖7 冰厚0.067 m 下的數(shù)值模擬Fig. 7 Numerical simulations with 0.067 m ice thickness

3 破冰船冰阻力模型試驗(yàn)

3.1 試驗(yàn)?zāi)Ｐ团c環(huán)境條件模擬

不論是破冰船的初步設(shè)計(jì)還是對(duì)其水動(dòng)力性能的研究，冰水池的船模試驗(yàn)都至關(guān)重要，船模試驗(yàn)?zāi)軌蜃顬檎鎸?shí)和準(zhǔn)確地反映航行過(guò)程中的受力狀態(tài)[18]。本文中的破冰船的船模試驗(yàn)在天津大學(xué)低溫水池實(shí)驗(yàn)室開(kāi)展，試驗(yàn)?zāi)Ｐ团c數(shù)值模擬模型參數(shù)一致，具體尺度如表3 所示。

天津大學(xué)低溫水池實(shí)驗(yàn)室在模型試驗(yàn)中采用第2 代模型冰——尿素冰，這種冰模型采用“噴霧引晶”技術(shù)模擬天然海冰生長(zhǎng)過(guò)程，與天然海冰具有完全一致的分層紋理結(jié)構(gòu)，表層為細(xì)小、緊密的粒狀結(jié)晶層，下層為垂向分布的柱狀結(jié)晶層，粒狀層與柱狀層的耦合關(guān)系與天然海冰保持一致。模擬環(huán)境條件斷面與天然結(jié)構(gòu)對(duì)比如圖8 所示。

圖8 模擬環(huán)境條件斷面與天然結(jié)構(gòu)對(duì)比圖Fig. 8 Comparison of the simulation condition and realistic condition

試驗(yàn)中利用高精度的圖像采集設(shè)備對(duì)作用過(guò)程進(jìn)行水上與水下的同步實(shí)時(shí)記錄，以實(shí)現(xiàn)對(duì)破壞模式、斷裂長(zhǎng)度的變化及運(yùn)動(dòng)進(jìn)程進(jìn)行精確的觀測(cè)。冰阻力模型試驗(yàn)簡(jiǎn)圖如圖9 所示。

圖9 冰阻力模型試驗(yàn)Fig. 9 Model test in the towing tank

3.2 試驗(yàn)結(jié)果與分析

以船模試驗(yàn)結(jié)果為基準(zhǔn)，對(duì)經(jīng)驗(yàn)公式法以及數(shù)值模擬結(jié)果進(jìn)行了比對(duì)分析。由于文中采用的3 種經(jīng)驗(yàn)公式中的參數(shù)為實(shí)船試驗(yàn)所得結(jié)果[19]，在以往的文獻(xiàn)中均用于實(shí)船冰阻力預(yù)報(bào)，故本文對(duì)比分析的結(jié)果均為外推至實(shí)船尺度的冰阻力值大小。詳細(xì)數(shù)據(jù)如表7 所示。

以船模試驗(yàn)結(jié)果為基準(zhǔn)，圖10 為經(jīng)驗(yàn)公式法、FEM、ICEM 與船模試驗(yàn)所得冰阻力結(jié)果的對(duì)比圖。從圖中可以看出，Vance 與Lewis 這2 種經(jīng)驗(yàn)公式對(duì)冰阻力結(jié)果的預(yù)報(bào)非常保守，誤差在47%以上。Lindqvist 預(yù)報(bào)結(jié)果較好，與數(shù)值方法所得結(jié)果接近，且誤差均在30%以內(nèi)。Lindqvist公式、船模試驗(yàn)法以及數(shù)值模擬方法所得冰阻力大小均隨航速與與冰厚的增加呈上升趨勢(shì)。

圖10 不同方法所得結(jié)果的對(duì)比圖Fig. 10 Comparison of the error based on different methods

對(duì)于數(shù)值方法而言，F(xiàn)EM 與ICEM 在厚度較小的情況下，冰阻力預(yù)報(bào)結(jié)果更準(zhǔn)確，在冰厚較大（如2.5 m）時(shí)，誤差較大在25% 左右。在冰厚較小航速較高的情況下，如1.5 和2 m 時(shí)，ICEM方法預(yù)報(bào)的冰阻力值較傳統(tǒng)有限元法更為準(zhǔn)確，與船模試驗(yàn)結(jié)果相比精度誤差在10%以內(nèi)。

表7 各方法所得冰阻力匯總Table 7 Results of ice forces obtained by different methods

4 結(jié) 論

本文采用經(jīng)驗(yàn)公式法、數(shù)值模擬法以及船模試驗(yàn)法對(duì)破冰船在平整冰中連續(xù)航行時(shí)的冰阻力進(jìn)行了預(yù)報(bào)。數(shù)值模擬法中分別采用傳統(tǒng)有限元方法與黏聚單元法構(gòu)建冰層，以船模試驗(yàn)結(jié)果為基準(zhǔn)，對(duì)比分析了數(shù)值模擬結(jié)果以及經(jīng)驗(yàn)公式法所得冰阻力預(yù)報(bào)結(jié)果。為了能夠最為真實(shí)和準(zhǔn)確地反映航行過(guò)程中的受力狀態(tài)，本文中的破冰船的船模試驗(yàn)在天津大學(xué)低溫水池實(shí)驗(yàn)室開(kāi)展。得出結(jié)論如下：

1） 冰阻力隨航速、船寬、冰厚以及彎曲強(qiáng)度4 個(gè)變量的增大均呈上升趨勢(shì)，其中冰厚對(duì)冰阻力的影響最大。3 種經(jīng)驗(yàn)公式中，Lindqvist 公式的預(yù)報(bào)結(jié)果與船模試驗(yàn)結(jié)果更為接近，而Vance和Lewis 公式更為保守。

2） 數(shù)值模擬法可以較為準(zhǔn)確地進(jìn)行冰阻力預(yù)報(bào)，預(yù)報(bào)值與船模試驗(yàn)結(jié)果較為吻合。在冰厚為1.5 m 時(shí)可以采用傳統(tǒng)有限元法對(duì)冰阻力進(jìn)行預(yù)報(bào)，其預(yù)報(bào)結(jié)果精度在可接受范圍內(nèi)且結(jié)果偏保守，適于工程應(yīng)用。在冰厚為2 m 時(shí)，可以采用黏聚單元法，因?yàn)槠漕A(yù)報(bào)精度最高，在3 種航速下誤差均在10%以內(nèi)。

本文驗(yàn)證了數(shù)值模擬法和經(jīng)驗(yàn)與解析公式預(yù)報(bào)的適用性和準(zhǔn)確性。說(shuō)明采用傳統(tǒng)有限元法與黏聚單元法構(gòu)造冰層數(shù)值模擬模型對(duì)破冰過(guò)程中的冰阻力預(yù)報(bào)較為可靠，且文中選取的冰材料參數(shù)也與實(shí)際參數(shù)相對(duì)符合，可用于后續(xù)的工程應(yīng)用。在實(shí)際的冰阻力預(yù)報(bào)中，可以結(jié)合經(jīng)驗(yàn)公式法與數(shù)值模擬的方法，以兼顧預(yù)報(bào)結(jié)果的準(zhǔn)確性與高效性。