水下爆炸近場峰值壓力的數(shù)值模擬*

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

（1.中國科學(xué)技術(shù)大學(xué)近代力學(xué)系，安徽 合肥 230027；2.北京應(yīng)用物理與計(jì)算數(shù)學(xué)研究所，北京 100094）

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

本文中利用SPH方法進(jìn)行水下爆炸近場壓力的數(shù)值模擬，對模擬結(jié)果進(jìn)行基于波系傳播規(guī)律的定性分析，總結(jié)近場內(nèi)峰值壓力的變化規(guī)律，以期提高對水下爆炸近場壓力變化規(guī)律的認(rèn)識。

1 SPH方法

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

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

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

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

在軸對稱的柱坐標(biāo)系統(tǒng)中，歐拉方程組的相應(yīng)形式[11]為：

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

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

2 CJ爆轟模型

在計(jì)算炸藥爆轟過程時(shí)，使用了經(jīng)典的CJ爆轟理論模型[12]。假定炸藥微元在通過爆轟波陣面時(shí)瞬間轉(zhuǎn)變?yōu)楸óa(chǎn)物并完成放能，在這種假定下，反應(yīng)區(qū)的寬度變?yōu)榱悖蔀橐员賯鞑サ膹?qiáng)間斷面。假設(shè)波陣面的傳播速度為D，產(chǎn)物的粒子速度為u1，在自持定常爆轟狀態(tài)下，產(chǎn)物離開波陣面的速度D-u1等于當(dāng)?shù)芈曀賑1：

通常認(rèn)為爆速D保持不變，需要由實(shí)驗(yàn)測定。

3 狀態(tài)方程

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

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

爆炸產(chǎn)物的狀態(tài)方程選取JWL形式：

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

4 自由場爆炸驗(yàn)證

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

Cole[1]在大量理論研究與試驗(yàn)數(shù)據(jù)的基礎(chǔ)上，提出了自由場中水下爆炸氣泡脈動第一個(gè)周期和最大半徑的經(jīng)驗(yàn)公式：

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

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

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

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

Cole[1]總結(jié)了經(jīng)驗(yàn)峰值壓力公式，Zamyshlyaev等[14]后來對Cole的結(jié)果進(jìn)行了深入和擴(kuò)充，給出的峰值壓力公式為：

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

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

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

從對自由場水下爆炸的計(jì)算結(jié)果看，本文使用的SPH方法得到的計(jì)算結(jié)果較好地與經(jīng)驗(yàn)公式相符合，使用SPH方法進(jìn)行水下爆炸近場內(nèi)的壓力研究是可行的。

圖2 峰值壓力計(jì)算結(jié)果與經(jīng)驗(yàn)公式的比較Fig.2 Comparison of peak pressure between calculated results and empirical formula

5 水下爆炸近場中的波系

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

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

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

6 峰值壓力公式

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

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

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

平面二維：

平面一維：

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

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

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

測定系數(shù)

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

式中：峰值壓力單位為GPa，擬合結(jié)果與文獻(xiàn)結(jié)果[4]是相近的。

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

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

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

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

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

圖7 近場中l(wèi)nPm-ln(R/a)曲線的一階和二階導(dǎo)數(shù)Fig.7 First- and second-order derivatives of lnPm-ln(R/a) curve in near-field

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

7 結(jié) 論

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

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

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

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