冰區船艏部連續破冰模式下破冰阻力數值計算研究

剛旭皓，田于逵，季少鵬，寇 瑩，國 威，余朝歌

（中國船舶科學研究中心，江蘇 無錫 214082）

0 引 言

冰區航行船舶是北極探索和開發過程中不可或缺的基本裝備。冰區船舶冰阻力模擬既能為冰區航行船舶的初步設計提供參考，又能為船-冰相互作用研究提供依據，避免船舶在航行過程中因主機功率不足而隨流冰漂移或傾覆，因此準確模擬船舶在冰區航行時的破冰過程并預報船舶冰阻力具有重要的理論和現實意義。

隨著計算技術的不斷進步，數值計算方法逐漸成為冰區船舶冰阻力性能研究的有效工具。Wang[1]將破冰過程理想化為連續的接觸-擠壓-彎曲破壞的循環，提出了時域內數值流程的解決方法；Sawamur[2]提出了一種用于循環計算冰阻力的連續接觸流程，模擬了碎冰受到外力作用后的動態響應；Su[3]將破冰的過程分為破冰和碎冰運動兩個部分，建立了冰載荷與船體運動的三自由度耦合方程，并將計算結果與破冰船AHTS/IB 試驗數據對比，顯示了很好的一致性。海冰的離散元模型從20世紀80年代發展起來并不斷完善，并在涉冰研究領域取得了一定的研究成果。Hopkins[4]等將冰塊離散為隨機尺寸的圓盤，模擬了海冰在斜坡式結構前的堆積過程和壓力冰脊的形成；Lau[5]在此基礎上，對平整冰與錐體結構相互作用中的非線性大變形和斷裂過程進行模擬，編寫了三維離散元程序DEMICE 3D；季順迎和狄少丞等[6-7]將基于GPU的高性能計算方法引入離散元數值模擬中，顯著提高了DEM方法計算規模和計算效率,并將其應用于破冰船的研究。除了上述數值計算方法之外，理論分析方法和邊界元方法也在船-冰相互作用領域有一定的應用，但由于船-冰相互作用過程的復雜性和隨機性較強，還需進一步完善。國內冰阻力預報數值研究已經取得了較大的進步，但由于模型試驗和實船驗證不充分，預報的準確性和適用度相對不足。只有將冰阻力預報與海冰物理力學特性及試驗驗證緊密結合，才能建立起有效的冰阻力預報方法，為冰區船航行性能與船型設計優化等工程應用提供扎實、有力的技術手段。

本文基于中國船舶科學研究中心小型冰水池（CSSRC SIMB）中開展的連續破冰模式下冰區船艏破冰阻力模型試驗[8]，根據離散元模型的基本原理，結合元胞數組方法，利用冰力學試驗數據精細化參數，建立離散元數值計算模型Bri-DEM，對船艏模型的破冰阻力及冰層破壞模式進行數值計算，分析船艏傾角和冰厚變化對破冰阻力的影響，獲得破冰阻力相應的變化規律。

1 基于離散元方法的船艏破冰阻力數值計算

1.1 基本原理

在離散元方法中，海冰單元的力學行為是牛頓第二定律與接觸處力-位移定律交替應用的結果：牛頓第二定律用于確定各質點在接觸力、阻尼力和物體力作用下的平動和旋轉，而力-位移定律則用于更新各接觸點相對運動產生的接觸力以及單元的位置信息。為模擬海冰的連續特性，相鄰單元之間采用平行粘結模型粘連。平行粘結模型的基本原理是設想一組具有恒定法向和切向剛度的線性彈簧，以兩單元的接觸點為中心，均勻分布在粘結圓盤上，完成單元間作用力和力矩的傳遞。單元之間的粘結失效通過梁理論進行求解：

式中，Fn、Fs分別為作用在粘結圓盤上的法向力和切向力，Mn、Ms分別為作用在圓盤上的法向力矩和切向力矩，σmax、τmax分別為作用在粘結圓盤上的最大的拉應力和剪應力，I、J分別為粘結圓盤的慣性矩和極慣性矩，A為粘結圓盤面積。當圓盤上的最大拉伸應力σmax超過法向粘結強度，或最大剪切應力τmax超過切向粘結強度時，顆粒間的粘結鍵將會被破壞，由此模擬內部裂紋產生的情況[9-10]。

在單元粘結處，單元受力和位移的關系可以通過以下參數來表示：法向剛度Kn和切向剛度Ks、摩擦系數μ、法向阻尼和切向阻尼，以及粒子重疊時在重疊區域中心形成的接觸區。

單元的法向接觸剛度和切向接觸剛度對單元間的受力具有重要的影響，單元的接觸剛度與接觸處的有效彈性模量相關。在冰單元與結構物相互作用時，單元與結構之間的有效彈性模量Ee可表示為

式中，Es為結構物的彈性模量，Ei為冰的彈性模量，v為粘性系數。冰單元的法向接觸剛度Kn可表示為

式中，hi、ρi表示冰厚和冰密度，l為冰特征長度，g為重力加速度。根據Ji[11]的研究，Ks與Kn之間存在線性關系，切向剛度Ks可表示為

式中，r根據彈性模量和切變模量的關系確定，取r=1/2（1+v）。

對接觸處進行力學分析，接觸力Fi可分解為切向力和法向力：

式中，Fn、Fs分別表示接觸力的法向和切向分量，ni、ti為單位向量。根據Hertz接觸理論，法向力Fn可以表達為

式中，xn表示單元之間的重疊量，vn為法向相對速度，Cn為單元間的法向阻尼系數，可表示為

式中，ξn為單元間的阻尼比，M為單元的等效質量，ei為回彈系數。ei可參考Ramírez[12]提出的不同碰撞速度下的回彈系數計算，如圖1所示。

圖1 不同碰撞速度下的回彈系數[11]Fig.1 Rebound coefficients at different speeds

單元間的切向力Fs一般以增量的形式進行計算。剪應力的增量為每個時間步內切向方向的變化量與切向剛度的乘積：

式中，Δxs和vs分別為兩個單元在接觸點處的切向位移增量和相對切向速度，Cs為切向阻尼系數。假設單元間的接觸力滿足摩爾-庫倫定律，則有

1.2 船艏及冰層離散化

1.2.1 船艏離散化

設計破冰型冰區船艏部模型，利用三維建模軟件UG進行數值化建模，利用網格劃分工具Hypermesh 將模型合理地劃分為一系列三角離散單元，圖2 是船艏模型離散化示意圖，表1 為船艏模型基本參數。提取網格節點信息導入數值計算模型，通過判斷冰單元與船艏模型網格單元之間的接觸來判斷冰單元受到的接觸力和力矩，以此對船艏受到的作用力進行求解，進而計算破冰阻力。

表1 船艏模型基本參數Tab.1 Basic parameters of bow model

圖2 船艏模型離散化示意圖Fig.2 Schematic diagram of bow model discretization

1.2.2 冰層離散化

利用矩形單元采用最密排列的方式對冰層進行構造。對冰層進行離散化后，將單元編號為k1、k2，…，kn，采用元胞數組儲存顆粒的位置和運動信息，隨著時間步更新，利用元胞數組特性對粒子的位置和運動信息進行實時更新，可實時對粒子間的相互作用力進行計算，簡化了顆粒搜索過程，同時也利于對冰單元kn與船艏相互作用力以及單元與單元之間相互作用力的計算。

冰單元kn與船艏模型相互作用時，作用力主要有單元與船艏之間的擠壓力以及單元在船艏表面滑移產生的阻尼力以及摩擦力。對冰單元之間力的傳遞方式進行設置，冰單元與船艏接觸的位置不同，力在冰單元之間的傳遞也會有所不同，圖3表示不同接觸位置接觸力的傳遞方式。當船艏與冰層左側的冰單元接觸時（k1、k2），冰單元之間的接觸力自左側向右傳遞和向后傳遞；當船艏與冰層右側的冰單元接觸時（k4、k5），接觸力自右側向左傳遞和向后傳遞；當船艏與中部的冰單元接觸時（k3），冰單元之間的接觸力自中間向兩側傳遞及向后傳遞。在完成第一步傳遞之后，接觸力逐漸拓展，直到覆蓋到整個冰層，計算每一個冰單元受到的力。通過對接觸力的不同傳遞方式進行模擬，以此來計算冰單元之間力的作用以及冰層的裂紋生成和斷裂。

圖3 不同接觸位置接觸力的傳遞方式Fig.3 Transmission mode of contact force at different contact positions

1.2.3 計算參數輸入

考慮到冰的力學特性對冰單元之間的粘結作用有較大的影響，為了確定計算參數，對模型冰的力學特性在不同加載速率下的變化進行了研究。研究發現：冰的彎曲強度隨加載速率的增加基本沒有發生變化；而冰的壓縮強度，在試驗速度范圍內呈現出隨著加載速率增加而逐漸減小。圖4為彎曲強度與壓縮強度隨加載速率的變化。為避免偶然性結果，對同一工況進行了重復性實驗，圖中使用不同形狀標記表示。根據不同速度下冰強度值確定單元間的粘結強度，將冰的物理力學參數和計算得到的粘結強度輸入計算模型，具體的計算參數如表2所列。

表2 船艏模型破冰阻力數值計算輸入參數Tab.2 Input parameters for ice-breaking resistance calculation of bow model

圖4 彎曲強度與壓縮強度隨加載速率的變化Fig.4 Variation of bending strength and compression strength with loading rate

本次數值計算冰厚取40 mm，船艏模型與冰相互作用的速度設置為0.01 m/s、0.05 m/s、0.10 m/s、0.15 m/s 和0.20 m/s，冰層的邊界條件設置為三邊約束、一邊自由的狀態，自由邊為船艏模型進入冰層的方向，與冰水池模型試驗的情況一致。在數值計算模型中，將船艏模型的網格模型導入計算程序中，設置船艏與冰層的相對位置，如圖5所示。

圖5 三維船艏與冰層的設置形式Fig.5 3D setting of the bow and ice sheet

1.3 計算結果分析

利用Bri-DEM 程序對船艏與冰相互作用過程進行模擬，獲得了船模破冰阻力，同時對破冰過程及碎冰運動進行了模擬，圖6為不同時刻冰單元的運動情況。

圖6 不同時刻冰單元的運動情況Fig.6 Movement of ice units at different times

破冰過程中，船體艏柱線首先與冰層接觸，對冰單元產生擠壓作用，引起冰單元旋轉，同時冰單元之間的相互作用力逐漸向后、向兩側延展，引起冰單元運動，直至冰單元之間的作用力超過粘結強度，引起冰單元斷裂。船艏與冰層之間的接觸隨著船艏前進，由單點接觸轉變為多點接觸，形成接觸-擠壓-破壞的周期性過程。由于未考慮水動力的影響，冰單元在冰層上斷裂之后，會直接向下運動，且不會與冰層再次粘結。

船艏在冰層中破冰航行時產生接觸-擠壓-彎曲破壞的循環過程，從船艏受力的時歷曲線中可看出，曲線具有明顯的脈動特性和周期性，圖7 為速度0.15 m/s 時船艏受力時歷曲線。冰單元的受力與船艏與冰的重疊距離有關，在單一時間步內，隨著船艏與冰單元之間的重疊面積增加，引起冰單元受力增加，船艏受力增加，但計算得到的重疊距離具有隨機性，因此不同時間步的力值出現較大的波動。

圖7 速度0.15 m/s時船艏受力時歷曲線Fig.7 Time-history curve of bow force at a speed of 0.15 m/s

船舶破冰速度的大小決定了船-冰相互作用的速度和頻率，是影響船體破冰阻力的一個重要因素，分別對不同速度的計算結果繪制力-速度曲線，如圖8所示。隨著破冰速度的增加，破冰阻力呈增大趨勢。為驗證計算模型準確性，將采用模型試驗數據和經驗公式對計算結果進行校驗。

圖8 不同速度下船艏破冰阻力Fig.8 Ice-breaking resistance of ship bow at different speeds

2 船艏破冰阻力數值計算結果驗證

在冰阻力的數值預報中，能否準確地預報冰阻力值以及破冰模式，是評價數值計算模型優劣的重要指標。將Bri-DEM 計算結果與模型試驗數據進行比較，可驗證數值計算模型的可靠性。在此基礎上，使用驗證后的模型進行拓展計算，進一步對船艏冰阻力進行研究。

2.1 破冰效果比較分析

計算模型可將冰單元的運動及位置信息輸出，獲得冰層斷裂和碎冰旋轉浸沒的過程。當船艏與冰層接觸時，船體艏柱與冰層最先產生接觸，冰層產生擠壓破壞，隨著船艏繼續破冰前進，垂向力逐漸增加，迫使冰層向下彎曲破壞，碎冰從冰層上脫落；同時會對其相鄰的冰單元產生影響，使其發生不同程度的變形，這種現象與試驗過程中船艏與冰層接觸產生擠壓破壞并產生裂紋、最終發生彎曲破壞的過程較為一致，見圖9。

圖9 彎曲破壞階段對比Fig.9 Comparison of flexural failure stages

隨著船艏繼續破冰前行，冰單元與船艏的接觸逐漸由單點接觸逐漸轉換為多點接觸，同時冰層內部單元之間的粘結作用會對航道兩側的冰層產生影響，如果航道兩側的冰層受力未達到斷裂強度值，船艏破冰后，將恢復到原來的狀態。與試驗照片進行比對，數值模擬能夠反映冰層的破壞與船艏對破冰航道外冰層的影響，如圖10所示。

圖10 冰層破壞模式對比Fig.10 Comparison of ice destruction patterns

由于在數值計算過程中未考慮水動力以及水的浮力的作用，冰單元斷裂之后會直接脫落，而在模型試驗過程中，碎冰會逐漸覆蓋船底，并隨船向后運動，如圖11所示。

圖11 冰層破冰效果對比Fig.11 Comparison of ice breaking effects

利用Bri-DEM 程序對冰層的破壞及碎冰的運動進行模擬，能夠較好地反映出冰層的破壞過程，并可模擬出冰層首先破壞的位置及其對周圍冰層的影響。

2.2 破冰阻力對比分析

在船艏與冰相互作用的過程中，隨著破冰速度的增加，船艏受到的破冰阻力也在不斷地變化。統計比較Bri-DEM 計算結果與模型試驗結果，如圖12 所示。在與模型試驗比較的同時，利用Jeong 經驗公式估算了相同工況下的破冰阻力。Jeong[13]公式是在Spencer公式的基礎上提出的一種冰阻力經驗估算方法，計算公式表達如下：

圖12 數值計算結果與試驗數據對比Fig.12 Comparison of numerical calculation results and test data

式中：Ri(v)、Ro(v)分別為冰阻力和敞水區域阻力；Cb、Ch和Cr分別為浮力系數、清冰力系數和破冰力系數；Fh和Sn分別為Froude 數和強度因子；p和q分別為Froude 數和強度因子的冪指數。可以看到，Bri-DEM 能夠模擬出破冰阻力的變化趨勢，即隨著速度的增加，破冰阻力逐漸增大，Bri-DEM 計算結果與模型試驗結果、Jeong 經驗公式估算結果大體相近，且Bri-DEM 計算結果與模型試驗結果二者較為一致。

3 多方案破冰阻力數值預報

冰區船在冰層中航行的破冰阻力受到諸多因素的影響，如冰的物理力學特性、船體形狀、船舶航態以及航速等，選取船艏傾角和冰厚作為典型影響因素，研究其對破冰阻力的影響。

3.1 船艏傾角的影響

船艏傾角大小對破冰模式以及破冰時主要受力類型有重要的影響。當船艏傾角較小時，破冰的主要類型為彎曲破壞，而船艏傾角較大時，破冰模式為擠壓破壞和彎曲破壞的共同作用。為了研究船艏傾角改變對冰阻力的影響，計算了不同船艏傾角下的破冰阻力，如圖13 所示。隨著船艏傾角的增加，船艏的破冰方式中擠壓破壞的比例逐漸升高，彎曲破壞的比例下降，船艏受到的破冰阻力逐漸增加。在現代破冰船舶設計中，船艏的外傾角增加和水線角與船艏傾角減小可以提高船舶的破冰性能，本文的模擬結果反映出這種隨船艏傾角減小破冰阻力相應減小的效果。

圖13 不同船艏傾角下破冰阻力比較Fig.13 Comparison of ice-breaking resistance at different bow angles

3.2 冰厚的影響

冰厚是海冰最重要的物理特性之一，不同厚度的冰層具有不同的力學性質，其抵抗破壞的能力也會不同。海冰條件對船舶冰阻力的影響，也以冰厚最為顯著。為了分析冰厚對冰阻力的影響，計算了不同冰厚下的破冰阻力，圖14 為不同冰厚下船艏破冰阻力比較。隨著冰厚的增加，船艏的破冰阻力增加，并且速度越高，這種趨勢越明顯。因此冰區船在不同冰區航行時，應根據所處冰區的冰厚，在考慮效率和推進功耗的基礎上，合理地選擇最佳航行速度。

圖14 不同冰厚下船艏破冰阻力比較Fig.14 Comparison of ice-breaking resistance of bow under different ice thicknesses

4 結 論

冰阻力預報問題一直是冰區航行船舶研究的熱點和難點。本文采用離散元方法編寫了冰阻力數值預報程序Bri-DEM，對船艏破冰阻力進行了數值計算并利用模型試驗結果進行了校驗，最后開展了多方案數值預報與分析，研究了船艏傾角和冰厚對破冰阻力的影響，主要得到了以下結論：

（1）基于離散元方法對船艏破冰阻力進行了數值計算，利用矩形單元對冰層進行離散，采用元胞數組對冰單元進行儲存，控制冰單元力的傳遞方式，使每一個冰單元在循環中均能夠實時地更新狀態，提高了計算效率；

（2）將數值計算結果與模型試驗結果進行比較，船艏破冰阻力呈現出隨船艏破冰速度的增加而增加的趨勢，與試驗結果趨勢一致，同時也基本能夠反映出冰層的破碎以及碎冰的旋轉浸沒的過程。

（3）利用經試驗驗證的數值計算模型研究了典型影響因素船艏傾角和冰厚對船艏破冰阻力的影響，表明隨船艏傾角和冰厚的增加，船艏受到的破冰阻力會逐漸地增加。

通過對冰區船艏模型破冰阻力及破冰模式的預報和分析，建立了破冰阻力預報的數值模型，可為冰區船舶破冰能力分析及船型優化提供技術支撐。