基于近場動力學的破冰阻力預報方法研究

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

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

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

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

本文參考ITTC規范，將提出的區域接觸判別方法應用到鍵型近場動力學方法中，編寫極地航行船破冰航行中破冰阻力求解程序，建立了船-冰作用過程中的破冰阻力計算模型。基于模型試驗結果驗證數值方法的有效性。在此基礎上，開展不同航速下的破冰阻力預報工作。

1 近場動力學理論基礎

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

圖1 物質點x的作用域模型

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

力密度函數表示為：

(1)

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

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

(2)

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

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

(3)

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

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

(4)

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

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

(5)

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

(6)

式中κ為體積模量。

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

(7)

式中s0可表示為：

(8)

式中G0為材料破壞時的能量釋放率。可以通過斷裂力學中張開型裂紋的公式進行推導，計算公式為[18]：

(9)

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

(10)

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

(11)

2 破冰阻力數值計算模型

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

2.1 數值模型基本假設

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

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

2.2 數值模型離散化方法

以某極地航行船為研究對象，將船體表面三維模型離散為一系列四邊形網格單元。船艏是與海冰發生碰撞的主要區域，且船艏區域復雜，曲率變化大，為了提高計算精度，同時優化計算的效率，針對船艏區域進行網格加密處理，其他區域采用適當的網格大小來劃分。劃分結果如圖3所示。

圖3 船體表面網格劃分結果

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

圖4 破冰阻力數值計算模型

2.3 數值模型動態接觸判別方法

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

(12)

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

(13)

即發生接觸條件，需重新分配該穿透的海冰物質點。

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

2.4 數值模擬破冰阻力計算

圖6 發生穿透物質點的重新定位

具體位置坐標為：

(14)

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

(15)

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

(16)

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

(17)

圖7給出采用接觸判別方法的數值模型模擬結果及局部細節。復雜的船艏區域未發生穿透現象，表明采用接觸判別算法的準確性。

圖7 采用接觸判別的計算結果

3 算例分析

3.1 數值模型參數設定

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

圖8 層冰阻力試驗場景

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

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

參照模型比尺的海冰參數設置如表2所示。

表2 海冰參數原型及模型值

參考文獻[22]收斂性分析計算結果，本文設定時間步長為0.0 005 s，海冰物質點直徑為L/900.0。

3.2 破冰阻力數值模擬現象

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

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

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

3.3 數值模擬破冰阻力計算結果

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

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

圖10 破冰阻力計算數據使用段

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

3.4 不同航速破冰阻力數值模擬

在此基礎下，進行不同航速的數值模擬，不同航速下破冰過程的模擬現象如圖12所示。

圖11 破冰阻力實時曲線

圖12 不同航速下破冰阻力模擬現象

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

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

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

4 結論

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

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