張 媛，王 超，葉禮裕，楊 波

（哈爾濱工程大學，哈爾濱 150001）

0 引 言

極地日益增長的航運活動、科學考察工作、資源開采需求以及軍事領域的應用價值極大提升了極地結構物的需求量，并對其極地航行船舶的性能提出了更高的要求。破冰船是極地航行必需的開辟航道、引領冰區(qū)航行、保證極地船舶安全航行的重要交通工具[1-2]。破冰能力是衡量破冰船性能的關鍵標準之一，而破冰船的船艏形狀特征參數(shù)在很大程度上決定了破冰效率、破冰模式，以及破冰第二階段的冰翻轉過程和浸沒過程的碎冰塊運動軌跡、運動模式和冰的浸沒阻力。為了深度掌握船艏型線特征對破冰船破冰能力的影響，探索各關鍵船艏特征參數(shù)與破冰過程的相關性機理，以幫助設計出型線更優(yōu)、性能更佳的破冰船艏，提高國家極地航行事業(yè)的競爭力，開展破冰船艏形狀特征的破冰模式差異預報分析研究是十分必要的。

破冰船艏部形狀可通過5 個特征參數(shù)來表達：外飄角、水線進流角、縱剖線傾角、首柱傾角、艏部長度[3]。其中，外飄角主要影響破冰船的破冰效率及破碎冰塊的下沉效率，水線進流角的影響主要體現(xiàn)在清除船艏前端及其兩側形成的堆積冰的速率，即清冰效率，縱剖線傾角和首柱傾角也會一定程度地影響到破冰和冰塊下沉效率。破冰船艏的設計主要是圍繞上述特征開展的，目前常見的破冰船艏主要有7種：具有平行縱剖線的直線艏、凹形艏柱艏（懷特艏）、高外飄角艏（米勒艏）、具有車刀的勺形艏、具有脊骨的半勺形艏、平板系艏和狄森-瓦斯艏[4]。其中，前三種破冰艏均保持了平滑的船體，敞水域阻力性能較好，通常被稱為傳統(tǒng)型破冰艏；而后四種因艏部型線非平滑過渡，通常被稱為非傳統(tǒng)型破冰船艏。

數(shù)值預報方法被認為是研究冰-船接觸過程中破冰模式以及破冰載荷的有效手段。在平整冰、冰脊和冰山等冰場中，常采用有限元方法（FEM）模擬局部冰載荷和全部冰載荷[5-6]。然而，在描述冰-船相互作用過程中的破冰特性時，由于離散元方法（DEM）既可以在微觀尺度上描述冰結構的離散性，又可在宏觀尺度上模擬船-冰相互作用過程中的冰破裂過程，DEM 表現(xiàn)出了其在冰荷載預報方面的固有優(yōu)勢[7-10]。文獻[11]采用SPH方法模擬了層冰工況下，冰-船接觸過程中的冰破壞過程并預報了破冰船的破冰阻力。近場動力學方法是另外一種具有模擬材料斷裂特性優(yōu)勢的無網(wǎng)格方法，現(xiàn)已成功應用于冰的模擬預報中[12-14]。

本文的主要目的是應用數(shù)值預報的手段來分析傳統(tǒng)型破冰船艏和非傳統(tǒng)型破冰船艏的破冰模式差異。為了更好地模擬破冰模式，文中采用了一種新型非局部無網(wǎng)格方法——常規(guī)狀態(tài)近場動力學理論（后文簡稱OSB-PD）來建立冰的彈脆性本構模型；建立了一種快速接觸檢測算法（后文簡稱FCDA）來處理基于粒子方法的固固耦合接觸判斷問題；嵌入了消息傳遞接口（MPI）的并行技術來提高無網(wǎng)格方法的計算效率；在FORTRAN 語言環(huán)境下編寫了破冰船破冰過程的數(shù)值預報程序；預報了兩種典型船艏形狀破冰船的破冰過程；最后總結分析了兩種船艏的破冰模式和破冰阻力的差異。

1 冰-船作用數(shù)值模型

1.1 冰的本構模型

文中冰的本構模型是基于OSB-PD 理論建立的各向同性、均質分布的標準彈脆性材料。在OSBPD 理論中，冰材料離散為無限個材料粒子點，通過關注每個粒子點的運動信息和物理信息，進而建立整個冰體的變形和破壞的數(shù)學模型。由于PD方法是非局部方法，因此每個材料點與其相鄰一定范圍Hx內的粒子存在相互作用關系，超出這個范圍的粒子與該粒子不存在任何作用關系。在笛卡爾坐標系下，粒子空間位置為x，其占據(jù)一定的空間體積Vx，其密度用ρ(x)表示，鄰域粒子的坐標為x'。當變形發(fā)生時，粒子發(fā)生位移u，且具有新的坐標位置y，同樣地，鄰域粒子的位置發(fā)生變化為y'。粒子間的相互作用用力密度T來表示，兩個粒子間的力密度方向相反，大小不同，方向分別指向對方粒子。由此可以得知，兩個粒子間的作用力是兩個不同的力密度，分別為[x,t]和[x',t]。OSB-PD 方法的最終控制方程如式（1）所示[15]：

通過應變能密度和力密度關系與經(jīng)典介質力學中對應關系的對照可推導出OSB-PD 方法的相關常量參數(shù)，代入公式（1）中，即可得到PD的詳細積分表達式：

式中，a、d和b為PD 常量，Λ為輔助參數(shù)，δ為近場域的半徑，θ和θ'分別為當前粒子點和當前作用鄰域粒子點的體積膨脹，s為粒子間相互作用的伸長量，b為外部作用力。s和θ的表達式分別為

材料破壞的標準通過伸長量s表示，當s超過材料的極限伸長量時，則材料發(fā)生破壞，這種破壞過程是不可逆的，即粒子間變形伸長量超過極限之后，相互作用會永久消失。由此可以簡單引入一個歷史變形狀態(tài)標量Ω來表示粒子間的作用關系：當Ω取值為0 代表作用關系消失；當Ω取值為1 時代表作用關系仍舊存在。

1.2 接觸檢測算法（FCDA）

在FCDA 方法中，船體離散為一系列足夠表述船體表面形狀的四邊形面元，則冰粒子與船體的碰撞檢測過程可以簡化為空間中點和面的相對位置與距離判斷的數(shù)學問題。檢測步驟如下：

（1）冰材料離散為粒子形式，船體劃分為接近平面的四邊形的面元；

（2）在t時刻，任何一個粒子都有與曲面接觸的可能性，為了實現(xiàn)高效的搜索，可以先排除該時刻完全不可能與曲面接觸的粒子以減少粒子搜索量，因此建立一個包含整個不規(guī)則目標碰撞結構的規(guī)則矩形體，只有進入這個矩形體的粒子才有可能與目標面發(fā)生碰撞，由此排除了大量不可能發(fā)生接觸的粒子，減少了不必要的粒子搜索過程；

（3）在t+1時刻冰粒子穿透矩形體和目標面，接下來只需關注進入矩形體的所有粒子；

（4）在這些穿透矩形體的粒子點里面找到每個粒子點接觸或者穿透的唯一面元；

（5）假設粒子的坐標為(x0,y0,z0)，對于目標撞擊體表面的所有四邊形面元，可以找到其四個點在x方向上的最小值xmin和最大值xmax，在y方向上的最小值ymin和最大值ymax，以及在z方向上的最小值zmin和最大值zmax。若面元與點關系為xmin＜x0＜xmax，且zmin＜z0＜zmax，或者xmin＜x0＜xmax，且ymin＜y0＜ymax，則可認為這些面元可能與該物質點發(fā)生碰撞，該面元為可能與粒子碰撞的元素；

（6）在上一步找到的可能碰撞面元中，開始判斷粒子是否穿透或者接觸面元，此時接觸檢測過程已經(jīng)簡化為空間點面位置關系的數(shù)學問題，應用空間點面距離公式即可，建立面元的空間方程Ax+By+Cz+D=0，則最終判斷粒子是否與面接觸的準則如下：

冰粒子與船體的接觸檢測過程完成后，對進入船體的粒子重新分配位置，并更新粒子的運動信息和計算冰載荷，該過程可直接參考文獻[15]第十章的處理過程，本文不再解釋。

1.3 信息傳遞接口并行技術（MPI）

本文編譯了基于MPI手段的OSB-PD 并行程序，可以在小內存的硬件上實現(xiàn)大數(shù)據(jù)量的運算，也可以為以后的超大數(shù)據(jù)量的并行擴展提供基礎技術保障。

如圖1（a）所示，模型冰層在長寬方向上遠遠大于厚度方向，因此針對冰層所在平面，將計算域平均分解為兩個維度的9 個子域并行進程，所占用線程從0 編號至8。該平整冰層x方向上的粒子數(shù)量為nx，y方向上的粒子數(shù)量為ny，z方向的粒子數(shù)量為nz，劃分在各個處理器后三個方向上的粒子數(shù)分別為、和，且=nz，由此每個處理器中的總粒子數(shù)為n=··nz。PD 方法中每個粒子與其鄰域粒子存在作用關系，鄰域的尺寸為δ=m·dx，dx為粒子間距，m為正整數(shù)，代表倍數(shù)關系。由此并行中每個處理器邊界處的粒子會與相鄰處理器中的粒子存在信息交互傳遞，要傳遞信息的粒子占據(jù)的寬度為m個dx，假設m=3，則可以設定不同相鄰處理器之間的信息交換的計算域寬度為3dx，如圖1（b）所示，則每個處理器最大粒子總容量為ntot=n+nz·2 +·nz·2 +nz·3·3·4。

圖1 并行計算域的劃分示意圖Fig.1 Schematic diagram of the parallel computing domain partition

通過上述的并行策略和計算域的劃分方法，使用MPI 常用指令，即可編譯基于MPI 方法的PD 并行程序。

2 破冰過程數(shù)值預報分析

基于上述建立的數(shù)值計算方法，本章分別建立了傳統(tǒng)型破冰船艏和非傳統(tǒng)型破冰船艏的破冰船連續(xù)破冰的數(shù)值模型。為了探究不同船艏類型在破冰過程中的差別，本文預報了相同工況下兩種破冰船船艏的破冰過程，并對比分析了破冰模式和破冰載荷。

2.1 兩種船艏形狀描述

傳統(tǒng)型的破冰船艏是具有平行縱剖線的直線艏型，該船型起源于1950 年Souiet和Finland 的破冰船，呈尖瘦型的船艏形狀，具有優(yōu)秀的破冰能力，是發(fā)展至今仍舊被廣泛應用和改進的破冰船艏型。非傳統(tǒng)型破冰船艏選用簡化的狄森-瓦斯型艏（Thyssen-Waas Bow），這種船艏和傳統(tǒng)型船艏相比，具有明顯的區(qū)別，該船型主要靠船艏最大寬度位置來造成冰面的剪切破壞進而依靠重力造成冰面的彎曲破壞，該船型具有極好的冰清除能力，可以開辟較寬且碎冰較少的破冰航道。兩種船型的模型圖如圖2所示，詳細的船型參數(shù)如表1所示。本文中兩種船型的計算縮尺比均為25。

表1 兩種船型的主尺度參數(shù)Tab.1 Principal dimensions of two kinds of icebreakers

圖2 兩種不同船艏形狀破冰船模型圖Fig.2 Icebreaker models with two different bow shapes

2.2 計算輸入

冰參數(shù)、離散信息和工況設定參見表2，該輸入信息對應于模型的縮尺比25。例如冰的實尺度厚度為1.0 m，則計算中輸入為0.04 m。參量的換算關系參考ITTC冰試驗的換算規(guī)則[16]。

表2 計算模型參數(shù)Tab.2 Calculation parameters of the model

2.3 計算結果

2.3.1 數(shù)值方法驗證

本小節(jié)首先對比了數(shù)值計算結果和試驗結果的破冰模式、數(shù)值計算結果的破冰載荷與試驗[17]和經(jīng)驗公式[18]的破冰載荷來驗證數(shù)值方法的可靠性，該實驗中的破冰船采用本文中傳統(tǒng)型船艏破冰船。破冰模式的結果對比如圖3 所示。為了和經(jīng)驗公式的結果做對比，破冰載荷依據(jù)ITTC 中冰區(qū)船舶換算公式換算為實尺度的冰載荷數(shù)據(jù)，對比圖如圖4所示。

圖3 破冰模式驗證示意圖（傳統(tǒng)破冰船艏；實船4節(jié)航速）Fig.3 Icebreaking pattern verification comparison snapshot(Conventional bow with prototype;speed:4 kn)

圖4 破冰阻力驗證對比圖（傳統(tǒng)破冰船艏；實船4 kn航速）Fig.4 Icebreaking resistance verification comparison(Conventional bow with prototype;speed:4 kn)

從破冰模式來看，環(huán)形裂紋的生成和擴展引發(fā)了冰層的彎曲破壞模式，整個破冰過程的破壞模式主要都體現(xiàn)為彎曲破壞模式。每一次冰層破壞模型的循環(huán)基本為：首先產(chǎn)生沿船長方向的環(huán)形裂紋并向船艏兩側不斷擴展，該類型裂紋總是從船艏的半船寬處開始成型；隨著船舶的前進，環(huán)形裂紋向部和船艉擴展，同時在船艏半船寬處產(chǎn)生放射狀的裂紋；緊接著，基本平行于船艏側邊的二次環(huán)狀裂紋產(chǎn)生并且擴展，伴隨著局部區(qū)域的冰擠壓破碎；船艏產(chǎn)生較短的環(huán)狀裂紋。數(shù)值模擬結果很好地捕捉到了上述在試驗中觀測到的破冰現(xiàn)象。

從破冰阻力均值來看，數(shù)值計算破冰力的均值為0.969 MN，試驗破冰力的均值為1.124 MN，經(jīng)驗公式計算得到的均值力為1.117 5 MN。由此可以得到結論：數(shù)值計算結果與試驗結果和經(jīng)驗公式的結果具有很高的吻合度，本文建立的數(shù)值計算方法能夠高效地模擬破冰過程和準確地獲取破冰阻力。

2.3.2 兩種船艏的破冰模式

本小節(jié)計算了兩種破冰船艏在實船航速為3 kn 的工況下的破冰模式和破冰載荷，分別如圖5和圖6所示。

從圖5 可以發(fā)現(xiàn)，兩種不同船艏的破冰船破冰過程表現(xiàn)出明顯不同的破冰模式。在裂紋初始階段，傳統(tǒng)尖瘦型的破冰船艏首先對冰層造成局部擠壓，形成船艏形狀的破冰缺口，進而在船艏兩側半寬處產(chǎn)生沿船長方向的環(huán)狀裂紋，向艉部和部演變發(fā)展，此時的破冰模式體現(xiàn)為船艏輪廓處的擠壓破壞和平行于船艏外側輪廓線的彎曲破壞，船艏兩側形成少許碎裂的冰塊；而非傳統(tǒng)型船艏在冰船接觸位置發(fā)生極小區(qū)域的擠壓破壞，形成輕微的類似于船艏形狀的弧狀冰層缺口，當破冰船進一步穿過冰層時，由于船艏沒有明顯的艏柱沖破冰層，而是平緩的過渡型船艏壓向冰面，形成了典型的放射狀裂紋和穿越放射裂紋端部的圓弧狀裂紋。在裂紋擴展演變階段，傳統(tǒng)型船艏破冰模式和冰層破壞的模式循環(huán)過程與2.3.1小節(jié)介紹的相同；而非傳統(tǒng)型船艏的破冰模式簡單地由放射狀裂紋和穿越放射裂紋的圓弧狀裂紋組成，冰層破壞模式的循環(huán)過程為放射狀裂紋和圓弧狀裂紋有規(guī)律地層層擴展，且每次破壞循環(huán)周期內的圓弧狀裂紋半徑都比上一個循環(huán)大。總體來看，兩種破冰模式中，彎曲破壞主導了整個破冰過程，冰層破壞的成型形狀基本與船艏輪廓相一致。

圖5 兩種類型船艏的破冰模式對照Fig.5 Comparison of icebreaking patterns of two bow types

從實尺度破冰航道來看，傳統(tǒng)尖瘦型破冰船的型寬為23.0 m，開辟的航道寬度為28.56 m。非傳統(tǒng)型船的型寬為20.0 m，破冰寬度為27.3 m。兩種船型的破冰航道寬度與船寬的關系分別為1.25B和1.36B。因此，對比這兩種破冰船艏可發(fā)現(xiàn)，非傳統(tǒng)型船艏破開航道寬度的能力更顯著。另外，傳統(tǒng)型船艏的航道周邊有明顯的冰尖形狀，而非傳統(tǒng)型的破冰航道邊緣更平滑。最后還可以發(fā)現(xiàn)，傳統(tǒng)尖瘦型船艏的破冰過程形成的碎冰塊普遍偏大，會對后續(xù)船舶的航道航行造成一定的影響，非傳統(tǒng)平緩船艏的破冰塊的形狀和大小分布較均勻，且尺寸偏小。

2.3.3 兩種船艏的破冰載荷

兩種破冰船艏破冰力和經(jīng)驗公式的對比如圖6 所示。從破冰載荷來看，兩種船艏的載荷趨勢和破冰模式相對應，傳統(tǒng)型艏的冰載荷具有連續(xù)性的特征，這是由于傳統(tǒng)型船艏破冰過程沒有明確的周期性分界線，往往下一個破冰循環(huán)過程和本次破冰循環(huán)過程的裂紋擴展同時進行；非傳統(tǒng)型艏的冰載荷具有明顯的周期循環(huán)特征，每個周期的冰載荷呈現(xiàn)先增加后減小的趨勢，增加過程對應于弧狀裂紋的產(chǎn)生和擴展，減小的過程對應于放射狀裂紋的擴展進而該區(qū)域冰全部破碎為小冰塊。

圖6 兩種船艏型的破冰載荷對照Fig.6 Comparison of icebreaking loads of two bow types

3 結 論

本文建立了基于粒子法碰撞問題的接觸檢測算法，編譯了基于MPI并行計算的數(shù)值模擬程序，預報了不同破冰艏的破冰過程和破冰載荷，并進行了對比分析，得到了以下結論：

（1）本文建立的數(shù)值模型能夠準確地模擬冰與船體的相互作用過程；

（2）傳統(tǒng)尖瘦型和非傳統(tǒng)平緩型破冰船艏具有不同的破冰模式，前者主要伴隨著環(huán)狀裂紋的發(fā)生擴展和少量放射狀的裂紋，后者主要呈現(xiàn)弧狀裂紋和放射狀裂紋的組合形式，裂紋擴展模式較簡單；

（3）傳統(tǒng)尖瘦型和非傳統(tǒng)平緩型破冰船艏具有不同趨勢的破冰載荷，前者呈現(xiàn)連續(xù)性的特征，后者呈現(xiàn)明顯周期性變化的特征。

冰塊碎裂后的浸沒、旋轉和沿著船體的滑動運動過程是分析破冰船性能另一關鍵因素。在后續(xù)的研究中，將增加碎冰塊的運動計算模型，建立破冰船的冰塊浸沒阻力預報方法。