基于近場動力學(xué)的破冰阻力預(yù)報方法研究

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

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

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

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

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

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

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

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

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

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

(1)

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

物質(zhì)點x的位置與時間t的物質(zhì)點運動方程為：

(2)

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

離散后每個物質(zhì)點x的運動方程可以表示為[17]：

(3)

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

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

(4)

式(4)中可以看出伸長率s的值僅取決于位移的大小，與方向無關(guān)，材料滿足各項同性的要求。破壞的判定條件是當(dāng)伸長率達(dá)到一定的值，認(rèn)為物質(zhì)點間發(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為材料破壞時的能量釋放率。可以通過斷裂力學(xué)中張開型裂紋的公式進(jìn)行推導(dǎo)，計算公式為[18]：

(9)

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

(10)

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

(11)

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

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

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

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

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

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

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

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

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

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

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

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

(12)

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

(13)

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

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

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

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

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

(14)

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

(15)

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

(16)

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

(17)

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

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

3 算例分析

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

圖8給出該極地船在層冰工況下試驗場景。

圖8 層冰阻力試驗場景

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

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

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

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

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

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

為了驗證數(shù)值模型的可行性，與模型試驗現(xiàn)象進(jìn)行比較，如圖9。通過數(shù)值模擬現(xiàn)象觀測，在破冰過程中，不斷有因彎曲破壞和擠壓破壞形成層冰的斷裂，在船肩處形成大量的彎曲破壞，因彎曲破壞形成的碎冰塊較大。而在船艏縱舯剖面附近與冰直接接觸的區(qū)域發(fā)生擠壓破壞，因擠壓破壞形成的碎冰塊小，破壞程度更高，為了現(xiàn)象的直觀清晰，本文剔除了因擠壓破壞導(dǎo)致破壞程度高于0.975的海冰物質(zhì)點。模型試驗和數(shù)值模擬彎曲破壞和擠壓破壞現(xiàn)象，如圖9(d)。通過比較圖9(e)，可以發(fā)現(xiàn)數(shù)值模擬的裂紋擴(kuò)展與試驗現(xiàn)象有一定的相似性，即在大面積彎曲破壞后，船體逐漸駛?cè)氡鶎?，船肩冰排開始出現(xiàn)局部的小型環(huán)形破壞，同時冰排在船體肩部稍后位置兩側(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ì)點具有一定的速度，并缺乏水的阻力作用，故不斷地沿著原有速度方向運動，逐漸遠(yuǎn)離船底。

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

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

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

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

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

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

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

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

圖11 破冰阻力實時曲線

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

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

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

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

4 結(jié)論

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

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