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

（哈爾濱工程大學(xué)，哈爾濱 150001）

0 引 言

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

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

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

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

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

1.1 冰的本構(gòu)模型

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

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

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

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

1.2 接觸檢測(cè)算法（FCDA）

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

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

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

（3）在t+1時(shí)刻冰粒子穿透矩形體和目標(biāo)面，接下來(lái)只需關(guān)注進(jìn)入矩形體的所有粒子；

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

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

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

冰粒子與船體的接觸檢測(cè)過(guò)程完成后，對(duì)進(jìn)入船體的粒子重新分配位置，并更新粒子的運(yùn)動(dòng)信息和計(jì)算冰載荷，該過(guò)程可直接參考文獻(xiàn)[15]第十章的處理過(guò)程，本文不再解釋。

1.3 信息傳遞接口并行技術(shù)（MPI）

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

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

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

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

2 破冰過(guò)程數(shù)值預(yù)報(bào)分析

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

2.1 兩種船艏形狀描述

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

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

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

2.2 計(jì)算輸入

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

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

2.3 計(jì)算結(jié)果

2.3.1 數(shù)值方法驗(yàn)證

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

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

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

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

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

2.3.2 兩種船艏的破冰模式

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

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

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

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

2.3.3 兩種船艏的破冰載荷

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

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

3 結(jié) 論

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

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

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

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

冰塊碎裂后的浸沒(méi)、旋轉(zhuǎn)和沿著船體的滑動(dòng)運(yùn)動(dòng)過(guò)程是分析破冰船性能另一關(guān)鍵因素。在后續(xù)的研究中，將增加碎冰塊的運(yùn)動(dòng)計(jì)算模型，建立破冰船的冰塊浸沒(méi)阻力預(yù)報(bào)方法。