王超， 曹成杰， 熊偉鵬， 葉禮裕， 汪春輝

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

冰阻力是影響極地航行船安全的主要因素。在航行過程中，打破有限厚度的海冰，保證安全的航線是其主要任務(wù)[1]。針對冰阻力問題，研究者進(jìn)行了大量的研究工作。實(shí)船試驗(yàn)在冰阻力值的測量和冰層破壞模式的觀測上有著不可取代的作用[2-3]，但代價(jià)太大，逐漸被日益成熟的、結(jié)果準(zhǔn)確、現(xiàn)象可靠的模型試驗(yàn)而取代[4-7]。但受限于試驗(yàn)條件，模型試驗(yàn)難以滿足全部研究者的需求。隨著計(jì)算機(jī)技術(shù)的逐漸成熟，數(shù)值模擬成本降低，參數(shù)設(shè)置靈活，結(jié)果較為準(zhǔn)確，使用范圍日益廣泛。

由于船-冰接觸過程會(huì)產(chǎn)生大量的冰塊斷裂、裂紋擴(kuò)展等非連續(xù)問題，采用傳統(tǒng)有限元方法難以準(zhǔn)確模擬冰的破壞形式及有效地處理冰的極端大變形問題，模擬出的現(xiàn)象往往不太理想[8-11]。而斷裂、破壞相關(guān)的數(shù)值模擬中采用較多的方法是基于離散元、光滑粒子流體動(dòng)力學(xué)及近場動(dòng)力學(xué)等無網(wǎng)格方法來實(shí)現(xiàn)的。其中，近場動(dòng)力學(xué)理論因其基于位移函數(shù)的積分形式來構(gòu)造基本的運(yùn)動(dòng)方程，構(gòu)建的運(yùn)動(dòng)方程不存在微分項(xiàng)，有效地避免類似連續(xù)介質(zhì)力學(xué)在處理不連續(xù)問題時(shí)存在發(fā)散的問題[12]。適用于計(jì)算均勻或非均勻材料結(jié)構(gòu)的斷裂、損傷、破碎和大尺度變形等問題。趙國良[13]應(yīng)用近場動(dòng)力學(xué)建立船艏-冰的相互作用模型，對連續(xù)破冰、碎冰航行等工況下的冰載荷影響因素進(jìn)行分析。陸錫奎[14]采用近場動(dòng)力學(xué)與有限元耦合方法對破冰船船艏連續(xù)破冰過程進(jìn)行了數(shù)值模擬。葉禮裕等[15-16]將近場動(dòng)力學(xué)方法應(yīng)用到冰-槳銑削等問題上。針對海冰破碎問題，近場動(dòng)力學(xué)方法有著獨(dú)特的適用性。

本文參考ITTC規(guī)范，將提出的區(qū)域接觸判別方法應(yīng)用到鍵型近場動(dòng)力學(xué)方法中，編寫極地航行船破冰航行中破冰阻力求解程序，建立了船-冰作用過程中的破冰阻力計(jì)算模型。基于模型試驗(yàn)結(jié)果驗(yàn)證數(shù)值方法的有效性。在此基礎(chǔ)上，開展不同航速下的破冰阻力預(yù)報(bào)工作。

1 近場動(dòng)力學(xué)理論基礎(chǔ)

近場動(dòng)力學(xué)將連續(xù)介質(zhì)離散為均勻或非均勻的物質(zhì)點(diǎn)。每一個(gè)物質(zhì)點(diǎn)都能承受一定的體載荷、速度、位移，并且會(huì)發(fā)生移動(dòng)和變形。理論上，任意兩物質(zhì)點(diǎn)間會(huì)存在相互作用，但當(dāng)超出一定的距離后，兩物質(zhì)點(diǎn)間的相互作用力較小。在近場動(dòng)力學(xué)里中，假定物質(zhì)點(diǎn)間的相互作用的近場域半徑為δ。每一物質(zhì)點(diǎn)的作用域受到近場域半徑大小δ影響。近場域半徑增大，作用范圍也增大。近場動(dòng)力學(xué)鄰域作用模型如圖1所示。

圖1 物質(zhì)點(diǎn)x的作用域模型

目前，鍵型近場動(dòng)力學(xué)材料中微觀彈脆性本構(gòu)模型發(fā)展較為成熟。在PMB鍵型本構(gòu)中，兩物質(zhì)點(diǎn)x和x′上的力是大小相等、方向相反作用在二者連線上的相互作用力。其鄰域內(nèi)作用模式如圖2所示。

力密度函數(shù)表示為：

(1)

式中:ξ表示兩物質(zhì)點(diǎn)的初始相對位置;η表示兩物質(zhì)點(diǎn)的相對位移量，η+ξ為任意時(shí)刻兩物質(zhì)點(diǎn)的相對位移。

物質(zhì)點(diǎn)x的位置與時(shí)間t的物質(zhì)點(diǎn)運(yùn)動(dòng)方程為：

(2)

式中：u為物質(zhì)點(diǎn)x的位移；u′為x′的位移；ρ為材料密度；Hx為物質(zhì)點(diǎn)x的近場域；V為近場域內(nèi)x′的體積；b(x,t)為物質(zhì)點(diǎn)受到的外力；式(1)、(2)共同構(gòu)成了材料的本構(gòu)關(guān)系。

離散后每個(gè)物質(zhì)點(diǎn)x的運(yùn)動(dòng)方程可以表示為[17]：

(3)

式中：n為時(shí)間步；下標(biāo)i和j代表物質(zhì)點(diǎn)的編號；i為要計(jì)算的物質(zhì)點(diǎn)；j為臨近xi的物質(zhì)點(diǎn)；Vj為物質(zhì)點(diǎn)j的體積。通過式(3)可以求得每個(gè)物質(zhì)點(diǎn)的加速度，從而得到每個(gè)物質(zhì)點(diǎn)的位移。

此外，為了描述材料的破壞情況，近場動(dòng)力學(xué)中提出極限伸長率的概念。即：

(4)

式(4)中可以看出伸長率s的值僅取決于位移的大小，與方向無關(guān)，材料滿足各項(xiàng)同性的要求。破壞的判定條件是當(dāng)伸長率達(dá)到一定的值，認(rèn)為物質(zhì)點(diǎn)間發(fā)生永久斷裂，該值稱為極限伸長率，記為s0。對于微觀彈脆性材料，力密度函數(shù)可以簡化為：

f(y(t),ξ)=g(s(t,ξ))μ(t,ξ)

(5)

式中g(shù)(s(t,ξ))為線性標(biāo)量函數(shù)，定義為：

(6)

式中κ為體積模量。

μ(t,ξ)為鍵斷裂判斷函數(shù)，定義為：

(7)

式中s0可表示為：

(8)

式中G0為材料破壞時(shí)的能量釋放率?？梢酝ㄟ^斷裂力學(xué)中張開型裂紋的公式進(jìn)行推導(dǎo)，計(jì)算公式為[18]：

(9)

式中KIC為斷裂強(qiáng)度因子，可以通過實(shí)驗(yàn)的方法來確定。季順迎等[19]基于海冰斷裂韌度實(shí)驗(yàn)，擬合了渤海域溫度變化對斷裂韌度的變化關(guān)系式。因此，G0可以表示為：

(10)

彈脆性材料在初始條件下是各向同性的，鍵斷裂之后會(huì)存在的各向異性情況。為反映變形后域內(nèi)鍵斷裂的程度，引入鍵的破壞程度參數(shù)，即：

(11)

2 破冰阻力數(shù)值計(jì)算模型

ITTC規(guī)范中極地航行船實(shí)際航行過程中阻力分為：冰阻力和水阻力。冰阻力又可詳細(xì)分為：冰層破壞引發(fā)的破冰阻力、碎冰隨海水沿船體滑動(dòng)引發(fā)的清冰阻力，以及被船艏擠壓在船體下方由于浮力引發(fā)的浸沒阻力[20]。本文針對極地航行船在層冰中航行的破冰阻力進(jìn)行研究，建立的破冰阻力數(shù)值計(jì)算模型基于以下幾個(gè)關(guān)鍵點(diǎn)進(jìn)行闡述。

2.1 數(shù)值模型基本假設(shè)

極地航行船舶實(shí)際航行時(shí)，冰面十分廣闊，但對航行性能影響較大的是位于航道及船側(cè)附近的冰面。數(shù)值模擬中理應(yīng)盡可能大地模擬層冰域，但限于方法和計(jì)算耗時(shí)的限制，本文參考ITTC冰阻力模型試驗(yàn)規(guī)范，建立數(shù)值模型層冰的計(jì)算域。

規(guī)范中，針對層冰阻力的模型試驗(yàn)給出的推薦長度是除了船長以外有1.5倍船長作為測量段[20]。而在數(shù)值模擬中，將層冰簡化為長方體，由于無需考慮加速、減速的長度，在考慮首尾邊界干擾，同時(shí)兼顧計(jì)算效率，計(jì)算域面積設(shè)定應(yīng)不小于船舶長寬的2倍。為了避免因?qū)挾仍O(shè)置的局限性，導(dǎo)致破壞后的冰層受船舶擠壓作用向外偏移，造成與真實(shí)情況不符情況，在冰層兩側(cè)添加了固定邊界以約束其運(yùn)動(dòng)。由于缺乏水環(huán)境的模擬，冰層破碎會(huì)后隨著水流運(yùn)動(dòng)，碎冰位置發(fā)生改變，通過添加重-浮力模型也不能準(zhǔn)確地實(shí)現(xiàn)浸沒阻力的模擬。故未添加重-浮力模型，以避免計(jì)算中產(chǎn)生的破冰阻力中含有浸沒阻力，從而影響破冰阻力計(jì)算的準(zhǔn)確性。此外，極地航行船對于船艏通常采用局部加強(qiáng)，在連續(xù)破冰航行中，破冰船艏部變形十分微小，本文忽略其變形帶來的影響，將船體當(dāng)作剛體處理。

2.2 數(shù)值模型離散化方法

以某極地航行船為研究對象，將船體表面三維模型離散為一系列四邊形網(wǎng)格單元。船艏是與海冰發(fā)生碰撞的主要區(qū)域，且船艏區(qū)域復(fù)雜，曲率變化大，為了提高計(jì)算精度，同時(shí)優(yōu)化計(jì)算的效率，針對船艏區(qū)域進(jìn)行網(wǎng)格加密處理，其他區(qū)域采用適當(dāng)?shù)木W(wǎng)格大小來劃分。劃分結(jié)果如圖3所示。

圖3 船體表面網(wǎng)格劃分結(jié)果

海冰模塊是由近場動(dòng)力學(xué)方法進(jìn)行建立。將層冰的計(jì)算域離散為物質(zhì)點(diǎn)形式。其中，鄰域半徑采用建議給出的δ=3.015 Δx為鄰域半徑[21]。船-冰相對位置依據(jù)船舶水線位置和不同冰厚情況下露出水面的厚度建立位置關(guān)系。建立的破冰阻力數(shù)值模型，如圖4所示。

圖4 破冰阻力數(shù)值計(jì)算模型

2.3 數(shù)值模型動(dòng)態(tài)接觸判別方法

在建立船-冰離散化數(shù)值模型后，為了保證接觸檢測的精度，通過考慮海冰物質(zhì)點(diǎn)和四邊形船體面元的相對位置關(guān)系進(jìn)行判斷。當(dāng)船體最近的面元中心駛向海冰物質(zhì)點(diǎn)達(dá)到一定距離后，開始進(jìn)入接觸識別區(qū)。過小的接觸識別區(qū)半徑會(huì)導(dǎo)致部分粒子穿透船體表面；過大的接觸識別區(qū)半徑會(huì)導(dǎo)致計(jì)算量的顯著增加，本文推薦的船舶面元的識別區(qū)半徑在1.5～2.5倍的單一步長船舶行進(jìn)距離為宜。對進(jìn)入接觸識別區(qū)的海冰物質(zhì)點(diǎn)，再通過當(dāng)前時(shí)刻和下一時(shí)刻與最近船體面元中心法向量的關(guān)系進(jìn)行接觸判別。船-冰數(shù)值模型接觸判別方法，如圖5所示。

(12)

即未發(fā)生接觸條件。而當(dāng)t+Δt時(shí)刻海冰物質(zhì)點(diǎn)穿透到船舶內(nèi)，如圖5(c)所示，則其滿足的關(guān)系為:

(13)

即發(fā)生接觸條件，需重新分配該穿透的海冰物質(zhì)點(diǎn)。

圖5 船-冰的接觸判別方法

2.4 數(shù)值模擬破冰阻力計(jì)算

圖6 發(fā)生穿透物質(zhì)點(diǎn)的重新定位

具體位置坐標(biāo)為：

(14)

重新分配后海冰物質(zhì)點(diǎn)i的速度為：

(15)

重新分配后海冰物質(zhì)點(diǎn)i受船舶的力為：

(16)

式中ρi和Vi分別代表海冰物質(zhì)點(diǎn)的密度和體積。將所有發(fā)生穿透的海冰物質(zhì)點(diǎn)i進(jìn)行累加，可以得到船舶所受冰阻力F為：

(17)

圖7給出采用接觸判別方法的數(shù)值模型模擬結(jié)果及局部細(xì)節(jié)。復(fù)雜的船艏區(qū)域未發(fā)生穿透現(xiàn)象，表明采用接觸判別算法的準(zhǔn)確性。

圖7 采用接觸判別的計(jì)算結(jié)果

3 算例分析

3.1 數(shù)值模型參數(shù)設(shè)定

圖8給出該極地船在層冰工況下試驗(yàn)場景。

圖8 層冰阻力試驗(yàn)場景

數(shù)值模擬采用與模型試驗(yàn)一致的縮尺比1∶40。滿足傅汝德數(shù)和柯西數(shù)相似準(zhǔn)則進(jìn)行縮尺。其中，冰強(qiáng)度、幾何長度、冰厚和冰彈性模量的縮尺比滿足1∶λ，時(shí)間和速度的縮尺比滿足1∶λ0.5，質(zhì)量和力的縮尺比滿足1∶λ3。參照模型試驗(yàn)比尺的船型參數(shù)，如表1所示。

表1 某極地航行船模型尺寸

參照模型比尺的海冰參數(shù)設(shè)置如表2所示。

表2 海冰參數(shù)原型及模型值

參考文獻(xiàn)[22]收斂性分析計(jì)算結(jié)果，本文設(shè)定時(shí)間步長為0.0 005 s，海冰物質(zhì)點(diǎn)直徑為L/900.0。

3.2 破冰阻力數(shù)值模擬現(xiàn)象

為了驗(yàn)證數(shù)值模型的可行性，與模型試驗(yàn)現(xiàn)象進(jìn)行比較，如圖9。通過數(shù)值模擬現(xiàn)象觀測，在破冰過程中，不斷有因彎曲破壞和擠壓破壞形成層冰的斷裂，在船肩處形成大量的彎曲破壞，因彎曲破壞形成的碎冰塊較大。而在船艏縱舯剖面附近與冰直接接觸的區(qū)域發(fā)生擠壓破壞，因擠壓破壞形成的碎冰塊小，破壞程度更高，為了現(xiàn)象的直觀清晰，本文剔除了因擠壓破壞導(dǎo)致破壞程度高于0.975的海冰物質(zhì)點(diǎn)。模型試驗(yàn)和數(shù)值模擬彎曲破壞和擠壓破壞現(xiàn)象，如圖9(d)。通過比較圖9(e)，可以發(fā)現(xiàn)數(shù)值模擬的裂紋擴(kuò)展與試驗(yàn)現(xiàn)象有一定的相似性，即在大面積彎曲破壞后，船體逐漸駛?cè)氡鶎?，船肩冰排開始出現(xiàn)局部的小型環(huán)形破壞，同時(shí)冰排在船體肩部稍后位置兩側(cè)出現(xiàn)一條向艏柱位置擴(kuò)展的環(huán)向裂紋；當(dāng)該環(huán)向裂紋接近船體舷側(cè)，沿船體擴(kuò)展的環(huán)向裂紋在船肩處發(fā)生局部破壞后，引發(fā)一條新的環(huán)向裂紋的形成與擴(kuò)展，顯現(xiàn)出由一系列環(huán)向裂紋所割裂出的大尺寸碎冰塊，如圖9(f)。

此外，通過觀察圖9數(shù)值模型的現(xiàn)象，可以發(fā)現(xiàn)，由于未添加重力和浮力模型，破碎后的海冰沒有貼近船體表面，而是由于接觸后重新分配的海冰物質(zhì)點(diǎn)具有一定的速度，并缺乏水的阻力作用，故不斷地沿著原有速度方向運(yùn)動(dòng)，逐漸遠(yuǎn)離船底。

圖9 數(shù)值模擬與模型試驗(yàn)現(xiàn)象(模擬航速3 kn，冰厚1.5 m)

3.3 數(shù)值模擬破冰阻力計(jì)算結(jié)果

在模型試驗(yàn)中，將航行阻力分解為冰破壞阻力和水下阻力2部分。通過觸覺式傳感器將整個(gè)船艏區(qū)域覆蓋測量冰破壞阻力和通過單項(xiàng)測力傳感器測量船舶總的航行阻力。并將總阻力與冰破壞阻力差值作為船舶水下阻力。本文與模型試驗(yàn)中冰破壞阻力進(jìn)行比較。

在數(shù)值模擬中，考慮邊界效應(yīng)的影響，在船舶進(jìn)入冰層和駛離冰層都留有0.25倍的船長作為邊界載荷干擾段，破冰阻力使用的數(shù)據(jù)段如圖10。

圖10 破冰阻力計(jì)算數(shù)據(jù)使用段

在連續(xù)破冰過程中，船-冰接觸表現(xiàn)為擠壓破壞和彎曲破壞反復(fù)循環(huán)過程，模型試驗(yàn)和數(shù)值模擬結(jié)果均顯示出船-冰接觸力變化劇烈，如圖11。在3 kn航速下，模型試驗(yàn)的冰破壞阻力均值為16.07 N，如圖11(a)；數(shù)值模擬的破冰阻力均值為14.45 N，如圖11(b)；破冰阻力均值誤差為10.08%。但數(shù)值模擬的阻力極大值較試驗(yàn)?zāi)M的更大，極小值較試驗(yàn)?zāi)M的更小，均值卻較模型試驗(yàn)小。造成該現(xiàn)象的原因可能在于：真實(shí)海冰受溫度、鹽度、孔隙度、內(nèi)部晶格排列方向等因素影響，模型試驗(yàn)制作的模型冰更貼近于真實(shí)海冰，而采用PMB本構(gòu)模型僅考慮了冰材料的脆性，忽視冰材料的塑性，導(dǎo)致冰材料斷裂時(shí)所受的阻力更大。此外，阻力成分劃分方式中，模型試驗(yàn)在有水環(huán)境下缺乏對水阻力成分進(jìn)行分離，現(xiàn)有的采集系統(tǒng)難以將水阻力和冰阻力進(jìn)行完全地分離采集，這可能是造成數(shù)值模擬破冰阻力更小的原因之一。

3.4 不同航速破冰阻力數(shù)值模擬

在此基礎(chǔ)下，進(jìn)行不同航速的數(shù)值模擬，不同航速下破冰過程的模擬現(xiàn)象如圖12所示。

圖11 破冰阻力實(shí)時(shí)曲線

圖12 不同航速下破冰阻力模擬現(xiàn)象

通過觀察不同航速現(xiàn)象，發(fā)現(xiàn)與船艏正面接觸的海冰斷裂形成的冰塊形狀、大小不一，有很大的隨機(jī)性，但航速越高，導(dǎo)致彎曲破壞形成的碎冰破壞程度更高，形成的碎冰塊尺寸更加細(xì)小，表明航速是影響連續(xù)破冰航行后形成碎冰尺寸的因素之一。此外，研究了不同航速下的破冰阻力，如圖13所示。

圖13 1.5 m冰厚條件下隨航速的變化的破冰阻力

隨著航速的增加，因彎曲破壞形成的碎冰尺寸更加細(xì)小，導(dǎo)致冰層斷裂破壞所需的能量增加，會(huì)引起船體在層冰中航行的破冰阻力增大。

4 結(jié)論

1)與試驗(yàn)航速進(jìn)行對比，發(fā)現(xiàn)海冰在受到航行船的作用后，會(huì)因彎曲、擠壓、裂紋擴(kuò)展等破壞形成的碎冰塊，與模型試驗(yàn)有一定的相似性；同時(shí)，破冰阻力計(jì)算結(jié)果與模型試驗(yàn)結(jié)果在量級趨于一致，誤差為10.08%，表明計(jì)算方法可行性；

2)在不同航速的數(shù)值模擬中，與船艏正面接觸的海冰斷裂形成的冰塊大小不一，有很大的隨機(jī)性，但航速越高，彎曲破壞形成的碎冰破壞程度越高，形成的碎冰尺寸越小，表明航速是影響海冰斷裂形成冰塊大小的因素之一；同時(shí)，船體在層冰中航行的破冰阻力隨航速的上升而增大。