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

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

0 引 言

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

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

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

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

1.1 基本原理

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

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

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

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

式中，Es為結(jié)構(gòu)物的彈性模量，Ei為冰的彈性模量，v為粘性系數(shù)。冰單元的法向接觸剛度Kn可表示為

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

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

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

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

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

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

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

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

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

1.2 船艏及冰層離散化

1.2.1 船艏離散化

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

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

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

1.2.2 冰層離散化

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

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

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

1.2.3 計算參數(shù)輸入

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

表2 船艏模型破冰阻力數(shù)值計算輸入?yún)?shù)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

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

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

1.3 計算結(jié)果分析

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

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

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

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

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

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

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

2 船艏破冰阻力數(shù)值計算結(jié)果驗證

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

2.1 破冰效果比較分析

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

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

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

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

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

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

利用Bri-DEM 程序?qū)Ρ鶎拥钠茐募八楸倪\動進行模擬，能夠較好地反映出冰層的破壞過程，并可模擬出冰層首先破壞的位置及其對周圍冰層的影響。

2.2 破冰阻力對比分析

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

圖12 數(shù)值計算結(jié)果與試驗數(shù)據(jù)對比Fig.12 Comparison of numerical calculation results and test data

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

3 多方案破冰阻力數(shù)值預(yù)報

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

3.1 船艏傾角的影響

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

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

3.2 冰厚的影響

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

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

4 結(jié) 論

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

（1）基于離散元方法對船艏破冰阻力進行了數(shù)值計算，利用矩形單元對冰層進行離散，采用元胞數(shù)組對冰單元進行儲存，控制冰單元力的傳遞方式，使每一個冰單元在循環(huán)中均能夠?qū)崟r地更新狀態(tài)，提高了計算效率；

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

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

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