不同真空度下空中爆炸近場(chǎng)特性的數(shù)值模擬研究

李科斌， 李曉杰,2， 閆鴻浩， 王小紅， 楊晨琛

(1.大連理工大學(xué) 工程力學(xué)系，遼寧 大連 116024；2.工業(yè)裝備結(jié)構(gòu)分析國(guó)家重點(diǎn)實(shí)驗(yàn)室，遼寧 大連 116024)

空中爆炸是爆炸力學(xué)研究中的一個(gè)基本現(xiàn)象，深入分析爆炸空氣沖擊波的傳播規(guī)律，特別是凝聚炸藥空中爆炸時(shí)的近場(chǎng)特性，對(duì)于研究武器試驗(yàn)、爆破工程、炸藥工業(yè)等領(lǐng)域的安全防護(hù)問題有重要意義，如露天爆破、爆炸焊接時(shí)沖擊波等有害效應(yīng)的控制，爆炸容器等結(jié)構(gòu)的抗爆防護(hù)設(shè)計(jì)等。

眾所周知，空中爆炸時(shí)，爆轟產(chǎn)物高速膨脹強(qiáng)烈壓縮周圍的空氣形成空氣沖擊波，而表征空中爆炸近場(chǎng)特性的特征參量包括峰值超壓Δpm、比沖量i、正壓作用時(shí)間t+、沖擊波到達(dá)時(shí)間t′等(如圖1所示)，為計(jì)算這些參量，學(xué)者們通過大量試驗(yàn)提出了各自的計(jì)算公式[1-8]，這些經(jīng)驗(yàn)公式具有不同的適用范圍，相互之間也存在一定的差異。

另一方面，當(dāng)外界空氣環(huán)境(如壓力、溫度、濕度、風(fēng)力等)改變時(shí)，上述特征參量的規(guī)律將會(huì)隨之發(fā)生變化[9]，即滿足薩克斯比例定律[10]，Dewey等[11]通過調(diào)節(jié)高度模擬器中的環(huán)境壓力和溫度，測(cè)定了多種尺寸裸彭托利特炸藥的爆炸波參數(shù)，得出的結(jié)論證實(shí)了薩克斯定律，該定律同樣由Ericsson等[12-13]的實(shí)驗(yàn)所證實(shí)；Veldman等[14]研究了環(huán)境壓力對(duì)反射沖擊波超壓和沖量的影響；Silnikov等[15]對(duì)比分析了高空中和常壓環(huán)境下爆炸沖擊波的參數(shù)變化及對(duì)飛行器的破壞影響。另一方面，Jack等[16]在炸藥近旁和模擬高空情況時(shí)的試驗(yàn)指出，上述條件下爆炸波的性質(zhì)將會(huì)發(fā)生改變，薩克斯定律此時(shí)不能用于特征參量的換算，如何準(zhǔn)確而又完整地描述爆炸近場(chǎng)特征參量的變化規(guī)律是本文將要解決的主要問題之一。

圖1 典型的空中爆炸p-t曲線Fig.1 A typical pressure-time curve for air explosion

AUTODYN有限元程序可很好地用于模擬爆炸沖擊波的傳播過程[17-20]，本文將基于該平臺(tái)建立空中爆炸的球?qū)ΨQ楔形模型，在與標(biāo)準(zhǔn)狀態(tài)經(jīng)驗(yàn)公式對(duì)比的基礎(chǔ)上，引入特征比例距離分析不同真空度下空中爆炸近場(chǎng)的特征參量變化規(guī)律，結(jié)論可為高空爆炸時(shí)目標(biāo)的毀傷效應(yīng)以及爆破工程的沖擊波防護(hù)提供參考。

1 計(jì)算模型與狀態(tài)方程

1.1 計(jì)算模型

本文采用AUTODYN程序中的球?qū)ΨQ楔形模型來模擬自由場(chǎng)中球形裝藥爆炸沖擊波的傳播，建立如圖2所示的計(jì)算模型。模型采用多物質(zhì)Euler算法，選用與實(shí)際對(duì)應(yīng)的ANFO炸藥，炸藥半徑與網(wǎng)格尺寸之比λ=100，根據(jù)文獻(xiàn)[18]和[19]的分析結(jié)果，所用計(jì)算模型和網(wǎng)格尺寸足以保證計(jì)算精度，從炸藥邊界開始設(shè)置一系列監(jiān)測(cè)點(diǎn)(編號(hào)1#～96#)用于特征參量的結(jié)果輸出。

圖2 自由場(chǎng)空中爆炸楔形計(jì)算模型Fig.2 Wedge model of free-space air explosion

1.2 參數(shù)設(shè)置

1.2.1 炸藥

采用AUTODYN標(biāo)準(zhǔn)材料模型庫里的ANFO炸藥，其狀態(tài)方程為經(jīng)典的JWL狀態(tài)方程，具體為

(1)

式中：P為爆壓，V為爆轟產(chǎn)物相對(duì)比容(爆轟產(chǎn)物體積與炸藥初始體積之比)，E0為炸藥單位體積的初始內(nèi)能，A、B、R1、R2、ω為炸藥的特征參數(shù)(常數(shù))，本文ANFO炸藥各參數(shù)取值如表1所示。

1.2.2 空氣

空氣的狀態(tài)方程近似采用理想氣體狀態(tài)方程，可描述為

表1 ANFO炸藥參數(shù)Tab.1 Parameters of ANFO

(2)

式中：P為空氣壓力，γ為理想氣體絕熱指數(shù)，取1.4；ρ為空氣密度；ρ0為標(biāo)準(zhǔn)狀態(tài)下的空氣密度，取1.225×10-3g/cm3；E0為空氣的初始比內(nèi)能，取2.068×105J/m3。本文通過改變空氣的初始密度近似實(shí)現(xiàn)不同的初始?jí)毫Γ磳?shí)現(xiàn)不同真空度的空氣環(huán)境。

2 爆炸相似律與經(jīng)驗(yàn)公式

為檢驗(yàn)?zāi)M結(jié)果的可靠性，通過對(duì)比經(jīng)驗(yàn)公式計(jì)算的結(jié)果分析數(shù)值計(jì)算的精度。根據(jù)薩克斯比例定律可知，包含超壓Δpm、比沖量i、正壓時(shí)間t+和某些空氣參數(shù)的無量綱群

特定地的是(Rp01/3/E1/3)的單一函數(shù)。其中，R為爆心距離，p0為外界空氣初壓，c0是環(huán)境空氣的聲速，E是炸藥的總能量(通常用裝藥質(zhì)量W代替)。當(dāng)炸藥種類、裝藥密度和空氣狀態(tài)不變時(shí)，即為一般形式的爆炸相似定律，也稱霍普金森比例定律，即

基于以上規(guī)律，學(xué)者們通過大量試驗(yàn)各自得到了計(jì)算沖擊波峰值超壓Δpm、正壓作用時(shí)間t+和比沖量i的經(jīng)驗(yàn)公式，但這些經(jīng)驗(yàn)公式的適用范圍和計(jì)算結(jié)果互有差異，以爆炸波峰值超壓為例，圖3給出了不同沖擊波峰壓經(jīng)驗(yàn)公式計(jì)算得到的結(jié)果與數(shù)值模擬結(jié)果之間的關(guān)系(R0為裝藥半徑)。

圖3 不同沖擊波峰壓經(jīng)驗(yàn)公式與模擬結(jié)果的關(guān)系

Fig.3 Comparison of simulated peak pressure with calculated results by empirical formulas

從圖3中可以看出：

① 各經(jīng)驗(yàn)公式在R>10～12R0后的計(jì)算結(jié)果趨于一致，因此經(jīng)驗(yàn)計(jì)算不適用于近場(chǎng)，其適用范圍一般為R>10～12R0；

② 數(shù)值計(jì)算結(jié)果與J. Henrych公式最為吻合，在R>6R0時(shí)，誤差在10%以內(nèi)，在R>12R0時(shí)，誤差小于5%；

③ 在裝藥近旁(R<1.6R0)，峰壓與距離之間存在其他對(duì)應(yīng)關(guān)系。

圖4和圖5則分別給出了利用不同經(jīng)驗(yàn)公式計(jì)算得到的正壓作用時(shí)間和比沖量。

圖4 不同經(jīng)驗(yàn)公式的正壓作用時(shí)間計(jì)算結(jié)果Fig.4 Calculated results of positive pressure duration by empirical formulas

可以看出，不同經(jīng)驗(yàn)公式關(guān)于正壓作用時(shí)間和比沖量的計(jì)算結(jié)果存在很大差別，無法進(jìn)行參照對(duì)比，因此這兩個(gè)參量的模擬結(jié)果不適合與經(jīng)驗(yàn)公式作比較，文獻(xiàn)[21]中通過數(shù)值計(jì)算得到了如圖6和圖7所示的不同密度PETN炸藥的比沖量和正壓作用時(shí)間變化曲線，實(shí)線1代表密度1.6 g/cm3的PETN，虛線2密度為0.4 g/cm3，此處Re=(E/pM)1/3，tM=(E/pM)1/3/(pM/ρM)1/2，iM=pMtM。圖中同時(shí)給出了本文常壓下ANFO炸藥的模擬結(jié)果，雖然與文獻(xiàn)選用的炸藥種類不同，但兩者在變化規(guī)律上基本一致。

圖5 不同經(jīng)驗(yàn)公式的比沖量計(jì)算結(jié)果Fig.5 Calculated results of specific impulse by empirical formulas

圖6 比沖量的數(shù)值計(jì)算結(jié)果對(duì)比Fig.6 Comparison of simulated specific impulse curves

圖7 正壓作用時(shí)間的數(shù)值計(jì)算結(jié)果對(duì)比Fig.7 Comparison of simulated positive pressure duration curves

通過以上分析和對(duì)比可知，本文采用的數(shù)值方法是可行的，利用AUTODYN有限元程序計(jì)算空中爆炸的結(jié)果具有較高置信度。

3 不同真空度下的模擬結(jié)果與分析

3.1 真空度對(duì)峰值超壓的影響

圖8 不同真空度下凝聚炸藥爆炸的峰值超壓曲線

Fig.8 Peak overpressure curves of condensed explosive at different vacuum degrees

通過觀察圖8各曲線發(fā)現(xiàn)，爆炸峰壓的變化呈現(xiàn)明顯的區(qū)間性，裝藥近旁和遠(yuǎn)處具有不同的變化規(guī)律，峰壓首先沿p0=0 atm(即完全真空)曲線變化，在某一距離后，爆炸峰壓的衰減即為空氣沖擊波的衰減，這是因?yàn)檠b藥近區(qū)的特征參量受爆轟產(chǎn)物的影響。圖9顯示了1標(biāo)準(zhǔn)大氣壓下距爆炸中心不同距離處壓力隨時(shí)間變化的曲線群，在沖擊波傳播的初始階段(4#Gauss點(diǎn))，壓力出現(xiàn)了兩個(gè)峰值，分別對(duì)應(yīng)爆轟產(chǎn)物-空氣界面前端的沖擊波峰值和界面后端的產(chǎn)物壓力峰值，在最初階段，空氣沖擊波峰壓小于爆轟產(chǎn)物峰壓(圖8中臨界點(diǎn)前的初始段)，隨著爆轟產(chǎn)物的急劇膨脹，沖擊波峰壓迅速超過產(chǎn)物峰壓，因此圖8中各真空度峰壓曲線的大部分體現(xiàn)的是空氣沖擊波峰壓的衰減過程。這樣，通過對(duì)爆轟產(chǎn)物自由膨脹曲線(p0=0 atm)和空氣沖擊波峰壓衰減曲線分別擬合，就可以得到爆炸場(chǎng)全范圍內(nèi)峰值超壓的表達(dá)式。

圖9 不同距離處壓力-時(shí)間變化曲線Fig.9 Pressure-time curve at different distances

從圖9中還可以看出，爆轟產(chǎn)物對(duì)壓力的影響逐漸減小，圖中56#監(jiān)測(cè)點(diǎn)距離爆炸中心12R0，該處的壓力衰減基本趨于平滑，此處也正是爆轟產(chǎn)物極限膨脹體積的位置，而大部分經(jīng)驗(yàn)公式的適用范圍也正是在極限膨脹體積之后。

為分析外界初始空氣狀態(tài)對(duì)爆炸特征參量的影響，我們定義特征比例距離

(3)

式中:fd=(ρ0/ρM)1/3=[(p0/pM)/(T0/TM)]1/3，稱為初始狀態(tài)比例系數(shù)，下標(biāo)0表示初始空氣狀態(tài)，M表示參考狀態(tài)(pM=0.101 325 MPa，TM=273 K)；R為爆心距離，m；W為炸藥TNT當(dāng)量，kg。

圖10給出了以特征比例距離Z表示的不同空氣環(huán)境的峰壓變化規(guī)律，計(jì)算結(jié)果與相關(guān)實(shí)驗(yàn)結(jié)果[22]符合良好。需要說明的是，由于本文旨在研究初始?jí)毫Φ挠绊?，因此?jì)算結(jié)果中空氣的初始溫度都取定為288 K。

圖10 以特征比例距離表示的沖擊波峰壓曲線

Fig.10 Peak overpressure of blast wave with the scaled distanceZas variables

考察圖10中本文計(jì)算結(jié)果可以發(fā)現(xiàn)，在該坐標(biāo)系中，對(duì)于同一炸藥，不同初始空氣壓力的沖擊波峰壓衰減趨于同一變化規(guī)律，擬合后的方程形式為

(4)

式中:Δpm為峰值超壓，MPa；p0為空氣域初始?jí)毫?，MPa。對(duì)于本文所選炸藥，α=5.565，β=1.078。

分析圖8中p0=0 atm的曲線，爆轟產(chǎn)物自由擴(kuò)散時(shí)峰壓的變化規(guī)律可擬合為

(5)

于是，我們就可以得到不同初始空氣壓力下全范圍峰壓的分段函數(shù)

(6)

其中，α、β、λ和ξ由炸藥的性質(zhì)所決定，對(duì)于AUTODYN程序提供的ANFO炸藥，上述系數(shù)分別為5.565、1.078、3.902和5.970。Rs為不同初始?jí)毫?duì)應(yīng)的臨界距離，它與真空度的關(guān)系如圖11所示。

圖11 不同真空度的峰壓臨界距離Fig.11 Critical distance vs vacuum degree

擬合后得

(7)

3.2 比沖量的變化規(guī)律

圖12 不同真空度下的比沖量變化Fig.12 Specific impulse curves of various vacuum degrees

為具體分析比沖量極值點(diǎn)的位置，我們以1標(biāo)準(zhǔn)大氣壓時(shí)的情況為例，給出了八個(gè)相繼時(shí)刻下20R0范圍內(nèi)的相對(duì)壓力分布(如圖13所示)。

圖13 1 atm下不同時(shí)刻的壓力分布曲線Fig.13 Pressure distributions at different times

首先，我們對(duì)t=0.085 ms時(shí)刻(分布4)進(jìn)行單獨(dú)考察。圖14給出了該時(shí)刻下的壓力分布曲線及云圖，圖中產(chǎn)物-空氣界面前是空氣沖擊波陣面，而在爆轟產(chǎn)物內(nèi)部，從界面反射回的稀疏波又從爆心產(chǎn)生了又一次反射，稀疏波尾部的壓力達(dá)到最小值，而主沖擊波后的壓力還沒能正常衰減到稀疏波尾部的狀態(tài)，因此在稀疏波尾部就形成了如圖所示的間斷面，該間斷面最初是在運(yùn)用特征線法進(jìn)行爆炸數(shù)值計(jì)算中發(fā)現(xiàn)的[23]，一些文獻(xiàn)中闡述此間斷面時(shí)稱其為“二次沖擊”(secondary shock)[24-25]。

圖14 不同真空度下的比沖量變化Fig.14 Specific impulse curves of various vacuum degrees

繼續(xù)觀察圖13，間斷面在初始階段由于爆轟產(chǎn)物的推動(dòng)作用由內(nèi)向外移動(dòng)，直到分布7(間斷面到達(dá)約10R0處)之后才開始向爆心運(yùn)動(dòng)。圖12中，曲線波谷對(duì)應(yīng)圖13中的分布4，該時(shí)刻稀疏波尾部的壓力等于初壓p0；而曲線的波峰對(duì)應(yīng)圖13中的分布7，該時(shí)刻則是間斷面?zhèn)鞑シ较虻霓D(zhuǎn)折點(diǎn)，即由向外移動(dòng)開始向爆心運(yùn)動(dòng)。其他真空度下的比沖量規(guī)律類似，這種非單調(diào)性變化都反映了爆轟產(chǎn)物中稀疏波尾部間斷面形成和傳播的波動(dòng)特性。

3.3 正壓作用時(shí)間

不同真空度下正壓作用時(shí)間的變化規(guī)律如圖15所示，可以發(fā)現(xiàn)，正壓作用時(shí)間隨比例距離Z呈分段變化規(guī)律，各折點(diǎn)位置同樣與稀疏波尾部的間斷面相關(guān)，其中A點(diǎn)對(duì)應(yīng)圖13中的分布4，即此時(shí)稀疏波尾部壓力等于p0；在此之前，真空度小的曲線始終位于真空度大的曲線上方，此后，各曲線重疊在一起，即可用特征比例距離Z為變量的同一方程來描述不同真空度下的正壓作用時(shí)間變化規(guī)律。折點(diǎn)B則對(duì)應(yīng)圖13中的分布7，即此時(shí)間斷面開始向爆心運(yùn)動(dòng)；折點(diǎn)C對(duì)應(yīng)爆轟產(chǎn)物膨脹的極限體積處，而在BC之間的范圍內(nèi)，正壓作用時(shí)間與特征比例距離呈線性關(guān)系。

3.4 產(chǎn)物界面和波陣面的運(yùn)動(dòng)

通過程序給定的監(jiān)測(cè)點(diǎn)的超壓數(shù)據(jù)可以得到峰壓的時(shí)程曲線，即為空氣沖擊波波陣面的運(yùn)動(dòng)軌跡，如圖16所示；而通過追蹤爆轟產(chǎn)物界面的運(yùn)動(dòng)軌跡，則可以得到產(chǎn)物界面運(yùn)動(dòng)的時(shí)程曲線，如圖17所示，兩者經(jīng)特征化后的曲線也相應(yīng)地給出，圖中tM含義與前述相同(能量E用質(zhì)量W代替)。

圖16 空氣沖擊波波陣面時(shí)程曲線Fig.16 Time-history curves of shock wave front

分析圖16和圖17可知，在常規(guī)坐標(biāo)系中，真空度越大，沖擊波陣面和產(chǎn)物界面運(yùn)動(dòng)速度越快，且在一定距離后，沖擊波速度各自趨于一恒定值，產(chǎn)物界面速度則趨于0，即爆轟產(chǎn)物膨脹至極限體積，而完全真空狀態(tài)時(shí)，爆轟產(chǎn)物膨脹速度為一恒定值(U∞≈6 600 m/s)。不同真空度下產(chǎn)物界面和沖擊波陣面時(shí)程曲線在相應(yīng)的特征坐標(biāo)系中都各自歸結(jié)為一條運(yùn)動(dòng)軌跡，即在以Z表示的系統(tǒng)中，盡管外界的真空環(huán)境不同，但產(chǎn)物界面和沖擊波陣面的運(yùn)動(dòng)規(guī)律都是相同的。

圖17 爆轟產(chǎn)物界面時(shí)程曲線Fig.17 Time-history curves of detonation product interface

不同真空度下爆轟產(chǎn)物的極限膨脹體積與初始空氣壓力的關(guān)系如圖18所示，可以推算，它的垂直漸近線為Y軸，水平漸近線為X軸；而極限體積與臨界距離之間呈線性關(guān)系(圖19)，即

(8)

圖18 極限膨脹體積與初始?jí)毫Φ年P(guān)系Fig.18 Ultimate expansion volume vs. initial pressure

圖19 極限膨脹體積與臨界距離的關(guān)系Fig.19 Ultimate expansion volume as a function of critical distance

4 結(jié) 論

本文基于有限元程序AUTODYN-2D，建立了空中爆炸的一維球?qū)ΨQ楔形模型，通過改變空氣域的初始參數(shù)，得到了不同真空度下爆炸近場(chǎng)特征參量的變化規(guī)律，主要結(jié)論如下：

(1)常壓下的模擬結(jié)果與峰壓經(jīng)驗(yàn)公式符合良好，尤其與J. Henrych最為吻合，在R>12R0時(shí)，誤差小于5%，因此數(shù)值計(jì)算具有較高的置信度。

(3)不同真空度下比沖量i和正壓作用時(shí)間t+的變化具有非單調(diào)性，它們的極值點(diǎn)和折點(diǎn)與稀疏波尾部間斷面密切相關(guān)，當(dāng)稀疏波尾部壓力等于p0時(shí)，i和t+達(dá)到極小值，而當(dāng)間斷面的傳播方向由外向爆心轉(zhuǎn)變時(shí)，i達(dá)到極大值，t+則出現(xiàn)折點(diǎn)。不同真空度在稀疏波尾部壓力等于p0之后具有相同的比沖量和正壓作用時(shí)間；此外，不同真空度下的爆轟產(chǎn)物界面和沖擊波陣面時(shí)程曲線各自重合，極限體積與臨界距離之間呈線性關(guān)系。