陳 娟 許永秋 姚熊亮 位 莎

1哈爾濱工程大學(xué) 船舶工程學(xué)院，黑龍江 哈爾濱150001 2浙江省臺(tái)州市海事局，浙江 臺(tái)州318000

水下爆炸沖擊載荷的一種新型算法研究

陳 娟1許永秋2姚熊亮1位 莎1

1哈爾濱工程大學(xué) 船舶工程學(xué)院，黑龍江 哈爾濱150001 2浙江省臺(tái)州市海事局，浙江 臺(tái)州318000

近年來(lái)，越來(lái)越多的研究者都重視使用計(jì)算機(jī)的數(shù)值模擬對(duì)水下爆炸進(jìn)行研究。水下爆炸包含了一系列復(fù)雜的物理過(guò)程，如大變形、高度非均勻性、可變形邊界和自由表面。對(duì)于傳統(tǒng)的基于網(wǎng)格的數(shù)值方法，水下爆炸問(wèn)題的模擬是一項(xiàng)非常具有挑戰(zhàn)性的工作。介紹了無(wú)網(wǎng)格光滑粒子流體動(dòng)力學(xué)（SPH）方法在水下爆炸模擬中的應(yīng)用，并且提出了處理固壁和物質(zhì)交界面的算法，給出了水下爆炸過(guò)程具有代表性瞬間的壓力分布、速度分布以及氣泡的演變和脈動(dòng)，為艦船水下爆炸沖擊載荷的計(jì)算提供了一種新型有效的研究方法。

水下爆炸；沖擊波；氣泡；壁面反射

1 引言

艦船水下爆炸近場(chǎng)問(wèn)題涉及水中沖擊波的入射、透射、反射以及氣泡的脈動(dòng)等復(fù)雜物理現(xiàn)象，使得氣、液、固態(tài)結(jié)構(gòu)之間的瞬態(tài)相互作用強(qiáng)烈而復(fù)雜，因此，艦船近場(chǎng)水下爆炸結(jié)構(gòu)動(dòng)態(tài)響應(yīng)問(wèn)題還有許多問(wèn)題有待研究和解決。然而對(duì)于傳統(tǒng)的基于網(wǎng)格的數(shù)值計(jì)算方法，水下爆炸問(wèn)題的模擬是一項(xiàng)非常具有挑戰(zhàn)性的工作。一方面，由于存在大變形、運(yùn)動(dòng)物質(zhì)交界面、可變形邊界和自由表面等特性，使得基于網(wǎng)格的數(shù)值方法難以處理；另一方面，在整個(gè)水下爆炸過(guò)程中，高能炸藥的起爆過(guò)程的復(fù)雜性和大尺度等問(wèn)題使基于網(wǎng)格的方法面臨更多的困難。

近年來(lái)，由于無(wú)網(wǎng)格技術(shù)的發(fā)展和應(yīng)用日趨成熟，人們開(kāi)始將無(wú)網(wǎng)格方法應(yīng)用到可壓縮多介質(zhì)流的數(shù)值模擬中，如氣泡上升［1，2］、深水強(qiáng)爆炸［3］、水介質(zhì)緩沖［4］等問(wèn)題。光滑粒子流體動(dòng)力學(xué)（SPH）方法［5，6］具有無(wú)網(wǎng)格性質(zhì)和拉格朗日粒子特性，它應(yīng)用離散化的粒子來(lái)表示物質(zhì)，能夠很自然地對(duì)多介質(zhì)流進(jìn)行模擬，而且由于不受網(wǎng)格劃分的限制，可以解決在基于網(wǎng)格的數(shù)值方法中由于高壓、高能、大變形等導(dǎo)致網(wǎng)格畸變而計(jì)算崩潰的問(wèn)題。因此本文用此方法對(duì)水下爆炸過(guò)程的機(jī)理進(jìn)行研究，并且提出了處理固壁和運(yùn)動(dòng)物質(zhì)交界面的點(diǎn)對(duì)點(diǎn)算法，旨在為艦船水下爆炸沖擊載荷的計(jì)算提供一種新型有效的研究方法。

2 水下爆炸理論基礎(chǔ)

水下爆炸是指在極短時(shí)間內(nèi)，在水下的極小體積內(nèi)或面積上發(fā)生極大能量轉(zhuǎn)換的過(guò)程。水下爆炸大體可分為3個(gè)階段：裝藥的爆轟、沖擊波的產(chǎn)生和傳播、氣泡的形成和脈動(dòng)。當(dāng)炸藥在水中爆炸時(shí)，其周?chē)橘|(zhì)直接受到具有高溫、高速、高壓的爆炸產(chǎn)物作用。在裝藥和介質(zhì)的界面處，爆炸產(chǎn)物以極高的速度向周?chē)鷶U(kuò)散，強(qiáng)烈地壓縮著相鄰的水，使其壓力、密度、溫度突躍式地升高，形成初始沖擊波。隨著沖擊波在水中的傳播，爆炸產(chǎn)物在水中以氣泡的形式存在并不斷膨脹與壓縮，形成氣泡脈動(dòng)現(xiàn)象。

2.1 狀態(tài)方程

狀態(tài)方程用來(lái)定義固體或流體在各種不同狀態(tài)下的壓力和密度以及比內(nèi)能之間的函數(shù)關(guān)系，正確選取狀態(tài)方程中的參數(shù)對(duì)于計(jì)算結(jié)果至關(guān)重要。狀態(tài)方程一般采用半經(jīng)驗(yàn)半理論的公式，方程中的主要參數(shù)由試驗(yàn)確定。

1）TNT爆炸產(chǎn)物狀態(tài)方程

本文中TNT爆炸產(chǎn)物采用標(biāo)準(zhǔn)JWL狀態(tài)方程［7］，JWL狀態(tài)方程是典型的動(dòng)力學(xué)狀態(tài)方程，它是一種不顯含化學(xué)反應(yīng)、由實(shí)驗(yàn)方法確定參數(shù)的經(jīng)驗(yàn)狀態(tài)方程，能比較精確地描述爆轟產(chǎn)物的膨脹驅(qū)動(dòng)做功，其形式為：

式中：p為爆轟產(chǎn)物的壓力；η為爆炸氣體密度與初始炸藥密度的比值，即η＝；A，B，R1，R2和ω為與炸藥狀態(tài)有關(guān)的常數(shù)；e為T(mén)NT單位質(zhì)量的內(nèi)能。

2）水的狀態(tài)方程

在高壓、高密度和高溫的沖擊載荷下，水的狀態(tài)方程也多種多樣，其中最常用的是Mie－Gruneisen狀態(tài)方程，其具體形式取決于水的狀態(tài)。在壓縮狀態(tài)下水的壓力為：

在膨脹狀態(tài)下水的壓力為：

式中：ρ0為水的初始密度；η＝為水?dāng)_動(dòng)前后的密度比；μ＝η－1，當(dāng)μ＞0時(shí)，水處于壓縮狀態(tài)，當(dāng)μ＜0時(shí)，水處于膨脹狀態(tài)；C為聲速；γ0為Gruneisen系數(shù)；a為體積修正系數(shù)；S1，S2，S3為擬合系數(shù)。

2.2 水下爆炸沖擊波的經(jīng)驗(yàn)計(jì)算公式

為了驗(yàn)證數(shù)值模擬結(jié)果的正確性，采用經(jīng)驗(yàn)公式進(jìn)行比較。水下爆炸沖擊波的經(jīng)驗(yàn)計(jì)算公式都是建立在爆炸相似律分析和對(duì)水下爆炸試驗(yàn)數(shù)據(jù)進(jìn)行擬合的基礎(chǔ)上的，由于試驗(yàn)數(shù)據(jù)的不同、擬合方法的不同，經(jīng)驗(yàn)公式的形式和使用范圍也不太一致，其中庫(kù)爾的經(jīng)驗(yàn)公式［8］是較為經(jīng)典的公式，其具體形式如下。沖擊波峰值壓力經(jīng)驗(yàn)公式為：

式中：pm為沖擊波峰值壓力，MPa；W為裝藥量，kg；R為爆距，m；k為實(shí)驗(yàn)系數(shù)，取52.4；α為壓力衰減系數(shù)，對(duì)標(biāo)準(zhǔn)TNT炸藥，取1.13。

3 SPH數(shù)值計(jì)算方法

3.1 SPH形式的水下爆炸控制方程

假設(shè)爆炸氣體和周?chē)橘|(zhì)是無(wú)粘性的，且整個(gè)水下爆炸過(guò)程被認(rèn)為是絕熱的，因此采用Euler方程作為控制方程。應(yīng)用SPH核近似和粒子近似可得到以下一系列SPH方程，且在動(dòng)量方程和能量方程中引入了Monaghan型人工粘性∏ij［9］，它不僅將動(dòng)能轉(zhuǎn)化為熱能，提供了沖擊波波面必不可少的耗散，而且防止計(jì)算過(guò)程中粒子之間的非物理穿透。

Monaghan型人工粘性∏ij表達(dá)式為：

光滑長(zhǎng)度h在SPH近似方法中非常重要，直接影響到求解的精度和效率。在早期的SPH應(yīng)用中，光滑長(zhǎng)度取決于系統(tǒng)的初始平均密度，而且在整個(gè)計(jì)算模擬過(guò)程中光滑長(zhǎng)度保持不變。但是對(duì)于水下爆炸問(wèn)題，由于其為高能高壓多介質(zhì)流動(dòng)問(wèn)題，應(yīng)用常值光滑長(zhǎng)度不僅難以對(duì)整個(gè)物理過(guò)程模擬再現(xiàn)，而且爆炸相關(guān)局部細(xì)節(jié)的信息也容易在計(jì)算中被遺漏。同時(shí)，在水下爆炸的模擬中使用常值光滑長(zhǎng)度不僅數(shù)值精度難以保證，而且穩(wěn)定性差，易于造成計(jì)算崩潰。因此，必須對(duì)光滑長(zhǎng)度進(jìn)行自適應(yīng)動(dòng)態(tài)變換，本文采用Benz［10］提出的動(dòng)態(tài)自適應(yīng)方法：

對(duì)上式進(jìn)行SPH近似，可得到以下形式：

下一個(gè)時(shí)間步長(zhǎng)的光滑長(zhǎng)度變?yōu)椋?/p>

3.2 固壁邊界處理

本文采用固壁粒子配合鏡像粒子的方法來(lái)模擬固定邊界。固壁粒子施加排斥力的方法防止粒子穿越固定邊界，排斥力的大小由式（11）確定，固壁粒子參與流體粒子密度及內(nèi)力的計(jì)算，本身的密度不斷更新，但位置和速度保持不變。

式中：r0為截止半徑，在本文算例中取為初始粒子間距，即1.0／100 m；rij＝ri－rj為粒子i和粒子j的距離；參數(shù)n1，n2分別取為6和4；參數(shù)D一般取與速度最大值的平方相等的量級(jí)，本文取為1.0e6。

圖1中粗直線表示剛性壁，圓表示粒子i的支持域。設(shè)j位于i的支持域內(nèi)，若實(shí)粒子j關(guān)于剛性壁對(duì)稱(chēng)的鏡像粒子j′也在i的支持域中，則粒子i鄰域內(nèi)的點(diǎn)也包括鏡像粒子j′，如果與k對(duì)應(yīng)的鏡像粒子k′落在粒子i的支持域外，則計(jì)算時(shí)不予考慮。由于粒子i也在其支持域內(nèi)，因此其鏡像粒子i′同樣需要考慮。

圖1 固壁粒子和鏡像粒子示意圖

3.3 交界面處理

在SPH中，物質(zhì)交界面處理是一個(gè)非常關(guān)鍵的問(wèn)題，因?yàn)镾PH方法具有拉格朗日性質(zhì)和粒子性質(zhì)，在整個(gè)演變過(guò)程中，來(lái)自不同介質(zhì)的相互接觸的粒子可能會(huì)隨著運(yùn)動(dòng)而分離，甚至有可能不再成為相鄰粒子。在水下爆炸模擬中的物質(zhì)交界面處理更為復(fù)雜，因?yàn)楦邏焊吣艿谋怏w和水之間的激烈相互作用會(huì)使粒子的運(yùn)動(dòng)相當(dāng)自由，本文采用的方法是，離得近但材料不同的兩粒子可以認(rèn)為是相鄰粒子，這樣能夠減少在邊界附近的粒子缺陷問(wèn)題，但為了防止粒子的非物理穿透或摻雜問(wèn)題，在交界面附近的不同粒子之間，當(dāng)其趨向于穿透時(shí)，即：pe＝≥1，施加一個(gè)排斥力PBij，其表達(dá)式為：

4 水下爆炸過(guò)程的數(shù)值模擬

4.1 數(shù)值模型

在本文算例中，為了簡(jiǎn)化物理問(wèn)題，采用了對(duì)稱(chēng)的二維水下爆炸模型，即正方形的TNT炸藥放置在一個(gè)正方形的四周充滿水的固壁箱內(nèi)，固壁邊界采用以上處理方法，炸藥尺寸為0.1 m×0.1 m，中心坐標(biāo)為（0 m，0 m），固壁箱尺寸為1 m×1 m，整個(gè)離散區(qū)域有100×100個(gè)粒子，其中100個(gè)TNT粒子，9 900個(gè)水粒子。初始幾何模型如圖2所示。炸藥從中心引爆后，產(chǎn)生初始沖擊波向水中傳播，并將能量傳入水中，推動(dòng)水介質(zhì)向外擴(kuò)散。TNT藥柱初始密度為1 630 kg／m3，其爆轟產(chǎn)物狀態(tài)方程取為JWL方程，即式（1），式中參數(shù)取值見(jiàn)表1；水的密度為1 000 kg／m3，狀態(tài)方程采用Mie－Gruneisen方程，即式（2）或式（3），式中參數(shù)取值見(jiàn)表2。

圖2 初始粒子分布圖

表1 JWL狀態(tài)方程中的參數(shù)及初始條件

表2 Mie－Gruneisen狀態(tài)方程中的參數(shù)

4.2 數(shù)值結(jié)果及分析

圖3為整個(gè)水下爆炸過(guò)程中，包括壁面反射之后的粒子速度分布圖。從圖中可以清楚地觀察到向外傳播的初始沖擊波、從固壁面反射回來(lái)的反射波、水的壓縮運(yùn)動(dòng)以及爆炸氣泡的膨脹和壓縮等現(xiàn)象。圖4為水下爆炸過(guò)程中具有代表性瞬間的壓力分布，本文采用了Delaunay三角化方法［11］將SPH法求解域中的離散點(diǎn)連接成背景三角形網(wǎng)格，然后用每個(gè)背景三角形網(wǎng)格3個(gè)頂點(diǎn)物理量的平均值所對(duì)應(yīng)的顏色填充該三角形區(qū)域，從而實(shí)現(xiàn)無(wú)網(wǎng)格法數(shù)值結(jié)果的云圖表征。從圖中可以看出，大約200 μs時(shí)沖擊波到達(dá)壁面，然后從壁面反射回來(lái)，反射波壓縮正在向外膨脹的氣泡，大約400 μs時(shí)，氣泡尺寸達(dá)到最大，然后開(kāi)始被壓縮，隨著時(shí)間的推進(jìn)，氣泡將不斷收縮，直至最小，此時(shí)氣泡完成第一次脈動(dòng)。這樣的向外膨脹和向內(nèi)收縮會(huì)連續(xù)出現(xiàn)多次，最終將達(dá)到一個(gè)平衡狀態(tài)，但在有限空間中，氣泡的脈動(dòng)周期將會(huì)非常短，并不像實(shí)際水下爆炸那樣會(huì)持續(xù)幾十毫秒，如圖5所示，其中氣泡半徑的大小取為當(dāng)前時(shí)刻爆炸氣體粒子區(qū)域的等效半徑。

圖3 不同時(shí)刻的粒子速度分布圖

圖4 不同時(shí)刻的粒子壓力分布圖

圖6給出了在沖擊波未到達(dá)壁面前，即自由場(chǎng)中，不同爆距處的壓力時(shí)間歷程，從圖中可以看出，在峰值壓力過(guò)后，SPH方法計(jì)算的壓力值還會(huì)有較大的擾動(dòng)，出現(xiàn)壓力的雙峰或多峰現(xiàn)象，而經(jīng)驗(yàn)公式計(jì)算的壓力則是平滑下降。試驗(yàn)所測(cè)得水下爆炸沖擊波壓力時(shí)域曲線一般是不太規(guī)則的強(qiáng)間斷曲線，且?guī)в卸喾瀣F(xiàn)象［12，13］。經(jīng)驗(yàn)公式中將壓力時(shí)歷曲線作為一個(gè)指數(shù)衰減函數(shù)處理帶有一定的近似性。從這一點(diǎn)上來(lái)說(shuō)，SPH數(shù)值模擬的這種脈動(dòng)更符合實(shí)際情況。

從圖6中還可看出，沖擊波壓力峰值隨距離的增加而減小，這與庫(kù)爾經(jīng)驗(yàn)公式是一致的。從表3中數(shù)據(jù)表明SPH方法計(jì)算結(jié)果與經(jīng)驗(yàn)公式（4）計(jì)算結(jié)果比較吻合，進(jìn)一步驗(yàn)證了SPH方法用于水下爆炸沖擊載荷計(jì)算的可行性。

5 結(jié)論

本文應(yīng)用SPH方法成功模擬了整個(gè)水下爆炸過(guò)程，并給出了不同時(shí)刻的粒子速度、壓力分布圖以及沖擊波超壓隨距離和時(shí)間的變化圖，所得結(jié)果通過(guò)與經(jīng)驗(yàn)公式的比較到達(dá)了較高的精度，為艦船水下爆炸沖擊載荷的計(jì)算提供了一種新型有效的方法。

圖5 氣泡脈動(dòng)圖

表3 不同距離處沖擊波壓力峰值計(jì)算結(jié)果比較

圖6 不同距離處的壓力時(shí)間歷程

［1］ COLAGROSSI A，LANDRINI M.Numerical simulation of interfacial flows by smoothed particle hydrodynamics［J］.Journal of Computational Physics，2003，191（2）：448-475.

［2］ 韓旭，楊剛，龍述堯.SPH方法在兩相流動(dòng)問(wèn)題中的典型應(yīng)用［J］.湖南大學(xué)學(xué)報(bào)，2007，34（1）：28-32.

［3］ LIU M B，LIU G R，LAM K Y.Smoothed particle hydrodynamics for numerical simulation of underwater explosion［J］.Computational Mechanics，2003，30（2）：106-118.

［4］ LIU M B，LIU G R，LAM K Y.Investigations into water mitigation using a meshless particle method ［J］.Shock Waves，2002，2（3）：181-195.

［5］ 韓旭，楊剛，強(qiáng)洪夫，譯.光滑粒子流體動(dòng)力學(xué)——一種無(wú)網(wǎng)格粒子法［M］.長(zhǎng)沙：湖南大學(xué)出版社，2005.

［6］ 張鎖春.光滑質(zhì)點(diǎn)流體動(dòng)力學(xué)（SPH）方法（綜述）［J］.計(jì)算物理，1996，13（4）：385－397.

［7］ 陳朗，龍新平，馮長(zhǎng)根，等.含鋁炸藥爆轟［M］.北京：國(guó)防工業(yè)出版社，2004.

［8］ 庫(kù)爾.水下爆炸［M］.羅耀杰，譯.北京：國(guó)防工業(yè)出版社，1960.

［9］ MONAGHAN J J.Simulating free surface flows with SPH.Journal of Computational Physics，1994，110：399－406

［10］ Benz W.Smoothed particle hydrodynamics：a review，NATO Workshop，Les；Arcs，F(xiàn)rance.

［11］ 田仲可.無(wú)網(wǎng)格法數(shù)值結(jié)果的云圖表征［J］.青島科技大學(xué)學(xué)報(bào)（自然科學(xué)版），2007，28（4）：332－336.

［12］ 池家春，馬冰.TNT／RDX（40／60）炸藥球水中爆炸波研究［J］.高壓物理學(xué)報(bào)，1999，13（3）：199－204.

［13］ 程素秋，寧永成，張臣，等.相似理論在水下爆炸模型試驗(yàn)中的應(yīng)用［J］.艦船科學(xué)技術(shù)，2008，30（3）：95－100.

A New Approach to the Calculation of Impulsive Load of Underwater Explosion

Chen Juan1Xu Yong-qiu2Yao Xiong－liang1Wei Sha1
1 College of Shipbuilding Engineering，Harbin Engineering University，Harbin 150001，China 2 Taizhou Maritime Safety Admimistration People's Republic of China，Taizhou 318000，China

Recently，more and more researches of underwater explosion of warships are focused on the numerical simulations by using computers.Underwater explosion consists of complicated sequence of physical processes，and involves large deformation，large inhomogeneities，moving interfaces and free surfaces.Simulation of underwater explosion problems is a big challenge for the conventional grid－based numerical methods.This paper presents the application of the meshfree smoothed particle hydrodynamics（SPH）method to simulate underwater explosion of warships.And methods are proposed to treat the wall and material interface.Simulation results show the pressure distribution，velocity distribution and the gas bubble evolution and pulsation at representative instants in the underwater explosion process.A new and effective method is provided for research of impulsive load of underwater explosion of warships.

underwater explosion；shock wave；bubble；rigid wall reflection

U661.4

A

1673－3185（2009）03－18－06

2008-12-08

國(guó)家自然科學(xué)基金資助課題（50779007）；國(guó)際科技合作基金資助課題（2007DFR80340）；高等學(xué)校博士學(xué)科點(diǎn)專(zhuān)項(xiàng)科研基金資助課題（50809018）

陳 娟（1985－），女，碩士研究生。研究方向：船舶與海洋工程結(jié)構(gòu)動(dòng)力學(xué)。E-mail：chenjuan＠hrbeu.edu.cn

許永秋（1966-），男，工程師。研究方向：海事管理

姚熊亮（1963－），男，教授，博士生導(dǎo)師。研究方向：船舶與海洋工程結(jié)構(gòu)動(dòng)力學(xué)，水下沖擊載荷作用下結(jié)構(gòu)響應(yīng)。E-mail：saibei8411@163.com