金 強，張佳寧，張 雷，田于逵

（1.大連海事大學(xué)船舶與海洋工程學(xué)院，遼寧大連 116026；2.中國船舶科學(xué)研究中心，江蘇無錫 214082）

0 引 言

近年來，隨著全球氣候變暖，北極航線受到更廣泛的關(guān)注，北極資源的開采也成為眾多環(huán)北極國家爭相開展的目標(biāo)，從而推動了極地船舶的建造與研究。對于極地船舶的研究開展最多的是對航行于層冰冰況下的船舶冰阻力性能預(yù)報，而對航行于極地碎冰區(qū)船舶冰阻力性能的探究卻并不多見。在現(xiàn)已開通的北極航線中，近70%航段為碎冰所覆蓋，船體所受冰阻力是冰區(qū)船設(shè)計要考慮的主要因素。

對于極地船舶冰阻力的研究，Jones 等[1]采用“Healy”號破冰船模擬層冰冰況下船舶破冰阻力性能試驗，并將試驗結(jié)果與實船航行數(shù)據(jù)進行對比。DuBrovin 等[2]在1950—1955 年間四次船舶冰池試驗數(shù)據(jù)基礎(chǔ)上總結(jié)推導(dǎo)出了針對預(yù)報船模冰池試驗中所受冰阻力的計算公式。Kim 等[3]使用LS-DYNA軟件模擬了具有破冰能力貨船在碎冰航道下航行的冰阻力性能，并將結(jié)果與冰池試驗結(jié)果進行了對比。Wang 等[4]研究了層冰冰況下“Terry Fox”號破冰船的破冰阻力以及自航試驗結(jié)果，并將試驗結(jié)果與實船航行資料進行了對比。Selvadurai等[5]運用二維離散元模型模擬移動層冰與冰區(qū)錨泊結(jié)構(gòu)物間的相互影響。國內(nèi)對于冰區(qū)船舶研究起步較晚，郭春雨等[6]采用非凍結(jié)模型冰試驗方法研究了碎冰條件下船舶冰阻力隨速度、密集度變化的規(guī)律，并采用Colbourne 實船轉(zhuǎn)化方法預(yù)報了不同冰厚下的冰阻力。李紫麟[7]采用離散元方法，用三維圓盤模型模擬碎冰，將數(shù)值模擬結(jié)果與非凍結(jié)模型冰試驗結(jié)果進行對比，分析影響冰阻力大小因素。本文根據(jù)冰池試驗規(guī)程及相應(yīng)的冰級規(guī)范要求，構(gòu)建浮/碎冰的離散元數(shù)值模型及用于模擬計算的數(shù)值冰水池，采用離散元方法（discrete element method，DEM）數(shù)值計算航行于碎冰航道中船舶所受的冰阻力，并將數(shù)值模擬結(jié)果與HSVA冰池試驗數(shù)據(jù)及Dubrovin經(jīng)驗公式方法的結(jié)果進行對比分析，進而驗證數(shù)值模擬方法的精度。

1 碎冰航道船舶冰阻力數(shù)值計算

1.1 數(shù)值計算理論基礎(chǔ)

離散元法方法由Cundall 提出，其核心為牛頓第二運動定律基本思想，將離散分布的介質(zhì)看作是一系列離散單元的集合[8-9]。極地碎冰區(qū)的海冰具有較強的離散分布特性，可以將碎冰區(qū)的海冰視為一系列離散分布的單元集合。對于單個碎冰的離散元模型采用若干個離散粒子填充的方式創(chuàng)建出具有一定物理形狀的碎冰模型，填充的基本粒子為球形顆粒。構(gòu)成碎冰的粒子單元i和單元j間簡化受力如圖1 所示。碎冰-船體間接觸時，假定與組成碎冰的球形單元間發(fā)生接觸的船體表面視為半徑R趨向于無窮大的球形單元。每個離散單元均滿足相應(yīng)的運動方程，一個單元求解的輸出作為下一單元求解的輸入，通過時間步迭代求解各個單元的運動方程進而求得碎冰的整體運動狀態(tài)。

圖1 碎冰單元i，j間接觸示意Fig.1 Schematic diagram of contact between brash ice units of i and j

離散元算法的關(guān)鍵在于對單元間接觸力的計算，構(gòu)成碎冰的球形單元i和單元j間接觸的等效物理模型如圖2所示。

圖2 碎冰單元i，j間接觸的等效物理模型Fig.2 Equivalent physical model of contact between brash ice units of i and j

式中，Kt和Nt分別表示切向彈簧剛度和阻尼，Cfs為摩擦系數(shù)，dn和dt分別表示在單元接觸點處法向和切向上的重疊量，vn和vt分別表示在接觸點處的速度沿法向和切向上的速度分量。

式中，Nn,damp和Nt,damp分別表示法向和切向阻尼系數(shù)，定義如下：

式中，Cn,rest和Ct,rest分別為法向和切向恢復(fù)系數(shù)[9]，Req、Meq、Eeq及Geq分別表示相對于兩接觸單元i和j的半徑、質(zhì)量、楊氏模量以及剪切模量。

式中，vi和vj分別表示接觸的海冰粒子單元i和j的泊松比；Ei和Ej分別表示相接觸的海冰粒子單元i和j的楊氏模量。

1.2 船-冰數(shù)值模型及數(shù)字碎冰場構(gòu)建

用于船模試驗的某冰區(qū)油船數(shù)值模型如圖3所示。

圖3 某冰區(qū)油船數(shù)值模型Fig.3 Numerical model of a tanker

根據(jù)芬蘭-瑞典冰級規(guī)范（FSICR）[14]中Class IA 冰級對碎冰冰厚的規(guī)定及試驗中生成的海冰厚度確定數(shù)值計算中創(chuàng)建的碎冰厚度。創(chuàng)建碎冰離散元數(shù)值模型，首先繪制碎冰模型草圖如圖4所示，初步設(shè)定碎冰的尺寸形狀為圖4（a）中0.08 m×0.08 m的方形，圖4（b）為基于圖4（a）中草圖繪制的不規(guī)則多邊形。根據(jù)所要創(chuàng)建的模型冰冰厚尺寸，以冰厚為37.1 mm 為例對草圖4（a）和（b）拉伸相應(yīng)的厚度尺寸得到碎冰模型的物理實體，如圖5所示。為在碎冰域中生成尺寸大小不同的碎冰模型，依據(jù)冰池試驗中碎冰尺寸大小，對于生成的離散元碎冰模型采用羅辛-拉姆勒隨機分布方法將模型冰尺寸在0.08 m×0.08 m 的基礎(chǔ)尺寸上在[0.04，0.1]m 內(nèi)隨機分布，采用不同大小的球形基礎(chǔ)顆粒對圖5 中生成的物理模型進行粒子填充，填充粒子數(shù)為50。生成的碎冰離散元模型如圖6所示，冰池試驗與數(shù)值計算中碎冰的分布如圖7和圖8所示。

圖4 創(chuàng)建的碎冰物理模型草圖Fig.4 Sketch of greated physical model of brash ice

圖5 生成的碎冰物理模型實體Fig.5 Generated brash ice physical model entity

圖6 碎冰的離散元模型Fig.6 Discrete element model of the brash ice

圖7 冰池試驗中碎冰形狀及分布Fig.7 Shape and distribution of brash ice in the ice tank

圖8 數(shù)字碎冰場中碎冰分布Fig.8 Distribution of brash ice in the numerical brash ice field

網(wǎng)格的劃分基礎(chǔ)尺寸選為0.215 m。為提高計算精度，在數(shù)值碎冰場自由液面處疊加三層網(wǎng)格加密體，并且只在垂直方向（Z方向）上對網(wǎng)格尺寸進行細化。三層網(wǎng)格控制體網(wǎng)格尺寸依次取基礎(chǔ)尺寸的10%、20%和40%，即三層網(wǎng)格控制體區(qū)域在垂直方向尺寸分別為0.021 5 m、0.043 m 和0.086 m。在船首尾部曲率變化較大部位以及船體兩側(cè)與碎冰相互作用的部位同樣添加控制體對網(wǎng)格進行加密處理。船首尾部及船體表面水線面處的控制體設(shè)置如圖9 所示，控制體區(qū)域內(nèi)網(wǎng)格加密方式設(shè)置為各向同性，對X、Y、Z三個方向進行相同比例的加密處理，網(wǎng)格尺寸取為基礎(chǔ)尺寸的10%，即0.021 5 m，同樣對于船體兩側(cè)水線面處網(wǎng)格的加密尺寸取基礎(chǔ)尺寸的10%。生成網(wǎng)格如圖10和圖11所示。

圖9 船首尾及水線面處網(wǎng)格加密控制體Fig.9 Encryption control body of the grid at the bow,stern and waterline surface

圖10 計算域中縱剖面處網(wǎng)格Fig.10 Grid at longitudinal section in the computational domain

圖11 船體局部網(wǎng)格加密處理Fig.11 Encrypted grid of the corresponding parts on the ship

根據(jù)試驗條件及ITTC冰池試驗相關(guān)規(guī)定[15]，對于數(shù)字碎冰場的構(gòu)造：船首前方為1.5Lpp，船尾后方為2.5Lpp，Lpp為船舶的垂線間長，碎冰航道寬度為兩倍船寬。航道兩側(cè)創(chuàng)建壁面區(qū)域模擬試驗中的層冰區(qū)域，整個冰場長、寬、高分別為46.0 m、16.0 m和10.0 m，如圖12所示，初始化的數(shù)字碎冰航道見圖13。

圖12 數(shù)字碎冰場的三維尺寸Fig.12 Three-dimensional size of numerical brash ice field

圖13 計算域初始化后的碎冰航道Fig.13 Brash ice channel after initialization of calculation domain

2 試驗?zāi)Ｐ团c方案設(shè)計

2.1 相似定律

冰池試驗中為確保船模與實船的相似性，船模試驗相似數(shù)（傅汝德數(shù)、雷諾數(shù)及柯西數(shù)）均需考慮，

式中，V表示速度，g表示重力加速度，L表示長度，ν表示運動粘度，ρ表示水密度，E表示楊氏模量。

由于在模型試驗中傅汝德數(shù)與雷諾數(shù)不能夠同時滿足，且對于冰池試驗而言流體粘性影響相對較小，因此冰模試驗中通常忽略雷諾數(shù)相似。根據(jù)傅汝德數(shù)與柯西數(shù)相似規(guī)律，船體、碎冰及流體主要的幾何及物理量相似關(guān)系滿足表1所示，其中λ為縮尺比。

表1 各物理量相似規(guī)律Tab.1 The similarity laws of physical quantities

2.2 試驗?zāi)Ｐ?/h3>

試驗采用某冰區(qū)油船為研究對象，測試冰況依據(jù)芬蘭-瑞典冰級規(guī)范（FSICR）[14]選取Ice Class1A冰級條件，碎冰航道寬度為兩倍船寬尺寸。船模首尾部如圖14 所示（本章節(jié)所用船模試驗說明圖示均引用自文獻[16]），船體主要參數(shù)見表2。

表2 某冰區(qū)油船主參數(shù)表Tab.2 Main parameters of a certain oil tanker

圖14 試驗?zāi)Ｐ痛孜膊渴疽鈭DFig.14 Schematic diagram of the bow and stern of the model ship

2.3 冰池試驗方案設(shè)計

2.3.1 凍結(jié)模型冰碎冰航道構(gòu)造

冰池試驗相對于一般敞水條件下的船模試驗的技術(shù)更加復(fù)雜，碎冰航道的生成是試驗前期準備最為耗時且技術(shù)要求復(fù)雜的工作。根據(jù)HSVA 標(biāo)準模型冰制造程序生成厚度滿足要求的層冰，在層冰冰厚達到設(shè)定要求后冰池室內(nèi)溫度開始上升，當(dāng)溫度接近-3℃時，切割出寬度等于兩倍船寬大小的航道并采用冰鑿等設(shè)備將航道內(nèi)的層冰破碎構(gòu)造初始碎冰航道。碎冰制備及生成的碎冰航道如圖15和圖16所示。

圖15 碎冰制備Fig.15 Generation of the brash ice

圖16 生成的碎冰航道Fig.16 Generated brash ice channel

2.3.2 試驗工況

試驗過程中流體設(shè)為靜止，即不考慮實際冰況下碎冰漂移速度及水流速對船舶所受阻力的影響。將模型船航速設(shè)為定值0.464 m/s，將碎冰冰厚設(shè)置為37.1 mm、39.8 mm、46.3 mm 和49.4 mm，探究船舶所受冰阻力隨冰厚的變化規(guī)律。試驗中采用拖曳小車的方式使模型船勻速航行于碎冰航道中，如圖17所示。

圖17 碎冰航道中的模型船F(xiàn)ig.17 Model ship in the brash ice channel

3 數(shù)據(jù)分析與處理

3.1 數(shù)值模擬過程與冰池試驗對比分析

數(shù)值模擬中采用相對坐標(biāo)系的方式將船體視為靜止，設(shè)置碎冰及水流以恒定速度V=0.464 m/s（實船航速5 kn）向船舶運動。以模型冰厚hice=39.8 mm 計算為例，碎冰航道中船舶航行如圖18 和圖19 所示。

圖18 試驗中碎冰航道中船舶航行狀態(tài)Fig.18 Navigation status of the ship in brash ice channel in test

圖19 數(shù)字碎冰場中船模航行狀態(tài)Fig.19 Navigation status of the ship in numerical brash ice field

船首及船體兩側(cè)水線部位為船冰間的主要接觸區(qū)域，航行于碎冰航道時船冰間主要發(fā)生三種相互作用形式：

（1）船首部位與碎冰間的碰撞、堆積與擠壓，這也是船體所受冰阻力產(chǎn)生的主要因素。

（2）船首前方碎冰被不斷推離至船體兩側(cè)，與船體兩側(cè)表面發(fā)生摩擦與擠壓作用并在船尾后方留下等同于船寬大小的無碎冰分布航道，試驗與數(shù)值計算中相應(yīng)的現(xiàn)象分別如圖20和圖21所示。

圖20 試驗中船體后方碎冰航道Fig.20 Brash ice channel behind the ship in the test

圖21 數(shù)字碎冰場中船尾后方碎冰航道Fig.21 Brash ice channel behind the ship in the numerical brash ice field

（3）碎冰自船首部位被壓入船體以下，由于浮力作用碎冰緊貼船體底部表面運動直至船尾處浮出水面，試驗與數(shù)值計算中相應(yīng)的現(xiàn)象分別如圖22和圖23所示。

圖22 冰池試驗中壓入船底的碎冰Fig.22 Brash ice pressed into the bottom of the ship in the test

圖23 數(shù)字碎冰場中壓入船底碎冰Fig.23 Brash ice pressed under the bottom of the ship in the numerical brash ice field

在數(shù)值模擬過程中，通過將船冰間三種典型的相互作用形式與冰池試驗過程中相應(yīng)的現(xiàn)象進行對比可知，采用粒子噴射的方式創(chuàng)建的數(shù)字碎冰場能夠有效地模擬出船冰間的各種相互作用，并且與試驗現(xiàn)象相吻合。

3.2 數(shù)值計算及經(jīng)驗公式結(jié)果與試驗數(shù)據(jù)對比分析

數(shù)值模擬計算中為得到穩(wěn)定結(jié)果，船模應(yīng)航行至少兩倍船長距離。選取t=200 s時刻不同冰厚條件下船舶所受阻力時歷曲線如圖24所示。

圖24 不同冰厚條件下船體所受冰阻力時歷曲線Fig.24 Time history curve of ice resistance on hull under different ice thicknesses

從圖24中可以明顯得出：

（1）不論在何種冰厚條件下，船舶所受冰阻力特點均是幅值隨時間波動較大。這是因為作用于船體表面的冰阻力并不是一個恒定的外力作用，隨著船舶在碎冰航道中的行進，船首部位與碎冰間的作用形式不斷變化，發(fā)生碰撞的碎冰數(shù)量也隨時間而改變。

（2）四種冰厚條件下，船舶所受冰阻力變化的總體趨勢均隨時間增長而增大，隨后又趨于平穩(wěn)。在船首初次接觸碎冰到整個船體進入碎冰航道的過程中，船冰間發(fā)生相互作用的船體表面面積與碎冰數(shù)量不斷增加，因而在這一過程中船舶所受冰阻力呈現(xiàn)不斷增長趨勢；而隨著船體整個進入碎冰分布均勻的航道后，與船體表面發(fā)生碰撞的碎冰數(shù)及船體表面積均為定值，從而船舶所受冰阻力數(shù)值也保持穩(wěn)定。

1970年DuBrovin提出針對碎冰試驗中模型尺度下船舶所受阻力的經(jīng)驗公式[2]：

式中，R表示模型船所受冰阻力；p1與p2分別為系數(shù)，取決于碎冰分布密度及航道寬度，F(xiàn)r為傅汝德數(shù)，n為功率系數(shù)，取決于船型，參數(shù)A和Φ定義[17]如下：

式中，r為碎冰分布密度，ρice為碎冰密度，f0為船冰間摩擦系數(shù)，α0和αH分別表示模型船水線面入射角的一半與船舶前體菱形系數(shù)。具體參數(shù)選取見表3。

表3 DuBrovin經(jīng)驗公式各參數(shù)選取Tab.3 The parameters for DuBrovin empirical formula

根據(jù)式（14）～（16）及表3 確定的參數(shù)可計算得出，在冰厚分別為37.1 mm、39.8 mm、46.3 mm 和49.4 mm條件下時，模型船相應(yīng)地所受阻力大小為49.3 N、51.51 N、56.57 N和58.91 N。

對數(shù)值計算所得冰阻力數(shù)值進行處理，得出四種不同冰厚條件下船舶所受冰阻力均值，并將試驗數(shù)據(jù)及經(jīng)驗公式所得結(jié)果整理，如表4所示。

表4 三種方法所得冰阻力數(shù)值Tab.4 Ice resistance values obtained by three methods

根據(jù)表4數(shù)據(jù)可知三種方法所得船舶冰阻力變化趨勢如圖25所示。

根據(jù)圖25及表4數(shù)據(jù)分析，數(shù)值計算結(jié)果與DuBrovin經(jīng)驗公式所得的結(jié)果更加接近，但后者所得冰阻力數(shù)值普遍高于前者，且在冰厚為37.1 mm、39.8 mm、46.3 mm 和49.4 mm 四種條件下與冰池試驗數(shù)據(jù)誤差分別為9.24%、8.85%、24.33% 和25.38%，而數(shù)值計算所得結(jié)果在相應(yīng)冰厚條件下與試驗結(jié)果誤差分別為12.60%、11.33%、28.64%和29.45%。從數(shù)據(jù)中可以看出，在冰厚尺寸較小條件下，經(jīng)驗公式方法與離散元算法計算所得阻力數(shù)值與試驗數(shù)據(jù)誤差較小，有較高的精度；隨著冰厚增大，這兩種計算方法的誤差開始增大。這是由于對于經(jīng)驗公式方法而言，式（14）～（16）為DuBrovin 總結(jié)1950—1955 年間的四次船模冰池試驗數(shù)據(jù)推導(dǎo)而得，原始數(shù)據(jù)相對較少并且從公式中可以看出冰厚因素hice對阻力計算并未起到主要影響作用，因此當(dāng)冰厚增大時阻力計算的誤差也相應(yīng)增大；就離散元算法而言，船冰間的接觸力主要取決于接觸力Fc的分量，即法向力Fn與切向力Ft的計算，而從理論部分的式（1）可知，對于碎冰單元間接觸力計算影響最大的為單元間法向與切向的剛度與阻尼系數(shù)，與碎冰厚度有關(guān)的質(zhì)量M對接觸力計算影響不大，因此隨著冰厚增加所得的阻力值變化并不明顯。

圖25 三種計算方法所得阻力變化趨勢Fig.25 Variation trend of resistance obtained by three methods

4 結(jié) 論

本文目的在于依據(jù)試驗中冰池構(gòu)造，采用粒子噴射的方式創(chuàng)建數(shù)字碎冰場，并且使用基于離散元算法的數(shù)值計算方法模擬船舶在冰厚為37.1 mm、39.8 mm、46.3 mm 和49.4 mm 四種條件下所受冰阻力，將計算結(jié)果及經(jīng)驗公式數(shù)據(jù)分別與模型船冰池試驗數(shù)據(jù)進行對比分析，驗證數(shù)值計算方法精度及DuBrovin經(jīng)驗公式方法估算船舶在碎冰區(qū)域所受冰阻力的可靠性，得出如下結(jié)論：

（1）對比冰池試驗，在數(shù)值模擬中船冰間相互作用形式與試驗中現(xiàn)象吻合一致。表明基于粒子噴射方式創(chuàng)建的數(shù)字碎冰場是合理的，能夠準確地模擬出船舶航行于碎冰區(qū)船冰間的物理現(xiàn)象。

（2）使用DuBrovin 經(jīng)驗公式計算四種冰厚條件下的冰阻力，與試驗結(jié)果相比較的誤差分別為9.24%、8.85%、24.33%和25.38%，數(shù)據(jù)誤差在一定的合理范圍之內(nèi)，從而表明該經(jīng)驗公式在冰阻力初步估算與預(yù)報中具有一定的準確度與參考意義。

（3）基于離散元算法的數(shù)值模擬方法與經(jīng)驗公式計算結(jié)果吻合較好，最大誤差僅為5.70%；在冰厚尺寸較小條件下與冰池試驗數(shù)據(jù)誤差也在合理范圍之內(nèi)，即在模型冰厚為37.1 mm 和39.8 mm 時誤差僅為12.60%和11.33%；而當(dāng)模型冰厚增大為46.3 mm 和49.4 mm 時，由于離散元算法機理導(dǎo)致數(shù)值計算誤差分別增至28.64%和29.45%。從而，采用基于離散元算法的數(shù)值計算方法，在冰厚較薄情況下計算船體所受冰阻力能夠達到較高的精度，而在碎冰厚度較厚時該方法的誤差偏大。