二維離散方法破冰船層冰破冰阻力數值預報

葛 媛，張佳寧，金 強，周陳炎，張 雷

(1. 大連海事大學 船舶與海洋工程學院，遼寧 大連 116026；2. 南通理工學院 電氣與能源工程學院，江蘇 南通 226000)

0 引 言

由于破冰過程的特殊性及層冰物理力學性質的復雜性，準確預測破冰船在層冰中的冰阻力十分困難。早期，Spencer[1]，Lindqvist[2]和Riska[3]等根據層冰與船舶相互作用的幾個階段，將總阻力劃分成不同的部分，結合大規模的模型試驗和實船數據，提出相應算法對層冰破冰阻力進行了研究。為了模擬船-冰相互作用，Wang[4]通過考慮破碎，彎曲及碎冰形成3個連續過程，基于力學冰破壞模型，采用幾何網格法對冰錐相互作用進行了數值模擬。Nguyen[5]將這一推導出的冰破壞模型應用于船-冰相互作用，模擬動態定位（DP）船在層冰下的運動。Su[6]進一步考慮了冰阻力與船舶運動的耦合問題，優化船-冰接觸模型，對連續破冰模式進行數值模擬。

本文根據簡化后的連續破冰過程，運用編程語言，模擬理想狀態下破冰船在層冰條件下的運動過程，計算破碎力，建立完整的破碎力算法流程圖。由于簡化后的連續破冰過程假設層冰脫落后馬上消失，并不考慮層冰破碎后碎冰在船首的運動，因此計算出的阻力值較實際值偏低。本文在理想假設的基礎上，考慮船首部分碎冰運動引起的浸沒阻力，對破碎力進行修正，計算修正后破冰阻力，得到破冰船在破冰過程中的冰阻力變化曲線。同時對不同工況下的船-冰運動進行數值模擬，分析船速、冰厚等因素對破冰阻力的影響。

1 破冰船連續破冰過程的數值模擬

根據Valanto[7]對連續破冰過程中船-冰相互作用的分析，當船首與層冰發生接觸時，接觸區域的層冰自由邊緣發生局部擠壓破碎，此時層冰破壞模式主要是擠壓破壞和剪切破壞。隨著船的前進，船首與層冰接觸面積增大，破碎力增加，層冰發生彎曲傾斜，當層冰內部的應力超過其應力極限，層冰發生彎曲破壞。發生彎曲破壞的層冰從冰層上斷裂下來，在船體的作用下進一步加速并旋轉，直至失去與船舶的接觸。船首與新的層冰邊緣接觸，開始新一輪循環，如圖1所示。

圖1 連續破冰過程的4個階段Fig. 1 Four stages of continuous ice-breaking process

在數值模擬過程中，假設：

1）船-冰相互作用是一個連續的過程，包括層冰的破碎，彎曲，斷裂的循環；

2）船-冰相互作用僅發生在船體水線面處，且忽略船在破冰運動過程中的垂直運動；

3）在擠壓斷裂過程中，接觸區域平面保持平整；

4）在發生彎曲斷裂時，冰面彎曲斷裂的形狀認為是圓形的。

1.1 層冰破壞模式分析

1.1.1 確定接觸面積

在數值模擬過程中，船體與層冰的接觸情況是判斷破冰過程破碎力大小的關鍵。層冰在發生彎曲破壞后，會產生環向裂縫和徑向裂縫。其中環向裂縫與層冰邊緣的交點是判斷彎曲破壞過程中冰層破損形狀的關鍵位置，在數值模擬中，環向裂縫被理想化為破碎半徑R的函數[1]：

式中：Cl和Cv為經驗參數；Cl為破冰半徑與特征長度的比例；Cv為浮冰的大小隨碰撞速度的變化而變化。這里，Cl為正值，Cv為負值。在程序中，Cl、Cv的值是可調控的，取值情況可根據實驗進行調整。vnrel為船體相對于層冰離散點的垂直破冰速度，l為層冰特征長度:

式中：E為彈性模量；hi為冰厚；ν為泊松比；ρw為海水密度；g為重力加速度。

基于上述理想化破冰過程，根據冰厚與破冰半徑的大小關系，接觸面積可以分為三角形接觸和四邊形接觸2種情況，如圖2所示。

圖2 船冰接觸面積的2種計算模式Fig. 2 Two cases for the calculation of contact area

2種接觸模式下的接觸面積計算公式如下：

情況1當時：

情況2當時：

式中：AC為接觸面積；Lh為冰層上表面的擠壓寬度；Ld為冰層上表面的擠壓深度。

1.1.2 計算破碎力Fcr及垂直分量力FV

破碎力Fcr與接觸面積AC的關系式如下：

式中：σC為層冰抗壓強度；AC為接觸面積。

破碎力Fcr分為水平方向分力FH和垂直方向分力FV。若考慮摩擦阻力[6]，則摩擦阻力也分為水平方向分力fH和垂直方向分力fV，fH與相對速度分量vtrel成正比，fV與相對速度分量成正比，如圖3所示。

圖3 力和速度分量Fig. 3 Force and velocity components.

由于在模擬過程中忽略船體的垂直運動，船體的摩擦力可根據船-冰相對速度來確定。其中，摩擦力分力與破碎力分力如下：式中：μi為層冰摩擦系數；和是船體接觸點處相對速度在法向和切向的分量，其中：

1.1.3 彎曲破壞分析

為了判斷層冰的彎曲情況，Kerr[8]引出失效載荷Pf的概念，即冰層的承載能力。

式中：θ為冰楔的開角，σf為層冰的彎曲強度，hi為層冰厚度，Cf為經驗系數，其變化會影響船速與阻力之間的關系，取值情況可根據實驗進行調整。在數值模擬過程中，通過比較破碎力在垂直方向上的分力與失效載荷的大小，判斷層冰是否發生彎曲失效。當垂直分量力FV＜Pf時，層冰僅在冰緣處發生擠壓破壞；當垂直分量力FV≥Pf時，層冰發生彎曲破壞.

破碎力在數值模擬中的計算流程圖如圖4所示。

圖4 破碎力計算流程圖Fig. 4 Flow chart of ice breaking force calculation

1.2 浸沒力計算及修正

浸沒力是指冰層發生彎曲失效后，斷裂的碎冰在船首處翻轉滑移所產生的阻力。在數值模擬中，假定船首發生彎曲破壞的滑冰脫落后消失，但在實際情況中，船首的碎冰會被淹沒并隨著船體滑行。根據Zhou[9]的模型試驗觀測，船舶模型底部幾乎沒有浮冰滑動，大部分浮冰以較低的漂移速度在冰緣附近橫向移動。根據Lindqvist算法對破冰阻力的劃分，考慮浮冰勢能的損失和船體與浮冰之間的摩擦之后可得浸沒阻力為[10]：

式中：ρw為海水密度；ρi為層冰密度；g為重力加速度；hi為冰厚；B為船寬；T為船長；μ為摩擦系數；Af為船首面積。

當船速增加時，碎冰與船體接觸碰撞更加頻繁，與船體產生的摩擦力增大，浸沒阻力增加。考慮浸沒阻力與船速之間的關系，修正后的浸沒阻力為：

式中：vrel為船體與層冰之間的相對速度，L為船長。

理想狀態下的數值模型不考慮碎冰在船首處的運動，在原有數值模型上增加對浸沒力計算，對破碎力進行修正，得到破冰力。修正后的計算模型的模擬結果在理論上比原有的數值模型更加精確。

2 經驗公式計算方法

2.1 Spencer算法

Spencer算法將總阻力RT分為敞水阻力ROW，冰浮阻力RB，冰清阻力RC和破冰阻力RBR，即RT=ROW+RB+RC+RBR。RB，RC，RBR即公式中系數為 CB，CC，CBR的部分，2010年，Jeong[11]對公式中的無量綱系數進行了修正。修正后的公式：

式中：CB為冰浮阻力系數；CC為冰清阻力系數；CBR為破冰阻力系數；α為作用力指數，β為傅汝德數指數，取值如表1所示。V為破冰船航行速度；B為船寬；T為吃水；ρi為海冰密度；Δρ為水冰密度之差；hi為冰厚；σf為冰的彎曲強度；g為重力加速度。

表1 Spencer算法中無量綱系數取值Tab. 1 Constants in Spencer formulation for ice resistance.

2.2 Lindqvist算法

Lindqvist[2]算法將總阻力分為3個部分：破碎阻力RC，彎曲阻力RB和浸沒阻力RS。各分力及總阻力表達式如下：

2.3 Riska算法

Riska算法[3]是根據波羅的海的一些全尺度試驗數據得出的，主要公式如下：

式中：L為垂線間長；B為船寬；T為吃水；hi為冰厚；Lpar為平行中體的長度；Lbow為水線處的船首長度；V為航行速度；?為艏傾角；f1，f2，f3，f4，g1，g2，g3為經驗系數，取值如表2所示。

表2 Riska算法中經驗系數取值Tab. 2 Constants in Riska formulation for ice resistance.

3 數值模擬計算結果驗證及參數分析

3.1 建立離散點模型

根據破冰船Icebreaker Research Vessel的船體型值，對船體水線進行離散化。由于破冰時船體與層冰的作用位置主要為船首部位，所以編程時對船體首部的離散點進行3次樣條插值，降低離散點之間的間隔，保證計算的精度。層冰理想化為無限大的平面，將平面離散化為無數個間隔相同的離散點，層冰離散點間隔為0.5 m。

圖5為船體水線與冰緣的接觸示意圖，定義向量xv∈RNv×2為船體節點的x，y位置，定義向量xi∈RNi×2為冰邊節點的x，y位置，當海冰受到垂直方向力FV大于海冰彎曲失效載荷Pf時，破碎半徑為R的層冰部分冰節點失效，圖5（b）中虛線部分即失效的冰網格節點，冰緣節點隨之更新。在任意時刻步長ti，檢測第j個冰節點與第k個船節點Djk之間的距離，檢測冰節點是否進入船體內部，判定船體與層冰邊緣是否有接觸，當冰節點與船體發生接觸，計算冰網格與船體型線2個交叉點之間的距離Lh及冰節點尖端到船-冰接觸點連線的垂直距離Ld，根據式（3）和式（4）得到接觸面積AC。

圖5 船體水線與冰緣的離散化Fig. 5 Discretization of ship waterline and ice edges

3.2 船體主尺度及層冰物理特性

在數值模型中，冰區長為220 m，寬為60 m，時間步長為0.05 s，Cf取值為3.1。破冰船主尺度及層冰物理力學性質參數如表3和表4所示。

表3 破冰船主尺度Tab. 3 Primary dimensions of Icebreaker Research Vessel

表4 層冰力學特性Tab. 4 Ice mechanical properties

3.3 數值模擬方法與經驗算法對比驗證

本文以破冰船Icebreaker Research Vessel為實例，計算3～5 kn船速、0.8～1.5 m冰厚下的冰阻力。以破冰阻力曲線對破冰時間t的平均值表示破冰阻力，并與3種經驗算法進行對比驗證。

表5 12種不同船速和冰厚下的算例Tab. 5 Case of different velocity and ice thickness

12種算例下冰阻力變化曲線圖如圖6所示。

其中，simulation1為修正后的數值模擬結果，simulation2為修正前的數值模擬結果。Lindqvist算法出現較早，雖然考慮了船體形狀及層冰力學性質對冰阻力的影響，但是由于理論依據不足，且試驗時采用的船型較小，所以對冰阻力的預估一般偏低。

由圖6數據及曲線走勢可以看出，修正前的數值模擬結果simulation2與經驗公式結果趨勢相同，但結果總體偏低。當船速為3 kn、冰厚為0.8 m時，修正后的模擬結果比修正前高10%左右；當船速為5 kn、冰厚為1.5 m時，修正后的模擬結果比修正前高出18%左右；當船速為5 kn、冰厚達到2 m時，修正后的模擬結果比修正前高出31%左右。考慮浸沒阻力后，數值模擬結果與經驗公式結果更為吻合，說明修正后的數值模型更加精確，模擬結果更加可靠。

3.4 冰厚和船速對冰阻力的影響分析

在數值模擬中，本文以1.0 m，1.2 m冰厚，3 kn，4 kn航速為例，比較分析破冰船在破冰時的冰阻力曲線圖，如圖7所示。

由圖7可知，船冰接觸力隨時間變化并沒有周期變化的規律，阻力值的不規則性是由于破冰模式性質的變化造成的。有時是破碎為主，有時是彎曲。當部分冰層被船體首部連續壓碎而不發生彎曲破壞時，破碎力會急劇增大。

圖6 冰阻力經驗公式值與數值模擬結果的比較Fig. 6 Comparison of ice resistance obtained by empirical formulas and numerical simulation.

當船剛開始進入冰層時，阻力是逐漸增加的，當船體首部完全進入冰層之后，雖然阻力值依然變化劇烈，但是總體均值趨于穩定。船體和層冰接觸時會產生很大的碰撞力，當層冰彎曲失效時，船舶會經歷一個短暫的“卸載”過程，在數值模型中，船首處只受浮碎冰產生的浸沒阻力，所以破冰阻力的曲線變化較為劇烈。

隨著船速和冰厚的增加，冰阻力總體均值呈增加趨勢。在船速增加時，冰阻力的總體均值增加，但曲線變化趨勢較為平緩；在冰厚增加時，冰阻力曲線的振蕩幅度增加，峰值的變化幅度也較大，說明在數值模擬結果中，冰阻力對冰厚變化的敏感度較高。

圖7 破冰阻力時歷曲線及平均值Fig. 7 Time history of ice resistance and mean values

圖8為不同船速及冰厚下，冰阻力的變化曲線，其中破冰阻力取破冰阻力曲線對破冰時間t的平均值。對圖中曲線進行多項式擬合，可得到對應的擬合方程式。

圖8 不同船速及冰厚下破冰阻力變化曲線圖Fig. 8 Ice resistance at different Ice thickness and speeds

其中：x為冰厚hi，y1-5為破冰阻力F在速度1～5 kn時的破冰阻力；x’為航速V，y’1-5為破冰阻力F在冰厚在0.5 m，0.8 m，1.0 m，1.2 m和1.5 m時的破冰阻力。

從圖8（a）可以看出，在冰厚超過0.8 m時，曲線斜率減小，平均變化率增加，冰厚從0.5 m增加至0.8 m時，冰阻力增加9.83%，冰厚從1.2 m增加至1.5 m時，冰阻力增加23.88%，當冰厚大于1.5 m時，冰阻力增加可達到50%。對式（24）來說，隨著冰厚的增加，曲線斜率增大，所以，隨著冰厚的增加，冰阻力變化速率增大。對圖8（b）來說，船速的增加雖然也會引起冰阻力的增加，但影響效果較小，隨著船速增加，冰阻力的增加約為6%～9%。由式（23）可以看出，～的二次方系數的差別很小，曲線的開口寬窄基本一致，也就是說，速度從1kn增加到5kn，冰阻力隨冰厚的增加趨勢基本一致。由此可以看出，冰厚對冰阻力的影響較為顯著。

4 結 語

本文根據連續破冰過程，運用編程語言建立理想狀態下的船-冰運動模型，對層冰冰況下破冰船的破冰運動進行數值模擬，并建立完整的算法流程圖。同時，考慮理想狀態下運動模型忽略船首碎冰這一現象，增加對浸沒阻力的計算，對數值模擬的破冰阻力進行修正，并簡單探討了航速和冰厚對破冰阻力的影響，可以得出結論：

1）增加浸沒阻力修正后的數值模型模擬結果與經驗公式計算結果更加吻合，對冰阻力的預報更加準確、可靠。

2）在冰厚、層冰彎曲強度及層冰摩擦系數等參數不變的情況下，航速的增加會引起冰阻力的增大，但增加趨勢較小，船速增加1 kn，冰阻力的增加約為6%-9%。當船速變化量Δv＜ 1 kn時，冰阻力的變化量Δy不會發生明顯變化。

3）相比于船速，冰厚變化所引起的冰阻力的變化跨度較大。隨著冰厚增加，冰阻力時歷曲線振蕩幅度增加，變化速率增大。當冰厚大于1.2 m時，冰阻力增加速率達到23%以上。

由于在模型中沒有考慮冰厚對層冰抗彎強度等性質的影響，模擬結果可能略有偏差，在接下來的研究計算中會進一步進行修正。