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

陳 娟 許永秋 姚熊亮 位 莎

1哈爾濱工程大學 船舶工程學院，黑龍江 哈爾濱150001 2浙江省臺州市海事局，浙江 臺州318000

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

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

1哈爾濱工程大學 船舶工程學院，黑龍江 哈爾濱150001 2浙江省臺州市海事局，浙江 臺州318000

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

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

1 引言

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

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

2 水下爆炸理論基礎

水下爆炸是指在極短時間內，在水下的極小體積內或面積上發生極大能量轉換的過程。水下爆炸大體可分為3個階段：裝藥的爆轟、沖擊波的產生和傳播、氣泡的形成和脈動。當炸藥在水中爆炸時，其周圍介質直接受到具有高溫、高速、高壓的爆炸產物作用。在裝藥和介質的界面處，爆炸產物以極高的速度向周圍擴散，強烈地壓縮著相鄰的水，使其壓力、密度、溫度突躍式地升高，形成初始沖擊波。隨著沖擊波在水中的傳播，爆炸產物在水中以氣泡的形式存在并不斷膨脹與壓縮，形成氣泡脈動現象。

2.1 狀態方程

狀態方程用來定義固體或流體在各種不同狀態下的壓力和密度以及比內能之間的函數關系，正確選取狀態方程中的參數對于計算結果至關重要。狀態方程一般采用半經驗半理論的公式，方程中的主要參數由試驗確定。

1）TNT爆炸產物狀態方程

本文中TNT爆炸產物采用標準JWL狀態方程［7］，JWL狀態方程是典型的動力學狀態方程，它是一種不顯含化學反應、由實驗方法確定參數的經驗狀態方程，能比較精確地描述爆轟產物的膨脹驅動做功，其形式為：

式中：p為爆轟產物的壓力；η為爆炸氣體密度與初始炸藥密度的比值，即η＝；A，B，R1，R2和ω為與炸藥狀態有關的常數；e為TNT單位質量的內能。

2）水的狀態方程

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

在膨脹狀態下水的壓力為：

式中：ρ0為水的初始密度；η＝為水擾動前后的密度比；μ＝η－1，當μ＞0時，水處于壓縮狀態，當μ＜0時，水處于膨脹狀態；C為聲速；γ0為Gruneisen系數；a為體積修正系數；S1，S2，S3為擬合系數。

2.2 水下爆炸沖擊波的經驗計算公式

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

式中：pm為沖擊波峰值壓力，MPa；W為裝藥量，kg；R為爆距，m；k為實驗系數，取52.4；α為壓力衰減系數，對標準TNT炸藥，取1.13。

3 SPH數值計算方法

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

假設爆炸氣體和周圍水介質是無粘性的，且整個水下爆炸過程被認為是絕熱的，因此采用Euler方程作為控制方程。應用SPH核近似和粒子近似可得到以下一系列SPH方程，且在動量方程和能量方程中引入了Monaghan型人工粘性∏ij［9］，它不僅將動能轉化為熱能，提供了沖擊波波面必不可少的耗散，而且防止計算過程中粒子之間的非物理穿透。

Monaghan型人工粘性∏ij表達式為：

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

對上式進行SPH近似，可得到以下形式：

下一個時間步長的光滑長度變為：

3.2 固壁邊界處理

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

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

圖1中粗直線表示剛性壁，圓表示粒子i的支持域。設j位于i的支持域內，若實粒子j關于剛性壁對稱的鏡像粒子j′也在i的支持域中，則粒子i鄰域內的點也包括鏡像粒子j′，如果與k對應的鏡像粒子k′落在粒子i的支持域外，則計算時不予考慮。由于粒子i也在其支持域內，因此其鏡像粒子i′同樣需要考慮。

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

3.3 交界面處理

在SPH中，物質交界面處理是一個非常關鍵的問題，因為SPH方法具有拉格朗日性質和粒子性質，在整個演變過程中，來自不同介質的相互接觸的粒子可能會隨著運動而分離，甚至有可能不再成為相鄰粒子。在水下爆炸模擬中的物質交界面處理更為復雜，因為高壓高能的爆炸氣體和水之間的激烈相互作用會使粒子的運動相當自由，本文采用的方法是，離得近但材料不同的兩粒子可以認為是相鄰粒子，這樣能夠減少在邊界附近的粒子缺陷問題，但為了防止粒子的非物理穿透或摻雜問題，在交界面附近的不同粒子之間，當其趨向于穿透時，即：pe＝≥1，施加一個排斥力PBij，其表達式為：

4 水下爆炸過程的數值模擬

4.1 數值模型

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

圖2 初始粒子分布圖

表1 JWL狀態方程中的參數及初始條件

表2 Mie－Gruneisen狀態方程中的參數

4.2 數值結果及分析

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

圖3 不同時刻的粒子速度分布圖

圖4 不同時刻的粒子壓力分布圖

圖6給出了在沖擊波未到達壁面前，即自由場中，不同爆距處的壓力時間歷程，從圖中可以看出，在峰值壓力過后，SPH方法計算的壓力值還會有較大的擾動，出現壓力的雙峰或多峰現象，而經驗公式計算的壓力則是平滑下降。試驗所測得水下爆炸沖擊波壓力時域曲線一般是不太規則的強間斷曲線，且帶有多峰現象［12，13］。經驗公式中將壓力時歷曲線作為一個指數衰減函數處理帶有一定的近似性。從這一點上來說，SPH數值模擬的這種脈動更符合實際情況。

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

5 結論

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

圖5 氣泡脈動圖

表3 不同距離處沖擊波壓力峰值計算結果比較

圖6 不同距離處的壓力時間歷程

［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方法在兩相流動問題中的典型應用［J］.湖南大學學報，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］ 韓旭，楊剛，強洪夫，譯.光滑粒子流體動力學——一種無網格粒子法［M］.長沙：湖南大學出版社，2005.

［6］ 張鎖春.光滑質點流體動力學（SPH）方法（綜述）［J］.計算物理，1996，13（4）：385－397.

［7］ 陳朗，龍新平，馮長根，等.含鋁炸藥爆轟［M］.北京：國防工業出版社，2004.

［8］ 庫爾.水下爆炸［M］.羅耀杰，譯.北京：國防工業出版社，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，France.

［11］ 田仲可.無網格法數值結果的云圖表征［J］.青島科技大學學報（自然科學版），2007，28（4）：332－336.

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

［13］ 程素秋，寧永成，張臣，等.相似理論在水下爆炸模型試驗中的應用［J］.艦船科學技術，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

國家自然科學基金資助課題（50779007）；國際科技合作基金資助課題（2007DFR80340）；高等學校博士學科點專項科研基金資助課題（50809018）

陳 娟（1985－），女，碩士研究生。研究方向：船舶與海洋工程結構動力學。E-mail：chenjuan＠hrbeu.edu.cn

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

姚熊亮（1963－），男，教授，博士生導師。研究方向：船舶與海洋工程結構動力學，水下沖擊載荷作用下結構響應。E-mail：saibei8411@163.com