基于SPH方法的水下爆炸初始爆轟過程研究

1 引 言

對于水下爆炸問題的數值模擬，國內外開展了許多研究，炸藥爆轟是極其短暫且復雜的過程，并且爆炸沖擊效應的計算一直是計算力學中的難題，所以炸藥的理論解法只局限于一些簡單情形，而試驗研究則需承擔風險并且費用昂貴，但有時試驗也不一定能捕獲與爆轟相關的一些物理現象。隨著計算機技術和計算方法的不斷發展，許多數值模擬方法應用到了炸藥的爆轟研究中。國內的爆轟數值模擬主要以有限差分為主，有限元法正處在起步階段，而要開發一個具有實用價值的程序系統不僅需要大量的人力、物力，更需要有相當數量的系統實驗作為儲備，以提供充足的物理參數，這一點國內是缺乏的。炸藥在爆炸過程中會產生大變形，令網格發生嚴重變形，從而導致無效的極小時間步長的產生，有時甚至導致計算崩潰，所以基于網格的有限元法或有限差分法數值模擬的方法在模擬分析炸藥爆轟時遇到了很多困難。而SPH方法[1,2]具有無網格性質和拉格朗日粒子特性，它應用離散化的粒子來表示物質，流體粒子的運動歷程就能很自然地被捕獲。SPH方法的無網格性能能夠克服在計算機中與大變形相關的困難，再加上顯式算法，這一切使SPH方法非常適宜處理發生在炸藥爆炸的極短瞬時具有的大變形和高度非均勻的動力學極端情形。

2 爆轟理論基礎

2.1 炸藥爆轟原理

炸藥爆炸所產生的化學反應過程能迅速將初始炸藥轉變為高壓氣態產物。典型的炸藥爆炸由兩部分組成：一是爆轟過程；二是爆轟后向外界周圍媒質釋放高壓氣體產物的膨脹過程。在爆轟過程中，炸藥的反應率實際上是無窮大的，若達到化學平衡，則稱爆轟處于穩定狀態。在穩態中，爆轟波的速度D為常量，一維炸藥爆轟過程如圖1所示。20世紀初，Chapman D L和Jouquet E等人不考慮爆轟的化學動力學過程，單純從流體動力學角度出發，將爆轟波看作未反應物質與反應物質之間的間斷面，運用質量、動量、能量守恒方程研究問題，并且提出爆轟波穩定傳播條件，從而發展成為爆轟波C-J理論。沖擊波波面以爆炸速度D沿高能炸藥向外推進，并壓縮炸藥，使高能炸藥由初始狀態點(p0,ρ0)沿著炸藥的Hugoniot曲線到達新的狀態點，且壓力上升到p1,如圖2所示。反應完成后，系統沿著氣體產物的曲線到達狀態點(p,ρ)。此時，氣體產物沿著通過點(p,ρ)的等熵線膨脹。

圖1 一維炸藥爆轟過程示意圖

圖2 Hugoniot曲線和Rayleigh曲線

2.2 爆轟控制方程

由于爆炸和膨脹的速度相當高，故可假設氣體生成物是非粘性的且爆炸過程是絕熱的。因此可利用如下的歐拉方程與恰當的狀態方程來模擬爆炸過程。

(1)

2.3 爆轟產物狀態方程

在控制方程中，質量守恒方程、動量守恒方程和能量守恒方程比較固定，而狀態方程則多種多樣，狀態方程一般采用半經驗半理論的公式，方程中的主要參數由試驗決定。本文中爆轟壓力、單位體積內能E和相對體積V的關系采用Jones-Wilkins-Lee(JWL)狀態方程[3,4]，JWL狀態方程是典型的動力學狀態方程，它是一種不顯含化學反應，由實驗方法確定參數的經驗狀態方程，能比較精確地描述爆轟產物的膨脹驅動做功。

JWL狀態方程于1965年由美國勞倫斯利弗莫爾國家實驗室的Lee E L在Jones和Wilkins的工作基礎上提出的。其形式為:

(2)

上式中，V=ρ0/ρ為爆轟產物的相對比容，其中ρ0為初始炸藥的密度,ρ為爆轟產物的密度；P為爆轟產物的壓力；E0為比熱力學能；A、B、R1、R2、ω為JWL狀態方程的5個待定的參數，這5個參數的確定依賴于圓筒實驗，該試驗測定金屬圓筒壁在待測炸藥爆轟產物驅動下的膨脹過程。

以上由質量、動量、能量3個守恒定律建立了爆轟波的3個基本方程式，再加上爆轟氣體產物的狀態方程，就有了4個方程，但要求解5個參數p1,ρ1,u1,e1和D，就必須再建立一個方程，Chapman D L和Jouguet E根據爆轟波在炸藥中的穩定傳播條件，建立了第5個方程，即爆轟波穩定傳播的C-J條件。

3 SPH模擬水下爆炸初始爆轟的過程

3.1 人工粘性

已有很多學者對人工粘性的形式提出了很多的建議，但是最常用的方法還是采用摩納罕提出的形式：

(3)

3.2 光滑長度

在SPH方法中，光滑長度h的選取非常重要，它直接影響計算的效率和結果的精度。若h太小，在以kh為半徑的支持域內沒有足夠的粒子對所求粒子施加力的作用，直接導致較低的精度;若h太大，則可能將粒子的所有細節信息和局部特性忽略，同樣影響精度。并且h太大，支持域內粒子過多，也降低了計算效率。

有很多方法可以動態調整光滑長度h,使得最近粒子數目保持相對穩定。但幾乎所有形式的光滑長度都是與所研究問題直接相關，很少有廣義的光滑長度。在水下爆炸爆轟模擬中，問題域內會產生極大的密度不均勻性，并且在爆炸過程中粒子運動非常激烈。本文采用了劉謀斌博士推導出的在時間和空間上展開的自適應模式的光滑長度。

最簡單的方法是按照平均密度來修正光滑長度h：

h=h0(ρ0/ρ)1/d

(4)

式中，h0,ρ0分別為初始的密度和光滑長度；d是維數。

3.3 近鄰粒子搜索

在SPH方法中，一個非常重要也很耗時的一個步驟就是近鄰粒子的搜索，因為每一個粒子在其以kh為半徑的支持域內的所有粒子，在粒子近似過程中都要被使用，且粒子的相對位置都是不斷改變的，所以，每經歷一個時間步，都需要找出每一個粒子在該時刻的近鄰粒子。本文所用的是最簡單的直接搜索法。如圖3所示，對一個給定的粒子i，計算它與整個問題域內其它所有N-1個粒子的距離(N為問題域內所有粒子的總數)，若其距離小于給定粒子支持域半徑kh，這些粒子即是粒子i的近鄰粒子，同時i也是這些粒子的近鄰粒子。這個過程在每個時間步對問題域內所有粒子都要進行一次，所以需消耗大量時間，一般只在粒子數較少的時候使用。

圖3 直接搜索法

3.4 SPH離散公式

同時將摩納罕提出的人工粘性Πij加入到SPH等式中的壓力項中，則爆炸控制方程的SPH離散形式為：

(5)

3.5 時間積分

SPH的偏微分方程組的數值解法主要有蛙跳法、四階龍格庫塔法和改進歐拉法。在實際應用中，由于蛙跳法對存儲的需求量低，而且計算效率高，因此本文使用蛙跳法[5]進行數值積分。

(6)

對時間步長的取值需要特別注意，因為步長取得過大，整個粒子系統就會發散，步長取得太小，會使計算量大大增加，計算精度也得不到明顯改善。為滿足蛙跳法的穩定條件，時間步長同時考慮介質的聲速和粘性耗散，時間步長[6]的形式為：

Δt=min[0.3h/(hi·vi+ci+1.2

(αΠci+βΠ|·vi|))]

(7)

4 SPH爆轟數值的模擬驗證

本節將應用SPH方法對一維板條TNT的水下爆炸初始爆轟過程[7]進行模擬計算，在該模型中，采用0.2 m長的板條TNT炸藥，從中間起爆向兩端爆燃，起爆前，粒子沿板條炸藥的幾何形狀分布。所取時間步長為1.0×10-9s，起爆后產生一個平面爆炸波，取爆速為6 930 m/s，整個爆轟過程在14.4 μs左右完成。

圖4給出了一維板條TNT炸藥爆轟過程中的壓力、能量、密度和速度的數值計算結果，爆轟產物在C-J點處的壓力、密度、速度理論近似值分別為19.57 GPa，2 173 kg/m3，1 733 m/s1，對一維板條TNT爆轟問題，已有學者用實驗方法求得壓力的C-J峰值為21 GPa。從數值結果可知，SPH的求解值比理論方法更接近這個值，可見SPH方法的求解結果是合理的，且應用SPH方法能較好地預測出爆轟波的大小和形狀,以及在爆轟過程中的壓力分布，因此SPH方法非常適合解決炸藥爆轟問題。

圖4 一維板條TNT炸藥爆轟過程中的壓力、能量、密度和速度分布曲線

5 結 論

應用具有人工粘性的SPH方法能較好地預測出爆轟波的大小和形狀,以及在爆轟過程中的壓力分布，求解結果已達到了較高的精度,SPH方法非常適宜處理炸藥水下爆炸的極短瞬時具有大變形和高度非均勻的動力學極端情形。

[1] 光滑粒子流體動力學——一種無網格粒子法[M].韓旭,楊剛,強洪夫,等譯.長沙：湖南大學出版社，2005，10.

[2] 張鎖春.光滑質點流體動力學(SPH)方法(綜述)[J].計算物理, 1996,13 (4):385-397.

[3] 張守中.爆炸與沖擊動力學[M].北京:兵器工業出版社，1993.

[4] 陳朗,龍新平,馮長根,蔣小華.含鋁炸藥爆轟[M].北京:國防工業出版社, 2004.

[5] LIU M B, LIU G R, ZONG Z, LAM K Y. Computer simulation of high explosive explosion using smoothed particle hydrodynamics methodology [J]. Computers & Fluids,2003,32:305-322.

[6] LIU M B, LIU G R, LAM K Y, ZONG Z. Smoothed particle hydrodynamics for numerical simulation of underwater explosion [J]. Computational Mechanics,2003,30:106-118.

[7] LIU M B, LIU G R, LAM K Y. A one dimensional meshfree particle formulation for simulating shock waves[J].Shock Waves, 2003,13(3):201-211.