潛艇破冰上浮近場(chǎng)動(dòng)力學(xué)模型

葉禮裕，王超，郭春雨，常欣

哈爾濱工程大學(xué)船舶工程學(xué)院，黑龍江哈爾濱150001

0 引 言

隨著全球氣候變暖，北極地區(qū)潛在的政治、經(jīng)濟(jì)、科技和軍事價(jià)值日益突顯，資源開采和航道開通的需求也在不斷提高。目前，極地地區(qū)已成為國(guó)際政治經(jīng)濟(jì)競(jìng)爭(zhēng)與合作的熱點(diǎn)地區(qū)，吸引著許多國(guó)家的關(guān)注。隨著人們對(duì)北極地理和氣候環(huán)境研究的深入，越來越深刻認(rèn)識(shí)到潛艇在北極地區(qū)所能發(fā)揮的政治和軍事價(jià)值。北極海冰的平均厚度可達(dá)到3 m，其特殊的環(huán)境為潛艇的活動(dòng)提供了特殊的隱蔽條件。基于北極地區(qū)特殊的地理位置，掌握潛艇在冰層下自由航行和活動(dòng)的技能，將形成特殊的地緣戰(zhàn)略控制能力，從而在大國(guó)之間的戰(zhàn)略競(jìng)爭(zhēng)中發(fā)揮特殊作用。要想勝任在北極地區(qū)的作戰(zhàn)任務(wù)，潛艇必須具有破冰上浮的能力。由文獻(xiàn)［1］可知，為了滿足北極的發(fā)展戰(zhàn)略，美、俄兩國(guó)現(xiàn)階段設(shè)計(jì)和建造的大部分核潛艇都需要滿足破冰要求。但是，潛艇在破冰上浮的過程中有可能會(huì)陷入危險(xiǎn)狀態(tài)，從而使?jié)撏У陌踩艿酵{。因此，開展?jié)撏票细∵^程的科學(xué)研究至關(guān)重要。

由于潛艇冰區(qū)航行的研究涉及到一個(gè)國(guó)家的關(guān)鍵技術(shù)，公開資料非常少，可查閱的資料均為國(guó)家層面的軍事競(jìng)爭(zhēng)。自二戰(zhàn)結(jié)束后，美、蘇兩國(guó)潛艇在冰面下的爭(zhēng)奪愈演愈烈。1946～1952年，美國(guó)派出了多艘常規(guī)動(dòng)力潛艇前往北極偵察和收集資料，均未能完成使命。常規(guī)動(dòng)力潛艇因續(xù)航力有限，需要頻繁浮出水面充電，故不適合在冰下長(zhǎng)期作戰(zhàn)，而核潛艇則具有無限續(xù)航能力，可長(zhǎng)時(shí)間潛航，因而成為角逐北冰洋的新戰(zhàn)略力量［2］。自核潛艇誕生以來，軍事家和戰(zhàn)略家們就開始探討將其部署到北冰洋以發(fā)動(dòng)戰(zhàn)略攻擊。1958年，美國(guó)第1艘核動(dòng)力潛艇“鸚鵡螺”號(hào)成功到達(dá)北極點(diǎn)。自此以后，美國(guó)海軍便開始強(qiáng)化其在北極的水下力量。2009年，美國(guó)海軍在北極地區(qū)進(jìn)行了特殊的軍事演習(xí)，訓(xùn)練潛艇在北極條件下的行動(dòng)能力。截至目前，美國(guó)海軍的多艘核潛艇已在北極軍事演練中從水下破冰而出。同樣，俄羅斯海軍在北極冰下也具有較強(qiáng)的實(shí)力，北極一直是其潛艇活動(dòng)的戰(zhàn)略性海域。1995年8月，俄羅斯海軍利用“臺(tái)風(fēng)”級(jí)核潛艇在北極冰下進(jìn)行了SS-N-20導(dǎo)彈發(fā)射試驗(yàn)。2006年，“臺(tái)風(fēng)”級(jí)核潛艇浮出冰面，沿著冰層裂縫破冰航行，并進(jìn)行了導(dǎo)彈試射。2009年，俄羅斯海軍“德爾塔-Ⅳ”型戰(zhàn)略核潛艇成功躲避美國(guó)的偵察，發(fā)射了2枚彈道導(dǎo)彈。

一直以來，人們都比較重視極地航行船舶的科學(xué)研究，已建立了可靠的實(shí)驗(yàn)和數(shù)值研究手段，但潛艇在極地地區(qū)運(yùn)行的安全性研究卻并未受到重視，到目前為止還未建立理論計(jì)算、數(shù)值模擬和實(shí)驗(yàn)測(cè)量等研究手段。對(duì)潛艇破冰上浮過程模擬的關(guān)鍵是對(duì)海冰破碎過程的模擬。目前，在海冰破碎過程模擬方法中，應(yīng)用最廣的是有限元法和離散元法。但由于海冰在不同應(yīng)變率下會(huì)呈現(xiàn)出2種不同的力學(xué)性質(zhì)，故采用有限元法難以準(zhǔn)確模擬冰的真實(shí)本構(gòu)以及處理冰的極端大變形問題。離散元法適用于非連續(xù)介質(zhì)力學(xué)問題的計(jì)算方法［3-6］，但在連續(xù)體階段的結(jié)果偏差制約了其在動(dòng)態(tài)破壞問題上的應(yīng)用。一些學(xué)者還將光滑粒子流體動(dòng)力學(xué)（Smoothed Particle Hydrodynamics，SPH）方法應(yīng)用于海冰動(dòng)力學(xué)問題的求解。盡管如此，采用數(shù)值方法預(yù)報(bào)和分析冰—艇體的相互作用問題依然是一個(gè)比較有效的方法，必須采用更為先進(jìn)的數(shù)值計(jì)算方法來解決此類問題。近場(chǎng)動(dòng)力學(xué)方法是一種無網(wǎng)格方法，在計(jì)算斷裂和大尺度變形問題方面具有獨(dú)特的優(yōu)勢(shì)［7-8］，非常適合海冰破碎問題［9］。近場(chǎng)動(dòng)力學(xué)方法與SPH方法均為粒子法，但兩者的求解思路不同：近場(chǎng)動(dòng)力學(xué)方法是借助物質(zhì)點(diǎn)法和分子動(dòng)力學(xué)思想，對(duì)連續(xù)介質(zhì)力學(xué)進(jìn)行分析，并以積分形式構(gòu)建物體的運(yùn)動(dòng)方程；而SPH方法則是通過在計(jì)算域內(nèi)布置粒子，將這些粒子的變量通過核函數(shù)進(jìn)行插值近似來實(shí)現(xiàn)積分計(jì)算。

本文將采用近場(chǎng)動(dòng)力學(xué)方法解決冰—艇體的接觸問題，并模擬海冰破碎過程。首先，將艇體表面看成是剛性固體且不發(fā)生變形，采用接觸檢測(cè)算法識(shí)別冰—艇體的接觸位置。然后，基于近場(chǎng)動(dòng)力學(xué)方法和接觸識(shí)別算法，開發(fā)潛艇破冰上浮求解程序，通過跟蹤每一個(gè)時(shí)刻的破冰狀態(tài)來計(jì)算海冰對(duì)潛艇殼體的動(dòng)態(tài)冰載荷力。最后，以美國(guó)DARPA潛艇模型SUBOFF為研究對(duì)象，開展?jié)撏票^程的數(shù)值預(yù)報(bào)和特性研究，并對(duì)計(jì)算結(jié)果進(jìn)行分析。

1 近場(chǎng)動(dòng)力學(xué)理論

1.1 基本方程

近場(chǎng)動(dòng)力學(xué)是一種非局部連續(xù)方法，其提供了一種更加實(shí)用的粒子方法，能從宏觀角度模擬大尺度破壞和變形問題。近場(chǎng)動(dòng)力學(xué)方法將連續(xù)介質(zhì)離散為均勻的物質(zhì)點(diǎn)。在參考坐標(biāo)系下，物質(zhì)點(diǎn)x的位置與時(shí)間t的物質(zhì)點(diǎn)運(yùn)動(dòng)方程為［8］

式中：Hx為由臨近物質(zhì)點(diǎn)x的所有其他物質(zhì)點(diǎn)所構(gòu)成的域；Vx′為物質(zhì)點(diǎn)x′的體積；u為物質(zhì)點(diǎn)x的位移；ρ為材料密度；f為物質(zhì)點(diǎn)x′與物質(zhì)點(diǎn)x之間相互作用的力密度；b為物質(zhì)點(diǎn)受到的外力。由式（1）可知，近場(chǎng)動(dòng)力學(xué)運(yùn)動(dòng)方程采用的是空間積分形式，故可在任何位置處進(jìn)行計(jì)算。

為表述方便，將2個(gè)物質(zhì)點(diǎn)之間的相對(duì)位置定義為ξ=x′-x，2個(gè)物質(zhì)點(diǎn)之間的相對(duì)位移定義為η=u(x′，t)-u(x，t)。于是，2個(gè)物質(zhì)點(diǎn)在 t時(shí)刻的相對(duì)位置可以表示為ξ+η。力密度可以表示為f(η，ξ)，其大小依賴于ξ和η。2個(gè)物質(zhì)點(diǎn)之間的相互作用可以稱為“鍵”，類似于一對(duì)相互作用的彈簧力，如圖1所示。近場(chǎng)動(dòng)力學(xué)定義了一個(gè)近場(chǎng)范圍尺寸δ，物質(zhì)點(diǎn)x只與距離δ半徑的球形范圍內(nèi)的物質(zhì)點(diǎn)x′相互作用，即

對(duì)于一個(gè)連續(xù)體，成對(duì)的力可以視為對(duì)物質(zhì)點(diǎn)之間的相互作用，并且滿足線動(dòng)量和角動(dòng)量的關(guān)系：

1.2 微觀彈脆性（PMB）模型

對(duì)于由微彈性材料構(gòu)成的物體，力密度函數(shù)可以表示為

式（1）和式（5）共同構(gòu)成了鍵型近場(chǎng)動(dòng)力學(xué)材料的本構(gòu)模型。為了描述物質(zhì)點(diǎn)之間鍵的伸長(zhǎng)程度，近場(chǎng)動(dòng)力學(xué)方法使用了非線性彈性材料的鍵伸長(zhǎng)率：

當(dāng)鍵處于拉伸狀態(tài)時(shí)，s為正值。由于s不依賴于ξ的方向，因此這種材料是各向同性的。為了在本構(gòu)模型中引入材料破壞的概念，使用了一個(gè)比較簡(jiǎn)單的判定條件，即當(dāng)某一時(shí)刻的鍵伸長(zhǎng)率超過一定值時(shí)，認(rèn)為物質(zhì)點(diǎn)之間的鍵永久發(fā)生斷裂。從此刻起，此對(duì)物質(zhì)點(diǎn)之間的力密度為0。該條件可稱為近場(chǎng)動(dòng)力學(xué)的破壞準(zhǔn)則。因此，對(duì)于典型的微觀彈脆性（Prototype Micro-elastic Brittle，PMB）材料，力密度函數(shù)可以簡(jiǎn)化為

式中，g為線性標(biāo)量函數(shù)，可定義為

式中，c為均勻各向同性材料的標(biāo)量微模量常數(shù)，可定義為

式中，κ為體積模量。

式（7）中的μ為歷史變形判斷標(biāo)準(zhǔn)，當(dāng)鍵伸長(zhǎng)率達(dá)到破壞準(zhǔn)則時(shí)，鍵發(fā)生斷裂，μ=0，反之，μ=1，可以定義為

式中，s0為鍵破壞時(shí)鍵伸長(zhǎng)率的極限值，可稱為極限伸長(zhǎng)率。對(duì)于典型的PMB材料，s0可由下式計(jì)算得到：

式中，G0為裂縫擴(kuò)展時(shí)的能量釋放率。盡管彈脆性材料在初始條件下是各向同性的，但是有些鍵斷裂之后會(huì)引起材料的各向異性。為此，引入鍵的破壞水平參數(shù)

式中，Vξ為物質(zhì)點(diǎn)ξ的體積。

1.3 方程離散化

為了便于開展數(shù)值計(jì)算，需要將連續(xù)體材料離散成物質(zhì)點(diǎn)，然后對(duì)每個(gè)物質(zhì)點(diǎn)進(jìn)行積分運(yùn)算［9］。離散后，每個(gè)物質(zhì)點(diǎn)x的運(yùn)動(dòng)方程可以表示為

式中：m為時(shí)間步；下標(biāo)i和j為物質(zhì)點(diǎn)的編號(hào)，其中i為要計(jì)算的物質(zhì)點(diǎn)，j為臨近x的物質(zhì)點(diǎn)；Vj為物質(zhì)點(diǎn)j的體積。

通過式（13），可以求得每個(gè)物質(zhì)點(diǎn)的加速度。為了得到每個(gè)物質(zhì)點(diǎn)的位移，可以通過中央差分法計(jì)算得到：

式中，Δt為時(shí)間步長(zhǎng)。

2 潛艇破冰上浮計(jì)算方法

2.1 問題描述

北極大部分地區(qū)被海冰覆蓋，極大地限制了水面艦船的航行，但潛艇卻可以潛入冰層之下自由航行，不僅能躲避反潛巡邏機(jī)和反潛水面艦船的偵察與攻擊，而且極地地區(qū)冰層的膨脹、破裂和因摩擦產(chǎn)生的環(huán)境噪聲還能掩蓋潛艇本身的噪聲。考慮到北極地區(qū)作戰(zhàn)任務(wù)的需要，要求潛艇不僅要能在冰下長(zhǎng)久航行，關(guān)鍵時(shí)刻還需要能夠破冰上浮，發(fā)射彈道導(dǎo)彈或者對(duì)外通信聯(lián)絡(luò)或者接受救援補(bǔ)給。

圖2所示為潛艇破冰上浮情況。在潛艇破冰上浮的過程中，主要受力部位是指揮室圍殼和艇的上層建筑，通常需要對(duì)這些部位進(jìn)行結(jié)構(gòu)優(yōu)化設(shè)計(jì)，以滿足強(qiáng)度要求。極地地區(qū)的冰面通常無限寬廣，潛艇在上浮過程中只能使殼體周圍的海冰破碎。

參考潛艇破冰上浮的實(shí)際情況，建立如圖3所示的數(shù)值模型。該模型假定潛艇以一定的速度上浮，并與上方靜止的冰面接觸。為便于數(shù)值計(jì)算，將平整冰簡(jiǎn)化為長(zhǎng)方體，而且長(zhǎng)度和寬度遠(yuǎn)大于厚度。由于實(shí)際的海冰尺寸相對(duì)于艇體要大得多，因此在進(jìn)行模擬時(shí)冰塊的長(zhǎng)度和寬度要足夠大，兩邊加固定邊界條件以代替無限大域的邊界問題。另外，考慮到潛艇在破冰上浮過程中流體對(duì)冰載荷的影響較小，忽略了流體的干擾作用。本文的主要目的是計(jì)算冰對(duì)艇體的接觸力并模擬海冰的破碎過程，故將潛艇殼體看成是剛性固體，暫不考慮殼體變形的影響。

為了考慮海冰粒子重力和浮力的影響，根據(jù)浮力和重力計(jì)算公式，當(dāng)物質(zhì)點(diǎn)位于冰排上方時(shí),式（1）中外力的取值為b(x，t)=ρVg；當(dāng)物質(zhì)點(diǎn)位于冰排下時(shí)，外力的取值為b(x，t)=0.1ρVg，其中V為物質(zhì)點(diǎn)的體積。

2.2 模型數(shù)值離散化

本文的海冰模型是采用近場(chǎng)動(dòng)力學(xué)方法進(jìn)行模擬，在數(shù)值計(jì)算過程中，需要將冰塊連續(xù)體離散為物質(zhì)點(diǎn)形式。但是，需要特別注意物質(zhì)點(diǎn)間距Δx和鄰域半徑δ的選擇，它們是影響計(jì)算精度的重要參數(shù)，需要綜合考慮計(jì)算精度和速度，以選取最優(yōu)值。當(dāng)鄰域半徑過小時(shí)，裂縫擴(kuò)展過程對(duì)網(wǎng)格有依賴性，難以模擬材料的大尺度變形，因此建議采用δ=3Δx作為鄰域半徑［10］。

對(duì)于艇體模型，本文將艇體簡(jiǎn)化為剛性表面，并用四邊形面元來劃分。本文以美國(guó)的SUBOFF潛艇模型為研究對(duì)象，該模型的主要參數(shù)見第3.1節(jié)。將艇體三維模型離散為一系列四邊形面元形式，面元越小，越能逼近艇體的真實(shí)形狀。圖4所示為艇體表面網(wǎng)格劃分結(jié)果。

2.3 海冰物質(zhì)點(diǎn)與潛艇殼體的接觸識(shí)別

當(dāng)海冰與潛艇殼體表面離散化后，可將海冰與殼體表面的接觸判斷轉(zhuǎn)換為物質(zhì)點(diǎn)和四邊形面元的相對(duì)位置判斷，從而將復(fù)雜的接觸過程簡(jiǎn)化為一系列點(diǎn)和面的數(shù)學(xué)問題，有效降低接觸檢測(cè)計(jì)算的復(fù)雜度。

為減少計(jì)算量和運(yùn)算的復(fù)雜度，同時(shí)保證接觸檢測(cè)的精度和效率，基于接觸檢測(cè)中包圍盒的思想［11］，本文給出了一種簡(jiǎn)單高效的方法來剔除不可能與艇體表面發(fā)生接觸的物質(zhì)點(diǎn)。建立一個(gè)長(zhǎng)方體將整個(gè)艇體包圍在內(nèi)，該長(zhǎng)方體的長(zhǎng)度等于艇長(zhǎng)Ls，寬度等于艇寬Bs，高度等于艇體最高點(diǎn)到最低點(diǎn)z方向（垂向）的距離。海冰的物質(zhì)點(diǎn)位于長(zhǎng)方體內(nèi)部被認(rèn)為可能與艇體接觸的物質(zhì)點(diǎn)。

對(duì)于每一個(gè)可能與艇體接觸的物質(zhì)點(diǎn)，若與艇體發(fā)生接觸，只能有一個(gè)四邊形面元與物質(zhì)點(diǎn)接觸。為了提高搜索效率，通過物質(zhì)點(diǎn)和四邊形4個(gè)點(diǎn)的坐標(biāo)位置關(guān)系來確定可能與該物質(zhì)點(diǎn)接觸的面元。假設(shè)粒子的坐標(biāo)為(x0，y0，z0)，對(duì)于艇體表面的所有四邊形面元，可以找到其4個(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ā)生碰撞。

對(duì)于每一個(gè)可能與艇體接觸的物質(zhì)點(diǎn)，可計(jì)算得到可能與其發(fā)生碰撞面元的控制點(diǎn)之間的距離，從而找到距離最短的面元。通過點(diǎn)與面的數(shù)學(xué)公式，判斷物質(zhì)點(diǎn)與其距離最短面元的位置關(guān)系，如式（15）所示。

式中，A1，B1，C1，D1為四邊形面元所在空間平面方程的系數(shù)。

于是，與艇體表面發(fā)生接觸的物質(zhì)點(diǎn)和與其發(fā)生接觸的面元即可確定。

2.4 接觸作用力計(jì)算

潛艇在上浮的過程中，海冰物質(zhì)點(diǎn)將與潛艇殼體發(fā)生接觸。接觸以后，物質(zhì)點(diǎn)就會(huì)滲入到艇體內(nèi)，如圖5（b）所示。為了反映真實(shí)的物理情況，滲入到艇體內(nèi)的物質(zhì)點(diǎn)需要進(jìn)行位置再分配，將該物質(zhì)點(diǎn)分配到臨近處的艇體表面，如圖5（c）所示。通過下式，可以計(jì)算出滲入到艇體內(nèi)的物質(zhì)點(diǎn)和新分配物質(zhì)點(diǎn)的距離為

新分配物質(zhì)點(diǎn)的坐標(biāo)可以通過下式計(jì)算：

式中：下標(biāo)k為物質(zhì)點(diǎn)標(biāo)號(hào)；v0為物質(zhì)點(diǎn)k的速度；n為接觸點(diǎn)處的法向量。

物質(zhì)點(diǎn)x(k)在t+Δt時(shí)刻的速度可用式（18）計(jì)算：

3 潛艇破冰上浮算例分析

3.1 計(jì)算模型及參數(shù)設(shè)置

海冰破碎方式的多樣性以及不同海域海冰的物理和力學(xué)性質(zhì)有很大的不同，導(dǎo)致海冰的物理和力學(xué)性質(zhì)不太容易確定，較好的方式是通過采樣相關(guān)海域的海冰來試驗(yàn)和測(cè)量其物理性質(zhì)，例如壓縮強(qiáng)度、彈性模量和密度等。由于本文關(guān)注的重點(diǎn)不在于海冰的物理性質(zhì)，而是對(duì)海冰破壞現(xiàn)象的模擬，因此將海冰簡(jiǎn)化為勻質(zhì)材料，而海冰的物理和力學(xué)參數(shù)與文獻(xiàn)［12］相同，即彈性模量E=1.8 GPa，泊松比υ=0.25，臨界伸長(zhǎng)率s0=0.02，密度ρ=900 kg/m3。

本文的研究對(duì)象選用的是美國(guó)DARPA SUBOFF實(shí)尺度潛艇模型，如圖6所示。模型的型值取自文獻(xiàn)［13］，全附體潛艇長(zhǎng)104.544 m，艇身最大直徑12.192 m，指揮室圍殼長(zhǎng)8.832 m，平行中體為53.496 m，尾翼的翼型后緣位于96.168 m處，4個(gè)尾翼剖面均為NACA 0020翼型，并呈對(duì)稱布置［14-15］。

3.2 計(jì)算模型的收斂性分析

在開展?jié)撏票细√匦匝芯恐埃紫柔槍?duì)物質(zhì)點(diǎn)間隔Δx和時(shí)間步長(zhǎng)Δt對(duì)計(jì)算結(jié)果的影響進(jìn)行分析。將冰層尺度定為長(zhǎng)度L=150 m、寬度B=36 m、厚度T=0.9 m的長(zhǎng)方體，分析潛艇以0.5 m/s速度破冰上浮過程中冰載荷的動(dòng)態(tài)響應(yīng)特性。由于潛艇在破冰上浮過程中z軸方向（垂向）的力遠(yuǎn)大于x軸方向（縱向）和y軸方向（橫向）的力，故本文重點(diǎn)分析z軸方向的冰載荷。

為了分析物質(zhì)點(diǎn)間距對(duì)計(jì)算結(jié)果的影響，在其他工況相同的情況下，將冰層離散成間距分別為Δx=L/250，L/500和L/750的物質(zhì)點(diǎn)，時(shí)間步長(zhǎng)控制為Δt=0.006 25 s。圖7給出了冰載荷Fz隨時(shí)間變化的曲線。隨著物質(zhì)點(diǎn)間距的減小，計(jì)算曲線逐漸收斂。物質(zhì)點(diǎn)間距Δx=L/250的冰載荷大小要比間距Δx=L/500和L/750的稍小，曲線振蕩也比較劇烈。由于冰和結(jié)構(gòu)物相互作用有很大的隨機(jī)性，冰載荷將發(fā)生劇烈的動(dòng)態(tài)變化，因此相同工況下的冰載荷時(shí)程曲線無法完全重合。由圖可知，物質(zhì)點(diǎn)間距Δx=L/500和L/750的計(jì)算曲線基本一致。因此可以認(rèn)為，當(dāng)物質(zhì)點(diǎn)間距小于Δx=L/500時(shí)，計(jì)算結(jié)果已能夠達(dá)到收斂。

同理，控制物質(zhì)點(diǎn)間距為Δx=L/500不變，時(shí)間步長(zhǎng)分別取為 Δt=0.000 948，0.000 625和0.000 312 s。圖8所示為不同時(shí)間步長(zhǎng)時(shí)冰載荷Fz隨時(shí)間變化的曲線。由圖8可知，盡管3條曲線在某些時(shí)間點(diǎn)上有一些振蕩，但總體來看基本還是一致的。因此可以認(rèn)為當(dāng)時(shí)間步長(zhǎng)小于0.000 948 s時(shí)，計(jì)算結(jié)果可達(dá)到收斂。

綜上所述，本文計(jì)算模型的收斂性可以得到保證。在后續(xù)的計(jì)算中，兼顧計(jì)算速度和精度，物質(zhì)點(diǎn)間距取為 Δx=L/500，時(shí)間步長(zhǎng)取為 Δt=0.000 625 s。

3.3 潛艇破冰上浮動(dòng)態(tài)特性分析與驗(yàn)證

本節(jié)以潛艇以0.5 m/s速度連續(xù)破冰上浮撞破0.9 m厚的冰層為例，來分析潛艇破冰上浮過程中海冰的破碎特點(diǎn)和冰載荷的動(dòng)態(tài)響應(yīng)特性。選擇的冰層是長(zhǎng)度L=150 m、寬度B=36 m的長(zhǎng)方體。將離散后的海冰物質(zhì)點(diǎn)之間的間隔取為Δx=L/500，將時(shí)間步長(zhǎng)設(shè)定為Δt=0.000 625 s。

圖9給出了不同時(shí)刻潛艇破冰上浮時(shí)海冰的擠壓和彎曲破壞過程。圖9中云圖的數(shù)值表示冰物質(zhì)點(diǎn)的破碎水平，顏色的深度表示冰粒子的損傷度，其中藍(lán)色表示冰粒子處于無損狀態(tài)，紅色代表冰粒子處于完全損傷狀態(tài)。由圖9可知，潛艇指揮室圍殼上部首先撞擊到冰面，由于它垂直于冰面，圍繞著指揮室圍殼的海冰破碎，指揮室圍殼首先破冰而出，并且其上部還出現(xiàn)了一塊獨(dú)立于冰面的海冰，伴隨著指揮室圍殼一起向上運(yùn)動(dòng)。隨后，潛艇尾翼與海冰接觸，沖破冰層露出冰面。隨著破冰過程的進(jìn)行，潛艇上層建筑也與海冰接觸，冰面出現(xiàn)了大面積的擠壓和彎曲破壞。

下面，將通過與網(wǎng)上可找到的美、俄核潛艇在北極冰下破冰而出的畫面進(jìn)行對(duì)比來驗(yàn)證本文計(jì)算模型的有效性。圖10所示為美國(guó)核潛艇在北極地區(qū)破冰而出的最終畫面，將其與圖9（f）進(jìn)行比較后發(fā)現(xiàn)，潛艇上層建筑的海冰分布基本一致。圖9（f）給出了海冰物質(zhì)點(diǎn)損傷度的分布，可以看出，潛艇上層建筑的海冰受到潛艇的沖擊作用發(fā)生了破碎。潛艇破冰上浮后，只有指揮室圍殼和尾翼露出冰面，而潛艇上層建筑則被掩蓋在海冰之下。

圖11所示為本文計(jì)算的潛艇指揮室圍殼和美國(guó)核潛艇指揮室圍殼破冰而出的局部放大圖。通過比較可知，兩者的海冰破壞形式基本一致，指揮室圍殼上部均出現(xiàn)了一塊獨(dú)立于冰面的海冰。由圖11（a）可知，在當(dāng)前計(jì)算工況下，指揮室圍殼沖擊海冰的過程中，其前端的海冰發(fā)生了縱向剪切破壞，而后端的海冰則以彎曲破壞為主。可以推測(cè)，在指揮室圍殼沖擊海冰的過程中，海冰的破壞模式與指揮室圍殼形狀有關(guān)，當(dāng)然也可能與潛艇的上浮速度和海冰的厚度等有關(guān)。

圖12所示為本文計(jì)算的潛艇尾翼和美國(guó)核潛艇尾翼破冰而出的局部放大圖。通過比較可知，在當(dāng)前計(jì)算工況下，兩者尾翼在沖擊海冰的過程中均以彎曲破壞為主，破壞形式基本一致，這種破壞模式不同于指揮室圍殼沖擊海冰的過程。這主要是由于與指揮室圍殼相比，尾翼更加窄小，降低了與海冰接觸的面積，因而受力面積更小，使得海冰更加易于破碎。

圖13所示為當(dāng)前計(jì)算工況下潛艇上浮過程中冰載荷Fz的時(shí)程曲線。在剛開始階段，由于指揮室圍殼對(duì)海冰的沖擊作用，潛艇受到的冰載荷出現(xiàn)小幅度的波動(dòng)。3.0 s以后，潛艇受到的冰載荷進(jìn)入一個(gè)比較平穩(wěn)的過程，大小約為168.3 kN。在該過程中，海冰受到指揮室圍殼的沖擊作用結(jié)束，指揮室圍殼破冰而出，此時(shí)潛艇主要受指揮室圍殼上部一塊獨(dú)立于冰面的海冰的作用，因此該過程中潛艇受到的冰載荷約等于這塊海冰的重力。通過測(cè)量，得到指揮室圍殼頂部面積S約為28.1 m2，而海冰的厚度T=0.9 m，從而可以大概估計(jì)出這塊海冰的重力為228.2 kN，稍大于該過程所受到的冰載荷，從而驗(yàn)證了本文接觸力計(jì)算方法的有效性。在8.5～9.0 s之間，冰載荷在波動(dòng)中上升，該時(shí)間段海冰受到尾翼和潛艇上殼體的共同作用，因此冰載荷要大于指揮室圍殼與海冰作用的峰值。在9.0～10.5 s之間，冰載荷迅速增大，該時(shí)刻曲線比較平滑，可以推測(cè)在該過程中潛艇與海冰主要發(fā)生了擠壓作用。在10.5 s以后，冰載荷開始出現(xiàn)劇烈的振蕩，這可能是由于在該過程中海冰不斷發(fā)生擠壓和彎曲破壞所引起，這將會(huì)導(dǎo)致艇體的激烈振動(dòng)，有可能會(huì)對(duì)艇體結(jié)構(gòu)安全性構(gòu)成威脅。當(dāng)潛艇上殼體上升到一定程度以后，冰載荷開始迅速減小，表明破冰過程將要結(jié)束。

為了更好地分析潛艇指揮室圍殼與海冰相互作用的冰載荷特性，單獨(dú)繪制了該過程冰載荷Fz的時(shí)程曲線，如圖14所示。指揮室圍殼上部與海冰發(fā)生接觸后，海冰在未發(fā)生破碎之前主要以擠壓破壞為主，冰載荷迅速增大，此時(shí)曲線比較平滑。隨后，海冰發(fā)生彎曲或者剪切破壞，冰載荷曲線發(fā)生劇烈振蕩。當(dāng)潛艇指揮室圍殼破冰而出時(shí)，此時(shí)海冰受到指揮室圍殼的沖擊作用結(jié)束，冰載荷迅速減小。

4 結(jié) 論

本文基于近場(chǎng)動(dòng)力學(xué)方法，建立了潛艇破冰上浮數(shù)值計(jì)算模型，考慮了海冰浮力和重力的影響，并將其加入到了近場(chǎng)動(dòng)力學(xué)基本方程中。通過運(yùn)用計(jì)算程序?qū)撏У钠票细∵^程和冰載荷動(dòng)態(tài)響應(yīng)進(jìn)行計(jì)算，分析得到了以下結(jié)論：

1）開展了潛艇破冰上浮過程中物質(zhì)點(diǎn)間隔和時(shí)間步長(zhǎng)的收斂性分析，計(jì)算結(jié)果表明，物質(zhì)點(diǎn)間隔小于Δx=L/500、時(shí)間步長(zhǎng)小于0.000 948 s時(shí)，計(jì)算結(jié)果可達(dá)到收斂。

2）基于本文建立的數(shù)值模型對(duì)潛艇破冰上浮過程中海冰的破壞過程進(jìn)行模擬，很好地再現(xiàn)了海冰的擠壓、縱向剪切破壞或彎曲破壞過程，證明了本文計(jì)算方法的有效性以及接觸識(shí)別算法的可靠性。

3）對(duì)潛艇破冰上浮過程冰載荷時(shí)程曲線的分析表明，在海冰擠壓破壞過程中，冰載荷迅速增大，曲線比較平滑，但海冰彎曲破壞過程一旦發(fā)生，冰載荷將出現(xiàn)劇烈的振蕩，對(duì)艇體結(jié)構(gòu)的安全性不利。

本文建立的潛艇破冰上浮過程計(jì)算模型，和對(duì)潛艇破冰上浮的動(dòng)態(tài)特性進(jìn)行的計(jì)算分析，可以用于指導(dǎo)潛艇殼體結(jié)構(gòu)的優(yōu)化設(shè)計(jì)。本文并未考慮流體與艇體和海冰之間相互作用的影響，計(jì)算結(jié)果可能存在誤差。下一步，將更全面地考慮各種因素的影響，建立更加可靠的數(shù)值預(yù)報(bào)方法。

參考文獻(xiàn)：

［1］劉征魯，李春迪.核潛艇破冰要有哪些“硬功夫”？［N］.中國(guó)國(guó)防報(bào)，2016-05-06（014）.

［2］陸俊元.北極地緣政治與中國(guó)應(yīng)對(duì)［M］.北京：時(shí)事出版社，2010.LU J Y.Geopolitics in the Arctic and China's response［M］.Beijing：Current Affairs Press，2010（in Chi?nese）.

［3］WILLIAMS J R，HOCKING G，MUSTOE G G W.The theoretical basis of the discrete element method［C］//Proceeding of the International Conference on Numeri?cal Methods in Engineering：Theory and Applications.Rotterdam，Netherlands，1985.

［4］LAU M，LAWRENCE K P，ROTHENBURG L.Dis?crete element analysis of ice loads on ships and struc?tures［J］.Ships and Offshore Structures，2011，6（3）：211-221.

［5］HAN K，PERIC D，CROOK A J L，et al.A combined fi?nite/discrete element simulation of shot peening pro?cesses-Part I：studies on 2D interaction laws［J］.Engi?neering Computations，2000，17（5）：593-620.

［6］HASHASH Y M A.Large-scale numerical simulations via parallel computing：an application using discrete element modeling of granular material［C］//NCSA Stra?tegic Applications Program Status Report 03/28/2006.Champaign，IL：Center for Super Computing Applica?tions，University of Illinois，2006.

［7］SILLING S A.Reformulation of elasticity theory for dis?continuities and long-range forces［J］.Journal of the Mechanics and Physics of Solids，2000，48（1）：175-209.

［8］SILLING S A，ASKARI E.A meshfree method based on the peridynamic model of solid mechanics［J］.Com ?puters&Structures，2005，83（17/18）：1526-1535.

［9］YE L Y，WANG C，CHANG X，et al.Propeller-ice con?tact modeling with peridynamics［J］.Ocean Engineer?ing，2017，139：54-64.

［10］MADENCI E，OTERKUS E.Peridynamic theory and its applications［M］.New York：Springer，2014.

［11］NEZAMI E G，HASHASH Y M A，ZHAO D W，et al.A fast contact detection algorithm for 3-D discrete el?ement method［J］. Computers and Geotechnics，2004，31（7）：575-587.

［12］ZHAO G L，XUE Y Z，LIU R W，et al.Numerical sim?ulation of ice load for icebreaker based on peridynam?ic［C］//Proceedings of the ASME 35th International Conference on Ocean，Offshore and Arctic Engineer?ing.Busan，South Korea，2006.

［13］GROVES N C，HUANG T T，CHANG M S.Geometric characteristics of DARPA SUBOFF models（DTRC ModelNos.5470 and 5471）：reportDTRC/SHD-1298-01［R］.Bethesda，USA：David Taylor Research Center，1989.

［14］操盛文，吳方良.尺度效應(yīng)對(duì)全附體潛艇阻力數(shù)值計(jì)算結(jié)果的影響［J］. 中國(guó)艦船研究，2009，4（1）：33-37，42.CAO S W，WU F L.Investigation of scaling effects on numerical computation of submarine resistance［J］.Chinese Journal of Ship Research，2009，4（1）：33-37，42（in Chinese）.

［15］吳方良，吳曉光，馬運(yùn)義，等.潛艇實(shí)艇阻力預(yù)報(bào)方法研究［J］. 中國(guó)艦船研究，2009，4（3）：28-32.WU F L，WU X G，MA Y Y，et al.Method of predict?ing the total submerged resistance of submarines［J］.Chinese Journal of Ship Research，2009，4（3）：28-32（in Chinese）.