1 引 言

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

2 爆轟理論基礎(chǔ)

2.1 炸藥爆轟原理

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

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

圖2 Hugoniot曲線和Rayleigh曲線

2.2 爆轟控制方程

由于爆炸和膨脹的速度相當(dāng)高，故可假設(shè)氣體生成物是非粘性的且爆炸過程是絕熱的。因此可利用如下的歐拉方程與恰當(dāng)?shù)臓顟B(tài)方程來模擬爆炸過程。

(1)

2.3 爆轟產(chǎn)物狀態(tài)方程

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

JWL狀態(tài)方程于1965年由美國勞倫斯利弗莫爾國家實(shí)驗(yàn)室的Lee E L在Jones和Wilkins的工作基礎(chǔ)上提出的。其形式為:

(2)

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

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

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

3.1 人工粘性

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

(3)

3.2 光滑長度

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

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

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

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

(4)

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

3.3 近鄰粒子搜索

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

圖3 直接搜索法

3.4 SPH離散公式

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

(5)

3.5 時間積分

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

(6)

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

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

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

(7)

4 SPH爆轟數(shù)值的模擬驗(yàn)證

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

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

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

5 結(jié) 論

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

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

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

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

[4] 陳朗,龍新平,馮長根,蔣小華.含鋁炸藥爆轟[M].北京:國防工業(yè)出版社, 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.