基于SPH方法的淺水爆炸研究

劉翠丹，明付仁，叢 剛，王 超

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

0 引言

爆炸沖擊波攻擊艦船和潛艇屬于水下爆炸在國防上的應用，水下爆炸產(chǎn)生的高壓沖擊波會使艦船發(fā)生毀滅性破壞。淺水爆炸由于邊界效應，沖擊波壓力十分復雜，對淺水環(huán)境中的艦船打擊作用顯著加強，因此研究淺水爆炸中沖擊波的特性對艦船受攻擊時的沖擊響應預報及艦船防護結構設計具有重要意義。

對于沉底和自由面爆炸，已經(jīng)有不少探索，但是對于淺水爆炸來說，目前可參考的資料還較少，這是由于淺水爆炸要考慮水底及水面雙重的影響，情況較復雜。在試驗方面，劉文華［1］等通過單個球形裝藥淺水爆炸研究了幾種參數(shù)對水中沖擊波峰值壓力的影響，探索了邊界條件影響下沖擊波的變化規(guī)律;高永軍［2］等對淺層水中爆炸沖擊波參數(shù)和空間流場的壓力分布進行分析，得出試驗范圍內(nèi)的沖擊波壓力的經(jīng)驗公式;張鵬翔［3］等分析淺水爆炸的水面切斷現(xiàn)象對沖擊波的正壓作用時間和比沖量的影響，對實驗數(shù)值進行回歸，得到精度較高的回歸公式。在數(shù)值模擬方面，顧文彬［4］等數(shù)值模擬研究水底和水面效應對沖擊波的影響;孫百連［5］等對沉底2個裝藥爆炸進行數(shù)值模擬，將其與單個裝藥爆炸和無限水域2個裝藥同時爆炸進行了對比，研究水底水面對沖擊波的影響。在理論方面，張鵬翔［6］等探討了淺水爆炸中的沖擊波切斷現(xiàn)象，主要分析了切斷現(xiàn)象的特性，產(chǎn)生機理以及對沖擊波的影響，對兩類切斷現(xiàn)象進行了探討。上述文獻中大多數(shù)采用試驗的方法進行分析，還有少數(shù)采用有限元分析。而傳統(tǒng)有限元方法分析爆炸問題存在弊端，特別是大變形、運動交界面問題，經(jīng)常導致計算崩潰甚至無法計算。近年來引起關注的是新一代計算方法——無 網(wǎng) 格 法。 SPH (Smoothed Particle Hydrodynamics)方法作為一種拉格朗日形式的無網(wǎng)格粒子法具有自適應、無網(wǎng)格、粒子形式等特性，在大變形的流體動力學問題上具有獨特優(yōu)勢，適合于解決水下爆炸產(chǎn)生的大變形問題。

本文通過SPH方法進行數(shù)值模擬，再現(xiàn)了淺水爆炸中的基本過程，分析沖擊波的壓力峰值和沖量特性、研究爆深和泥沙密度等參數(shù)對淺水爆炸的影響，旨在揭示近邊界爆炸沖擊波傳播特征，為艦船水下爆炸的載荷確定、結構響應預報及防護結構設計提供參考。文中使用的程序?qū)ψ杂蓤鲋械乃卤M，并將結果與經(jīng)驗公式對比，驗證了程序的精度，限于篇幅，這里不予呈現(xiàn)。

1 基本理論

1.1 SPH控制方程

拉格朗日描述下的流體控制方程可寫作一系列偏微分方程，也就是著名的Navier-Stokes方程，其中包括質(zhì)量守恒方程、動量守恒方程和能量守恒方程。SPH方法是應用核近似和粒子近似對N-S方程進行離散，得到以下SPH形式的控制方程［7］:

以上分別為連續(xù)性方程、動量守恒方程、能量守恒方程。Π采用Monaghan型人工粘度項［8］，下標i，j表示相互作用的1對粒子，上標α，β表示沿坐標軸的方向。

1.2 狀態(tài)方程

1.2.1 TNT狀態(tài)方程

TNT爆炸產(chǎn)生的爆轟產(chǎn)物采用Jones-Wilkins-Lee狀態(tài)方程［9］:

式中:η為爆轟產(chǎn)物密度與初始炸藥密度的比值，即η=ρ/ρ0;A，B，R1，R2和ω為與炸藥狀態(tài)有關的常數(shù)。

1.2.2 水的狀態(tài)方程

水在爆炸沖擊波傳播階段采用Mie-Grüneisen狀態(tài)方程［10］，其具體形式取決于水的狀態(tài)。

在壓縮狀態(tài)下水的壓力為

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

式中:μ為壓縮比;e為單位質(zhì)量內(nèi)能;μ=η-1，當μ＞0時，水處于壓縮狀態(tài)，當μ＜0時，水處于膨脹狀態(tài);C0為聲速;γ0為 Grüneisen系數(shù);a為體積修正系數(shù);S1;S2;S3為試驗擬合系數(shù)。

1.2.3 泥沙的狀態(tài)方程

泥沙的狀態(tài)方程引用文獻 ［11］中的經(jīng)驗形式。

多項式狀態(tài)方程是流體物態(tài)方程的經(jīng)驗形式，常用其擬合Mie-Grüneisen物態(tài)方程。多項式狀態(tài)方程由下式給出

式中:a1，a2，a3，b0，b1，b2為物態(tài)方程經(jīng)驗系數(shù);μ，e與前同。而Grüneisen系數(shù)表示為

1.3 淺水爆炸的基本理論

同時受水面水底影響的裝藥水中爆炸定義為淺層水中爆炸［12］。淺水爆炸沖擊波主要受水面、水底2種界面的影響，是一個涉及水、目標運動及邊界效應的復雜爆炸動力學問題。

稀疏波，該波經(jīng)過一段時間后到達測點，使測點處的壓力受到削弱，沖擊波在水底反射，當反射沖擊波遇上入射波，會使得測點處的壓力峰值比無限場中大。水面反射波到達測點，與入射沖擊波或反射沖擊波作用，使得合成壓力低于水介質(zhì)的靜壓力，表現(xiàn)在壓力曲線上波形被削去一截，稱作水面切斷現(xiàn)象［3］。

水底對爆炸效果的影響是多方面的。爆炸沖擊波到達水底后，產(chǎn)生反射沖擊波在水中傳播，同時產(chǎn)生縱波和橫波。當這2種彈性波與水底的分界面作用時產(chǎn)生瑞利波。瑞利波在水底的不規(guī)則反射形成前驅(qū)波，使壓力增加10%～20%。反射波追上入射波，二者融合形成馬赫波。與無限水域相比馬赫波使壓力增加3%～7%，可見水底對沖擊波的影響不可小視［12］。隨著測點高度或測點距爆心距離增加，水底影響減?。?］。裝藥量和水底材料一定時，爆深不同，測點壓力有很大不同。測點相同，隨著爆深的增加，相同比例距離上的峰值壓力增加［13］。

2 淺水爆炸數(shù)值模擬

2.1 數(shù)值模型

本文采用SPH方法模擬淺水爆炸，利用Fortran編程建模，模型圖見圖1。

圖1 淺水爆炸模型Fig.1 Shallow water explosion model

圖中上面為水域，下面為泥沙水底，正方形TNT炸藥位于水域的中央。模型中常量主要有:水深a=0.600 m，水底泥沙尺寸b=1.500 m，c=3.000 m，裝藥 TNT質(zhì)量 m=0.360 kg，直徑 d=0.045 m。模型中的變量有炸藥中心距自由面距離高度h、測點距爆心的水平距離e和垂向距離f以及泥沙的密度ρ。淺水爆炸模型中變量h=0.300 m，ρ=1 800 kg/m3。為了研究淺水爆炸時自由面和泥沙水底的影響，本文將其計算結果和相同藥量無限水域爆炸的計算結果進行對比。無限水域爆炸模型長3.0 m，寬2.1 m，炸藥位于水域中心，TNT尺寸與淺水爆炸計算模型相同。

2.2 基本物理現(xiàn)象

炸藥水中爆炸時，在裝藥和介質(zhì)的界面處，爆炸產(chǎn)物以極高的速度向周圍擴展，強烈地擠壓周圍的水介質(zhì)，使其壓力、密度、溫度突然升高，形成初始沖擊波。

當水域比較淺時，底部效應與自由面效應都不能忽略，自由面效應和底部效應共同對直達波產(chǎn)生作用。圖2顯示了不同時刻下淺水爆炸計算模型的壓力云圖和粒子分布圖，從圖中可知，TNT引爆后，會產(chǎn)生壓力極高的強沖擊波在水域中傳播，在70 μs時是沖擊波到達邊界之前的云圖，此時沖擊波還沒有受到邊界的影響;在130 μs時，沖擊波傳播至自由面處，反射一列稀疏波，使稀疏波的影響區(qū)域壓力降低，由于自由面附近爆炸會削弱氣泡對舷側的攻擊效果［14］，因此可以得出結論，若武器在淺水中靠近自由面爆炸時，會削弱爆炸產(chǎn)物的攻擊效果。此時沖擊波傳到水底在泥沙中產(chǎn)生圓形的沖擊波向外傳播，泥沙中的沖擊波壓力高于水中的沖擊波壓力。至200 μs時，2種反射波相向傳播，即將相遇;至380 μs時，2種反射波相遇，稀疏波的影響范圍逐漸增大，涵蓋近自由面的大部分區(qū)域，此時，水下爆炸沖擊波的高壓高能作用使自由面發(fā)生了水冢現(xiàn)象。在水和泥沙的交界面處，由于自由面反射的稀疏波致使泥沙中沖擊波壓力迅速降低，逐漸和水中壓力相同。同時沖擊波在水底的反射沖擊波遇上入射波，在近水底處產(chǎn)生前驅(qū)波，此波與入射波疊加，使測點處壓力變大。至540 μs時，水冢不斷地增大，由于爆炸沖擊的作用，鄰近自由面的水粒子產(chǎn)生凸起的現(xiàn)象，隨著爆炸產(chǎn)生的氣體接近底部泥沙，氣體的非球形就逐漸變得明顯，這里的現(xiàn)象與文獻 ［14］中的結論一致，即隨著氣泡距離壁面越近，氣泡非球狀出現(xiàn)的越早［14］。

圖2 淺水爆炸粒子分布圖及壓力云圖Fig.2 Particle distribution and explosion pressure pictures of shallow water explosion

2.3 壓力及沖量特性

圖3顯示了e=0.40 m，f=0.25 m，f=0.01 m，f=-0.24 m時，淺水爆炸和無限水域爆炸的壓力時程對比曲線。

圖3 不同測點的壓力時程曲線Fig.3 Pressure process curve of different point

當f=0.25 m時，由于測點位置靠近自由面，此時沖擊波壓力主要受自由面影響，受底部影響較小，自由面使沖擊波產(chǎn)生“切斷”效應，即把沖擊波的波尾截掉，減小沖擊波的沖量和能量，并使這部分能量轉(zhuǎn)化為流體向上運動的動能。

當f=-0.24 m時，由于測點位置靠近水底，此時沖擊波壓力主要受底部影響，從0.2 ms開始，水底反射波傳到測點，淺水爆炸壓力值大于無限水域壓力值。這是由于泥沙的反射作用使沖擊波變強。到0.45 ms，自由面反射波到達測點，使得壓力值降低，小于無限水域的壓力值。

當f=0.01 m時，測點位于爆心的水平左側，受到自由面反射波和水底反射波的雙重作用，壓力曲線既出現(xiàn)了“截斷”效應，又出現(xiàn)了壓力峰值交錯。

沖擊波沖量是丈量沖擊波強弱的標志，表1中給出了2種工況下的測點沖量。

表1 測點沖量Tab.1 Point impulse

自由面會使沖擊波的沖量減少，而底部會使沖擊波的沖量增大，分析上表可知，當f=0.25 m時，測點以自由面影響為主，沖量減少了36%;當f=-0.24 m時，測點以底部影響為主，沖量增加了14%;當f=0.01 m時，測點同時受到自由面和底部的影響，作用效果存在一定的抵消。

2.4 不同爆深的影響

不同爆深對淺水爆炸的壓力峰值有影響，為了研究該問題，本文分別取爆深h=0.1 m，h=0.3 m，h=0.5 m，其余的參量都與淺水爆炸模型相同。經(jīng)過程序計算，得到對比云圖 (圖4)和壓力時程曲線(圖5)。

從圖4可以看出，h越大，淺水爆炸產(chǎn)生的水冢現(xiàn)象就越不明顯，若爆深足夠大，則不會出現(xiàn)水?，F(xiàn)象。

圖5顯示了3種爆深的壓力時程曲線。這里由于爆深差別較大，3個測點集中選擇在爆心的左側附近。

對于不同爆深下相對位置相同的測點，當h=0.1 m時，爆心靠近水面，壓力曲線主要受到近水面反射波的影響，出現(xiàn)明顯的截斷，壓力峰值也迅速下降，壓力曲線尾部出現(xiàn)截斷效應。當h=0.5 m時，爆心接近水底，壓力曲線主要受底部影響，脈寬變大，作用時間變長。而h=0.3 m時，測點同時受到自由面和底部的影響，壓力曲線先是受水面影響出現(xiàn)截斷現(xiàn)象，然后受底部影響延長了作用時間，并且在尾部出現(xiàn)壓力波動的現(xiàn)象。

對于同一爆距下的3條曲線，當f=-0.09 m時曲線波動的幅度較其他2條明顯。這是由于此時測點距水底較近，流場區(qū)域內(nèi)受底部泥沙影響較大。但是壓力峰值比前二者中爆炸的峰值小。這是由于泥沙土壤在高壓爆炸產(chǎn)物的作用下，迅速變形成坑，消耗了一部分爆炸能量，引起水底附近沖擊波壓力的削弱，所以此時水底附近的壓力要比前二者小。

沖擊波的沖量是壓力曲線與坐標軸圍成的面積，分析圖5可以看出，隨著h的降低，沖擊波的沖量減少，表2給出了各工況下的沖量值，從數(shù)值上也說明了這點。自由面會使沖擊波的沖量減少，而底部會使沖擊波的沖量增大。分析表2可知，當h=0.1 m時，測點以自由面影響為主，沖量較小;當h=0.5 m時，測點以底部影響為主，沖量較大;當h=0.3 m時，測點同時受到自由面和底部的影響，作用效果存在不同程度的抵消?？偨Y上述分析可知，隨著h的增加，沖量會增加。

表2 測點沖量Tab.2 Point impulse

2.5 不同底部泥沙密度的影響

在淺水爆炸中，水底泥沙的密度也是影響壓力的因素，通常泥沙的密度范圍是1 500～2 700 kg/m3。為了研究該問題，本文分別選取變量泥沙密度為1 500 kg/m3，1 800 kg/m3和2 400 kg/m3這3 種工況，經(jīng)過編程計算，得到對比云圖(圖6)和壓力時程曲線(圖7)。

從圖6可以看出，泥沙密度對壓力值的影響，隨著泥沙密度的增大，壓力峰值也相應的增大，增幅約為1.1倍。圖6的現(xiàn)象是由爆轟氣體(入射波)、水面對沖擊波的反射波以及水底對沖擊波的反射波三者相互作用形成的，觀察圖6可知，隨著泥沙密度的增大，使得3種波相互作用的現(xiàn)象更加明顯，當ρ=2 400 kg/m3時，藥包下方的壓力明顯減小了，與水接近，這是由于泥沙密度增大使水底反射波的壓力變大，其與爆轟產(chǎn)生的氣體波作用時相互抵消的結果。

圖7中分別給出了3個測點的壓力時程曲線。

改變泥沙密度主要影響的是水底反射波的特性，從而影響整個沖擊波。

當f=0.2 m時，由于測點處于近水面，離水底較遠，幾乎不受水底反射波的影響，只有爆轟氣體和近水面反射波作用，因此3條曲線輪廓幾乎重合。而且可以看出，受近水面反射波的影響，壓力曲線出現(xiàn)截斷現(xiàn)象。

當f=0 m時，從180 μs開始，爆轟氣體傳到測點，開始產(chǎn)生壓力，此時測點位于水域中間，同時受水面和水底的影響。到340 μs，3種波同時到達測點，曲線從此時開始出現(xiàn)不同的現(xiàn)象。設p1500，p1800，p2400為某一時刻的壓力值，從圖7中可以看出，當t＞340 μs時，p2400＞p1800＞p1500，水底反射波使壓力值變大，隨著水底泥沙密度的增大，壓力增大的比例也相應變大。

當f=-0.2 m時，從200 μs開始，爆轟氣體傳到測點，開始產(chǎn)生壓力，此時測點靠近水底，以底部影響為主，脈寬變大。到260 μs時，水底反射波和爆轟氣體同時到達測點，曲線開始出現(xiàn)不同的現(xiàn)象。從圖中可以看出，受水底反射波影響，當 260 μs＜ t＜520 μs時，p2400＞p1800＞p1500，隨著水底泥沙密度的增大，壓力值增大。到520 μs時，水面反射波到達測點，測點此時受3種波作用，從壓力曲線可以看出，測點的壓力受到削弱。與無限水域相比，無論密度多大，受淺水效應的影響，到540 μs時壓力都等于0。

3 結語

淺水爆炸問題目前主要集中在試驗研究，而本文針對此類問題采用無網(wǎng)格SPH方法對其進行數(shù)值模擬，從不同爆深、不同底部泥沙密度等幾個方面進行了研究和思考，得出如下結論:

1)淺水爆炸中截斷與增強效應共存，這是由于水面和底部泥沙共同作用的結果。一方面水面效應能夠反射稀疏波，減小沖擊波的沖量和能量，并使這部分能量轉(zhuǎn)化為流體向上運動的動能，產(chǎn)生截斷效應和水冢現(xiàn)象;另一方面，爆炸沖擊波與水底沖擊作用，到達底部發(fā)生反射，使水底快速變形，壓實，由此產(chǎn)生的底部效應能夠增大沖擊波的沖量和能量，使壓力作用時間增長，表現(xiàn)為增強效應。

2)不同爆深時，淺水中壓力和沖量受水底和水面影響程度不同。當爆深較小時，受水面影響較大，壓力曲線出現(xiàn)截斷的效應。當爆深較大時，壓力曲線主要受底部影響，作用脈寬變大，峰值增大。在爆深從小逐漸增大的過程中，受水面和水底效應共同影響，沖量增大，壓力曲線分別經(jīng)歷峰值和脈寬的削弱、峰值和脈寬增大的過程，并且在尾部出現(xiàn)壓力波動的現(xiàn)象。隨著爆深的增加，水冢現(xiàn)象逐漸變得不明顯直至消失。

3)不同底部泥沙密度時，淺水中壓力與沖量呈增強趨勢。隨著泥沙密度的增加，壓力峰值增大，3種泥沙密度在380 μs增幅約為1.1倍。隨著泥沙密度的增大，底部反射波的增強效應逐漸大于水面反射波的截斷效應，3種波(入射波、水面反射波以及水底反射波)相互作用的現(xiàn)象更加明顯。

［1］劉文華，羅松林，顧文彬，等.單個球形裝藥淺層水中爆炸沖擊波特性的研究［J］.工程爆破，1999，5(3):1-5.LIU Wen-hua.LUO Song-lin.GU Wen-bin，et al.Research on specific properties of shock wave from explosion of single spherical charge in shallow water［J］.Engineering Blasting，1999，5(3):1-5.

［2］高勇軍，陳小波，等.單個球形裝藥淺層水中爆炸的試驗研究［J］.工程兵工程學院學報，1998，14(1):96-100.GAO Yong-jun.CHEN Xiao-bo，et al.Experimentation and research on a spherical-charge explosion under shallow water［J］.Journal of Nanjing Engineering Institute，1998，14(1):96-100.

［3］張鵬翔，顧文彬，等.淺層水中爆炸水面切斷現(xiàn)象的實驗研究［A］.第七屆爆破學術會議論文集［C］.成都:2001:353-357.ZHANG Peng-xiang.GU Wen-bin，et al.Experimental researchesofcutoffphenomenon by watersurface subjected to shock waves generated by explosion in shallowlayer water［A］.The 7th blasting academic conference proceedings［C］.Chengdu:2001:353-357.

［4］顧文彬，孫百連，等.淺層水中沉底爆炸沖擊波相互作用數(shù)值模擬［J］.解放軍理工大學學報，2003，4(6):64-68.GU Wen-bin.SUN Bai-lian，et al.Numerical simulation of explosive shock wave interaction in shallow-layer water［J］.Journal of PLA University of Science and Technology，2003，4(6):64-68.

［5］孫百連，顧文彬，等.淺層水中沉底的兩個裝藥爆炸的數(shù)值模擬研究［J］.爆炸與沖擊，2003，23(5):460-465.SUN Bai-lian.GU Wen-bin，et al.Numerical simulation of explosion shock wave interaction in shallow-layer water［J］.Explosion and Shock Waves，2003，23(5):460-465.

［6］張鵬翔，顧文彬，等.淺層水中爆炸沖擊波切斷現(xiàn)象淺探［J］.爆炸與沖擊，2002，22(3):221-228.ZHANG Peng-xiang.GU Wen-bin，et al.Discussions of blasting shock waves cutoff in shallow-layer water［J］.Explosion and Shock Waves，2002，22(3):221-228.

［7］LIU G R，LIU M B.光滑粒子流體動力學——一種無網(wǎng)格粒子法［M］韓旭，楊剛，強洪夫，譯.長沙:湖南大學出版社，2005.1-432.LIU G R，LIU M B.Smoothed particle hydrodynamics:A mesh free particle method［M］.Changsha:Hunan University Press，2005.1-432.

［8］MONAGHAN JJ，On the problem of penetration in particle methods［J］.Journal of Computational physics，1989，82:1-15.

［9］DOBRATZ B M.LLNL Explosive Handbook［R］.UCRL-52997，Lawrence Livermore National Laboratory，Livermore，CA.1981.

［10］STEINBERG D J.Spherical explosions and the equation of state of water［R］.UCID-20974，Lawrence Livermore National Laboratory，Livermore，CA.1987.

［11］湯文輝，張若棋.物態(tài)方程理論及計算概論(第二版)［M］.北京:高等教育出版社，2008.278-300.TANG Wen-hui.ZHANG Ruo-qi.Overview of theory and calculation of the forms equation.［M］.Beijing:Higher Education Press，2008.278-300.

［12］顧文彬.葉序雙，等.淺層水中爆炸水底影響的試驗研究［J］. 解放軍理工大學學報，2001，2(2):55-58.GU Wen-bin.YE Xu-shuang，et al.Experimental studies of bottom influence in shallow-layer water explosion［J］.Journal of PLA University of Science and Technology，2001，2(2):55-58.

［13］顧文彬.葉序雙，等.界面對淺層水中爆炸沖擊波峰值壓力影響的研究［J］.解放軍理工大學學報，2001，2(5):61-63.GU Wen-bin.YE Xu-shuang，etal.Peak pressure investigation of exploding wave influenced by interfaces in shallow-layer water［J］.Journal of PLA University of Science and Technology.2001，2(5):61-63.

［14］張阿漫，王詩平，等.不同沖擊環(huán)境下氣泡脈動特性實驗研究［J］.力學學報，2011，43(1):71-83.ZHANG A-man，WANG Shi-ping，et al.Experimental study on bubble pulse features under different circumstances［J］.Chinese Journal of Theoretical and Applied Mechanics，2011，43(1):71-83.