水下爆炸近場峰值壓力的數值模擬*

宮翔飛，劉文韜，張樹道，楊基明

（1.中國科學技術大學近代力學系，安徽 合肥 230027；2.北京應用物理與計算數學研究所，北京 100094）

壓力狀態研究是水下爆炸研究中的重要內容。以往水下爆炸壓力狀態研究的距離范圍主要集中于中遠場，現有的中遠場水下爆炸研究結果與實驗符合較好。水下爆炸近場壓力則具有典型的高壓高頻特點，壓力峰值至少在GPa量級，持續時間則為ms量級，給實驗測量造成很大困難[1]。近期，學者們對近場峰值壓力進行了試驗和數值研究[2-8]，認為在R/a≥1的一定范圍內峰值壓力的變化更符合指數規律，而在1≤R/a≤6范圍進行擬合時則都使用了冪次形式，最后的擬合公式也不盡相同。這些試驗和計算工作都證實，水下爆炸沖擊波中存在相似原理，因此用數值模擬對水下爆炸近場中不同波系的影響進行一些規律性的分析是可行的。

本文中利用SPH方法進行水下爆炸近場壓力的數值模擬，對模擬結果進行基于波系傳播規律的定性分析，總結近場內峰值壓力的變化規律，以期提高對水下爆炸近場壓力變化規律的認識。

1 SPH方法

SPH方法[9-11]對物理量進行局部空間的積分近似和粒子近似，是一種拉氏無網格粒子類方法，在水下爆炸領域得到了越來越廣泛的應用。本文計算中忽略黏性、熱傳導，使用的SPH控制方程為歐拉方程組：

式中：ρ、p、e分別為粒子密度、靜水壓和單位質量的內能，m為 粒子質量，v為粒子速度，下標i、j為粒子編號，vij=vi-vj，N為粒子總數目，x表示坐標軸，α ,β=1,2,3，xα或xβ分別表示不同的坐標軸方向,光滑函數Wij(R,h)有兩個參數，h為光滑長度，R為i、j粒子間的距離與h的比值。

光滑函數W取為最常用的三次樣條函數形式：

式中：bd在不同維數中取不同的值，一維b1=1/h，二維b2=15/(7πh2)，三維b3=3/(2πh3)。

在軸對稱的柱坐標系統中，歐拉方程組的相應形式[11]為：

式中：s(r,z)為粒子在柱坐標中的位置，。

計算中沒有考慮重力，為了模擬水深，在水的初始狀態中增加了初始壓力，比如若水深H=0，則取初始壓力p0=1.01×105Pa。

2 CJ爆轟模型

在計算炸藥爆轟過程時，使用了經典的CJ爆轟理論模型[12]。假定炸藥微元在通過爆轟波陣面時瞬間轉變為爆炸產物并完成放能，在這種假定下，反應區的寬度變為零，成為以爆速傳播的強間斷面。假設波陣面的傳播速度為D，產物的粒子速度為u1，在自持定常爆轟狀態下，產物離開波陣面的速度D-u1等于當地聲速c1：

通常認為爆速D保持不變，需要由實驗測定。

3 狀態方程

眾多學者對水下爆炸中適用的水的狀態方程進行了數值模擬與試驗的對比研究，認為在PH＜32.0 GPa，ρ＜2.02 g/cm3范圍內[13]，多項式形式的狀態方程能夠適用。本文中計算的TNT藥包水下爆炸算例中，水中的壓力峰值也在此范圍內，因此水的狀態方程選用了多項式形式：

式中：ρ0為初始密度，μ為壓縮度， μ= ρ/ρ0-1，EM是單位質量內能增量，A1、A2、A3、T1、T2為壓強量綱擬合系數，B0、B1為無量綱系數。其具體參數如下：

爆炸產物的狀態方程選取JWL形式：

式中：E為比內能，ρ0為炸藥初始密度，V=ρ0/ρ，A、B、R1、R2、ω為常系數。對TNT炸藥，常用值為A=373.8 GPa、B=3.747 GPa、R1=4.15、R2=0.9、ω=0.3。炸藥起爆方式使用時間起爆，爆速取D=6 930 m/s。

4 自由場爆炸驗證

先將自由場水下爆炸中氣泡脈動周期、最大半徑和中場峰值壓力的計算結果與經驗公式相比較，驗證SPH計算中相似率的符合，結果表明了本文數值計算的可靠性。

Cole[1]在大量理論研究與試驗數據的基礎上，提出了自由場中水下爆炸氣泡脈動第一個周期和最大半徑的經驗公式：

式中：T為氣泡第一次脈動周期，s；Rm為氣泡最大半徑，m；W為炸藥質量，kg；H為爆炸深度，m。對TNT炸藥，一般取系數值k1=2.11、k2=3.8。為了加快計算速度，取水深H=40 m，根據經驗公式，初始半徑a=0.1 m(W=6.83 kg)的球形TNT炸藥，氣泡第一次振動的周期最大半徑1.94 m。

使用二維軸對稱SPH算法，水中初始壓力為50.33 m高水柱在重力作用下產生的壓力值，進行自由場球形藥包中心起爆水下爆炸氣泡的計算。得到水下爆炸氣泡第一次振動的周期和最大半徑分別為T=0.159 s、Rm=1.98 m。計算結果與經驗公式相比的誤差分別為4.6%和2.1%，是較為接近的。

Cole[1]還進行了水下爆炸中相似率的理論論證，并用實驗數據作了簡單而充分的證明。根據這一理論，以(R/a,t/a)作為參數，不同裝藥量的同種炸藥具有相同的壓力分布。

本文中對半徑為0.1 m和0.2 m(W=54.62 kg)的球形TNT炸藥進行計算，取相同比距離R/a=10處的壓力曲線，a=0.2 m的模型中，炸藥半徑為另一模型的2倍，根據相似率其時間尺度也是另一模型的2倍。從圖1中可見，以t/a作為參數時，兩種計算模型的壓力衰減曲線基本吻合，不同計算模型中，相同R/a處的峰值壓力也是基本一致的，說明相似率符合很好。

Cole[1]總結了經驗峰值壓力公式，Zamyshlyaev等[14]后來對Cole的結果進行了深入和擴充，給出的峰值壓力公式為：

圖1 不同模型壓力的相似率Fig.1 Similarity of pressure in various models

式中：a為球形炸藥的初始半徑，R為測點到炸藥中心的距離。

本文中采用SPH二維軸對稱格式，取H=0，對球形裝藥進行爆轟模擬。比較6≤R/a≤12內的計算結果與經驗公式，令ρ=1 630 kg/m3，將經驗公式中球形藥包的質量與半徑進行了換算，如圖2所示，兩者的結果最大差別不超過2%。

從對自由場水下爆炸的計算結果看，本文使用的SPH方法得到的計算結果較好地與經驗公式相符合，使用SPH方法進行水下爆炸近場內的壓力研究是可行的。

圖2 峰值壓力計算結果與經驗公式的比較Fig.2 Comparison of peak pressure between calculated results and empirical formula

5 水下爆炸近場中的波系

爆炸產物內的波系對水下爆炸近場中的壓力有重要的影響，不同波的傳播速度和衰減規律并不相同，以比距離(R/a)作為參數，同一波形在不同范圍區間上的規律也是有差異的，多種復雜的衰減規律綜合在一起就造成了水下爆炸近場壓力狀態的復雜性。

水下爆炸近場內的波系結構主要決定于炸藥內傳播的爆轟波，使用CJ爆轟模型計算半徑0.1 m (W=6.83 kg)的TNT球形炸藥在中心起爆的波系，爆轟波由前導沖擊波和之后的稀疏波區組成，由計算結果結合波的傳播理論，得到圖3中的波系結構。需要注意的是，圖中以實線表示的激波、爆炸產物內的稀疏波區在向左的中心稀疏波到達之前，以及反射稀疏波波頭IU線為計算所得，其他虛線表示的稀疏波區則為示意圖。也就是說，圖3中I、M、U、S的位置為計算所得，爆炸產物和水界面上的K、L點的確切位置則是未知的，而J點的實際橫坐標R/a＜1.1，為了從波系圖中清晰描述反射壓縮波的匯聚過程，才標示在圖中位置。沖擊波OI在I點(爆炸產物和水的界面)形成向水內的透射沖擊波IC和向爆炸產物傳播的反射中心稀疏波區IUS，能夠匯聚到O點的左傳稀疏波會形成反射稀疏波向右傳播。爆炸產物內的稀疏波區OIK在界面IK段上都會形成向水內的透射稀疏波，而由于中心稀疏波區IUS的作用，界面左側爆炸產物內的聲阻抗迅速降低，由高于界面右側水的聲阻抗降低到低于水的聲阻抗，因此界面IK段會分為兩段，IJ段上反射波為壓縮波，JK段上反射波為稀疏波。IJ段上的反射壓縮波區在追趕稀疏波區IUS時，壓縮波的波頭會不斷進入稀疏波的波尾形成一條包絡線IS，同時壓縮波的波尾特征線會不斷匯聚到前方的壓縮波特征線上也形成一條包絡線JS，IS和JS的距離(壓縮波寬度)越來越小，形成二次沖擊波，二次沖擊波匯聚到O點后反射形成反射沖擊波SM。

圖3 水下爆炸近場內的波系Fig.3 Waves in near-field underwater explosion

6 峰值壓力公式

本文中取水深H=0，使用平面一維、平面二維和二維軸對稱格式的SPH方法計算半徑a=0.1 m的炸藥水下爆炸，平面一維、平面二維和二維軸對稱不同坐標系中水下爆炸的峰值壓力曲線如圖4所示。雖然維數不同，沖擊波傳入水中的初始壓力較為相近，都在15 GPa左右，峰值壓力隨距離增大而降低，峰值壓力在同等距離上降低的數值大體上隨距離增大而減小。由于幾何效應，峰值壓力衰減在高維中比低維更快。

對于6≤R/a≤12范圍的峰值壓力，本文中已經將二維軸對稱峰值壓力的計算結果與經驗公式進行了比較，符合冪次關系。由于維數效應，越高維中壓力衰減越快，各波形影響區域越窄，三維峰值壓力中以冪次變化的區域起點R/a=6，R/a跨度6～12，因此低維中相應波形影響的峰值壓力冪次變化區域起點R/a≥6，跨度不小于6。對R/a≥6的平面一維、平面二維峰值壓力計算結果進行統計擬合，可以發現，的確有這樣每點誤差都很小的區間：

圖4 峰值壓力曲線Fig.4 Peak pressure curves

平面二維：

平面一維：

將平面一維的擬合范圍10≤R/a≤30與圖4進行對比可知，R/a=10正對應于壓力峰值曲線中平臺的末端，也就是圖3中的L點，據此認為在平面二維、二維軸對稱中，冪次擬合范圍的起點R/a=6對應于L點也是合理的。由此可進行以下推論：在平面二維、二維軸對稱水下爆炸問題中，R/a＜6的近場區域內，存在1＜β＜6，在1＜R/a＜β的范圍內，壓力峰值受IJ段影響，降低速度要快于6＜R/a＜12區域，β＜R/a＜6范圍，則符合維數效應的冪次關系。

趙繼波等[6]的研究結果認為，在1＜R/a＜4.7(140 mm/30 mm)的近場內采用指數形式能夠較好地描述峰值壓力的規律，與本文中認為的在1＜R/a＜β的一定范圍內，峰值壓力的降低速度要高于冪次形式是相符合的。而在1＜R/a＜6的區域上采用統一形式進行分段擬合，還需要考慮到β＜R/a＜6(文獻[6]中β=4.7)范圍內冪次形式與指數形式的較大差別。

本文中二維軸對稱計算結果使用線性回歸方法擬合曲線，根據最小二乘法，得到lnPm-ln(R/a)曲線擬合直線y=r+sx的系數估計r=22.92，s=-2.08，同樣方法得到lnPm-R/a擬合直線的系數估計r=23.48，s=-0.84。在線性回歸理論中，用R2表示直線擬合曲線的符合程度，R2越接近于1表示曲線越近似為直線。

測定系數

式中：樣本平均值為樣本數。通過計算可得到lnPm-ln(R/a)的擬合直線R2=0.983，lnPm-R/a的擬合直線R2=0.901。lnPm-ln(R/a)的R2更接近于1，說明在R/a＜6的近場范圍內，采用冪次形式的擬合公式能夠比指數形式的擬合公式產生更小的偏差。其擬合直線為：

式中：峰值壓力單位為GPa，擬合結果與文獻結果[4]是相近的。

在圖5中lnPm-ln(R/a)計算曲線與擬合直線相比具有一定的下凸，這也意味著總體的衰減速度要快于冪次規律。對lnPm-ln(R/a)曲線進行冪次擬合能得到更加接近于計算曲線的擬合結果，圖6中新的曲線形式為：

圖5 近場中的lnPm-ln(R/a)和lnPm-R/a曲線Fig.5 lnPm-ln(R/a) curves and lnPm-R/a curves in near-field

前文已經指出，R/a＜6范圍內的擬合曲線需要分為兩段，而圖6中的新結果在ln(R/a)＜0.5范圍與數值結果曲線有較大差異。圖7中分別為lnPmln(R/a)曲線的一階和二階導數，從曲線中可以看到在ln(R/a)=0.45，也就是R/a=1.57附近為明顯的拐點，這與文獻[4]中選擇的R/a=1.52作為分界點是較為相近的。

圖8以R/a=1.5為界，進行了分段擬合，最后的擬合形式能夠使擬合曲線與計算曲線更加符合。對1＜R/a＜1.5范圍的擬合，一方面計算中這部分的誤差較大，計算結果只能作為參考，另一方面在擬合形式上也是為了與1.5＜R/a＜6內的擬合形式保持一致，因而沒有采用指數形式。

圖6 近場中lnPm-ln(R/a)的擬合曲線Fig.6 Fitted curve of lnPm-ln(R/a) in near-field

圖7 近場中lnPm-ln(R/a)曲線的一階和二階導數Fig.7 First- and second-order derivatives of lnPm-ln(R/a) curve in near-field

圖8 近場中lnPm-ln(R/a)的分段冪次擬合曲線Fig.8 Fitted curves of piecewise power of lnPm-ln(R/a) in near-field

7 結 論

用SPH方法進行了水下爆炸的數值模擬，通過自由場氣泡、相似率以及中場峰值壓力與試驗結果的符合，驗證了程序的有效性。進行了水下爆炸近場內波的傳播分析，通過波系與數值結果的對比，分析了水下爆炸近場內的壓力狀態并采用新的形式進行了公式擬合。得到以下結論：

(1)R/a=6是水下爆炸近場波形變化的一個分界點。

(2)R/a＜6區域的峰值壓力衰減速度分為明顯的兩段，R/a＜1.5的范圍內，水下爆炸峰值壓力降低的速度非常快。

(3) 使用冪次函數對曲線進行分段擬合，能夠在水下爆炸近場峰值壓力公式擬合中得到與數值結果更加符合的曲線。