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

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

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

0 引 言

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

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

1 破冰船連續(xù)破冰過程的數(shù)值模擬

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

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

在數(shù)值模擬過程中，假設：

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

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

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

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

1.1 層冰破壞模式分析

1.1.1 確定接觸面積

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

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

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

基于上述理想化破冰過程，根據(jù)冰厚與破冰半徑的大小關(guān)系，接觸面積可以分為三角形接觸和四邊形接觸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的關(guān)系式如下：

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

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

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

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

1.1.3 彎曲破壞分析

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

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

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

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

1.2 浸沒力計算及修正

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

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

當船速增加時，碎冰與船體接觸碰撞更加頻繁，與船體產(chǎn)生的摩擦力增大，浸沒阻力增加?？紤]浸沒阻力與船速之間的關(guān)系，修正后的浸沒阻力為：

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

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

2 經(jīng)驗公式計算方法

2.1 Spencer算法

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

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

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

2.2 Lindqvist算法

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

2.3 Riska算法

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

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

表2 Riska算法中經(jīng)驗系數(shù)取值Tab. 2 Constants in Riska formulation for ice resistance.

3 數(shù)值模擬計算結(jié)果驗證及參數(shù)分析

3.1 建立離散點模型

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

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

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

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

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

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

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

3.3 數(shù)值模擬方法與經(jīng)驗算法對比驗證

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

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

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

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

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

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

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

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

圖6 冰阻力經(jīng)驗公式值與數(shù)值模擬結(jié)果的比較Fig. 6 Comparison of ice resistance obtained by empirical formulas and numerical simulation.

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

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

圖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）可以看出，～的二次方系數(shù)的差別很小，曲線的開口寬窄基本一致，也就是說，速度從1kn增加到5kn，冰阻力隨冰厚的增加趨勢基本一致。由此可以看出，冰厚對冰阻力的影響較為顯著。

4 結(jié) 語

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

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

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

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

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