周豐峻，鄭磊，孫云厚，張鏖，夏明

（1. 中國(guó)人民解放軍軍事科學(xué)院國(guó)防工程研究院，北京 100850；2. 中國(guó)人民解放軍92611部隊(duì)，廣州 510715）

一、TNT裝藥強(qiáng)爆炸自模擬解

點(diǎn)源爆炸理論是研究球形裝藥爆炸沖擊波最簡(jiǎn)單的模型之一，其假定爆轟產(chǎn)物質(zhì)量無(wú)限小，而裝藥爆炸釋放的能量有限，忽略運(yùn)動(dòng)中靜止空氣的影響，因此，不存在任何一個(gè)可以表征爆炸的初始線尺度參數(shù)，任何一個(gè)物理量在空間上的分布隨時(shí)間的變化是自相似的，即離爆心不遠(yuǎn)處沖擊波后的氣體運(yùn)動(dòng)滿足自相似條件 [1]。Taylor [2]、Von Neumann [3]、Bach和Lee [4]以及Sedov [5]各自獨(dú)立得到了點(diǎn)源強(qiáng)爆炸自模擬的理論解，與實(shí)驗(yàn)比對(duì)可以很好地描述爆炸初始階段的流場(chǎng)特征 [6]。本文基于點(diǎn)源爆炸理論，對(duì)平面、柱面以及球面的一維變幾何對(duì)稱(chēng)面（面源、線源和點(diǎn)源）裝藥強(qiáng)爆炸解的自模擬解流場(chǎng)進(jìn)行了分析。

（一）控制方程

爆炸沖擊波的傳播是一個(gè)非常復(fù)雜的不定常過(guò)程，但滿足質(zhì)量守恒、動(dòng)量守恒和能量守恒關(guān)系。由于爆炸力學(xué)所考慮的運(yùn)動(dòng)變化快，以至流線之間不會(huì)發(fā)生明顯的動(dòng)量和能量遷移，因此可以忽略粘性和熱傳導(dǎo) [7]。

1. 歐拉坐標(biāo)下的控制方程

歐拉坐標(biāo)下流場(chǎng)的狀態(tài)量主要包括壓力P、密度ρ、質(zhì)點(diǎn)流動(dòng)速度u、內(nèi)能E以及溫度T，這些量值均為自變量x, y , z , t的函數(shù)。

質(zhì)量守恒：

動(dòng)量守恒：

式（2）中，F(xiàn)為質(zhì)量力。

能量守恒：

物態(tài)方程：

對(duì)于一維對(duì)稱(chēng)TNT裝藥爆炸可將介質(zhì)的運(yùn)動(dòng)看作是一維絕熱運(yùn)動(dòng)，即介質(zhì)中各量值只依賴(lài)時(shí)間和某單一幾何坐標(biāo)而變化，去掉質(zhì)量力，即F = 0，可將其在平面直角坐標(biāo)系、柱坐標(biāo)系和球形坐標(biāo)系下建立的歐拉坐標(biāo)下的控制方程組統(tǒng)一寫(xiě)為：，u , v , w分別為

式（5）中，r為歐拉位置坐標(biāo)；t為時(shí)間；N為常數(shù)，等于0，1和2時(shí)分別對(duì)應(yīng)一維平面對(duì)稱(chēng)、一維柱面對(duì)稱(chēng)和一維球面對(duì)稱(chēng)情形。

2. 拉格朗日坐標(biāo)下的控制方程

為了提高沖擊波陣面局部計(jì)算精度，可將方程組寫(xiě)為拉格朗日形式，控制方程為：

再補(bǔ)充一個(gè)運(yùn)動(dòng)關(guān)系式：

式（6）中，R為拉格朗日位置坐標(biāo)；u為質(zhì)點(diǎn)速度；P為壓力；E為內(nèi)能；ρ為密度；ρ0為爆炸波前初始空氣密度；t為時(shí)間；N為相應(yīng)平面、柱面和球面坐標(biāo)指數(shù)。

為了計(jì)算方便，引入質(zhì)量坐標(biāo)方程式：

則可將拉格朗日坐標(biāo)方程組轉(zhuǎn)化為無(wú)量綱形式：

式（9）中，N為常數(shù)，等于0，1和2時(shí)分別為包含一維平面、柱面和球面的一體無(wú)量綱方程組。

（二）真實(shí)空氣條件下TNT裝藥爆轟產(chǎn)物狀態(tài)方程

1. 真實(shí)空氣無(wú)量綱形式狀態(tài)方程

試驗(yàn)表明，當(dāng)溫度大于等于2 000 K或沖擊波超壓大于等于4 MPa時(shí)，實(shí)際空氣不再滿足理想氣體狀態(tài)方程 [7]。隨著溫度與壓力的不斷升高，空氣中各種分子之間會(huì)發(fā)生化學(xué)反應(yīng)，如氮氧分子開(kāi)始離解，經(jīng)反應(yīng)生成NO和NO2?？諝饨M分在不同溫度下是不同的，當(dāng)溫度高于8 000 K后，空氣中各種原子會(huì)發(fā)生電離反應(yīng)，空氣組分變得更加復(fù)雜；當(dāng)溫度達(dá)到300 000 K后，任何密度下氣體的電離過(guò)程都可以完成。在溫度升高過(guò)程中發(fā)生的空氣中各原子電離、分子離解及激發(fā)的各種分子振動(dòng)均需要消耗能量，因此，在高溫條件下，空氣的狀態(tài)方程與理想氣體狀態(tài)方程有顯著不同 [8,9]。俄羅斯學(xué)者Kuznetzov在試驗(yàn)基礎(chǔ)上，采用統(tǒng)計(jì)物理學(xué)和沖擊波氣體動(dòng)力學(xué)的理論方法得到了空氣熱力學(xué)函數(shù)表和沖擊絕熱線的計(jì)算表格，Brode [10]對(duì)表格數(shù)據(jù)進(jìn)行擬合得到了真實(shí)空氣狀態(tài)方程的擬合公式 ：

式（10）～式（16）中，e，θ，π，η分別為無(wú)量綱內(nèi)能、無(wú)量綱溫度、無(wú)量綱壓力和無(wú)量綱密度。

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

Brode [10]從能量角度出發(fā)，對(duì)TNT裝藥爆轟產(chǎn)物的狀態(tài)方程進(jìn)行了詳細(xì)研究，并給出了與真實(shí)空氣狀態(tài)方程形式相同的描述方程：

式（19）中，α1= 1.76，b1= 52.4892，c1= 7.3。

（三）面源、線源和點(diǎn)源TNT裝藥爆轟階段流場(chǎng)參數(shù)計(jì)算

利用量綱分析方法中的Π定理，可得到由主定參量初始密度ρ0、大氣聲速C0、爆炸所釋放的能量E、球形坐標(biāo)r和時(shí)間t構(gòu)成的面源、線源和點(diǎn)源TNT裝藥強(qiáng)爆炸的獨(dú)立自模擬解變量分別為ξ =若沖擊波陣面的位置為rs，則ξ取值為ξ0，則可得到面源、線源和點(diǎn)源TNT對(duì)稱(chēng)裝藥強(qiáng)爆炸自模擬解的控制方程組分別如下：

面源裝藥：

線源裝藥：

點(diǎn)源裝藥：

邊界條件：在波陣面上，f(1) = g(1) = h(1) = 1；在中心處，f(0) = 0。其中，f，g和h分別為無(wú)量綱速度、無(wú)量綱密度和無(wú)量綱壓力。

式（21）～式（23）為封閉方程組，可利用其進(jìn)行求解。

通過(guò)計(jì)算可得到面源、線源和點(diǎn)源裝藥強(qiáng)爆炸自模擬數(shù)值解，結(jié)果如圖1～圖3所示。圖中橫、縱坐標(biāo)的數(shù)值區(qū)間（0～1）均為無(wú)量綱形式，其中橫坐標(biāo)的相對(duì)位置，即r/rs，為無(wú)量綱形式的距爆炸中心的距離。

從圖1～圖3可以看出，面源、線源和點(diǎn)源裝藥強(qiáng)爆炸解參數(shù)隨波陣面相對(duì)位置分布的特征基本相同，均呈現(xiàn)類(lèi)似的變化趨勢(shì)。質(zhì)點(diǎn)速度隨著r/rs下降基本上呈線性下降趨勢(shì)；壓力在r/rs從1降低到0的過(guò)程中先下降比較快，當(dāng)r/rs降到0.5之后直到爆炸中心，壓力基本穩(wěn)定在0.36～0.37；密度隨r/rs降低很快接近至0。

試驗(yàn)表明，自模擬解 [7]可以很好地描述爆炸的初始階段。因此，利用真實(shí)空氣條件下面源、線源和點(diǎn)源裝藥強(qiáng)爆炸的自模擬解，作為T(mén)NT裝藥爆炸近區(qū)流場(chǎng)特性數(shù)值計(jì)算的初始條件是精確合理的。

二、球形TNT裝藥爆炸近區(qū)流場(chǎng)特性

（一）爆炸流場(chǎng)物理模型

選取爆轟波傳至裝藥表面的時(shí)刻為初始計(jì)算時(shí)刻，爆轟波的波陣面將整個(gè)流場(chǎng)分為兩個(gè)部分（見(jiàn)圖4）。

根據(jù)守恒定律，在激波陣面上必須滿足如下關(guān)系：

質(zhì)量守恒：

動(dòng)量守恒：

圖1 面源裝藥強(qiáng)爆炸解參數(shù)分布圖

圖2 線源裝藥強(qiáng)爆炸解參數(shù)分布圖

圖3 點(diǎn)源裝藥強(qiáng)爆炸解參數(shù)分布圖

能量守恒：

式（24）～式（26）中，P0、u0、ρ0、E0分別表示激波前未擾動(dòng)空氣的壓力、質(zhì)點(diǎn)速度、密度和內(nèi)能；用P、u、ρ、E表示激波陣面后擾動(dòng)空氣的壓力、質(zhì)點(diǎn)速度、密度和內(nèi)能；D為沖擊波速度。

圖4 初始計(jì)算時(shí)刻流場(chǎng)分布示意圖

對(duì)式（24）～式（26）無(wú)量綱化后，可以解出單值對(duì)應(yīng)的激波關(guān)系。對(duì)偏微分方程組（9）進(jìn)行差分近似，為了保持接觸間斷面處和激波跳躍間斷面處計(jì)算的貫通，采用人工粘性法，在出現(xiàn)壓力項(xiàng)處附加人工粘性項(xiàng)ζ，并將其無(wú)量綱化：ρ = ζ/η，于是將方程組（9）（N = 2時(shí)為球形裝藥爆炸）變?yōu)闊o(wú)量綱差分格式：

式（27）中，l為差分網(wǎng)格整格點(diǎn)，n為時(shí)間整步長(zhǎng)，e為空間半格l–1/2和時(shí)間整格點(diǎn)n + 1處。

其中粘性項(xiàng)為：

關(guān)于能量方程的差分，由于參數(shù)比較復(fù)雜，因此，采用新的方程式，從無(wú)量綱狀態(tài)方程出發(fā)，將e進(jìn)行代換可得到新的形式：

式（28）中，ξ = lnη。

從而得到能量的差分格式：

因上式是一個(gè)隱式的差分格式，因此采用二次迭代進(jìn)行求解。建立激波后裝配方法：首先進(jìn)行激波探查，用三次樣條插值求ξ的最大值點(diǎn)，將此位置定為激波陣面位置λs，求得激波陣面上的壓力，然后調(diào)用波陣面參數(shù)的求解子程序求得波陣面上的其他參數(shù)。數(shù)值計(jì)算的時(shí)間步長(zhǎng)由Courant條件和激波區(qū)域的差分方程穩(wěn)定性條件共同確定。

（二）算例計(jì)算與結(jié)果

1. 空氣和TNT裝藥參數(shù)

靜止空氣的參數(shù)見(jiàn)表1。

表1 靜止空氣參數(shù)

TNT裝藥參數(shù)見(jiàn)表2。

表2 TNT裝藥參數(shù)

2. 結(jié)果

采用本文建立的計(jì)算方法可得到下述結(jié)果：

（1）超壓

圖5為無(wú)量綱峰值超壓ΔPs、峰值動(dòng)壓qs隨無(wú)量綱激波陣面距離位置的變化曲線，從圖中可看出，本文的計(jì)算結(jié)果與Baker [11]的計(jì)算結(jié)果吻合度較好。

（2）接觸間斷面和主激波陣面的時(shí)–空曲線

圖6為主激波（MS）、二次激波（2S）以及接觸間斷面（CS）無(wú)量綱位置λ隨無(wú)量綱時(shí)間τ的變化規(guī)律。從圖中可看到Friedman [12]發(fā)現(xiàn)的二次激波現(xiàn)象，而出現(xiàn)的多個(gè)后續(xù)激波是由于波的反射產(chǎn)生的。從圖6可得到?jīng)_擊波的影響范圍遠(yuǎn)遠(yuǎn)大于爆轟產(chǎn)物接觸間斷面的擴(kuò)散范圍。

（3）峰值超壓、峰值動(dòng)壓隨激波陣面距離變化曲線

圖7和圖8分別為不同裝藥量球形TNT裝藥在真實(shí)空氣中爆炸產(chǎn)生沖擊波的主激波峰值超壓及峰值動(dòng)壓隨激波陣面距離的變化曲線圖（采用的是指數(shù)坐標(biāo)）。由圖7、圖8可知，對(duì)于同一裝藥量，球形TNT裝藥隨著距離的不斷增大，超壓峰值和動(dòng)壓峰值不斷下降。不同裝藥量的炸藥爆炸超壓、動(dòng)壓隨距離變化符合相似律關(guān)系。隨著裝藥量的增加，相同位置處的超壓峰值和動(dòng)壓峰值相應(yīng)增加。

圖5 峰值超壓、峰值動(dòng)壓隨激波半徑的變化曲線（實(shí)心點(diǎn)為本文結(jié)果，空心點(diǎn)為Baker的計(jì)算結(jié)果）

圖6 沖擊波與接觸間斷面的時(shí)–空曲線

（4）沖擊波到達(dá)時(shí)間隨激波陣面距離變化的曲線

圖9為不同裝藥量球形TNT裝藥在真實(shí)空氣中爆炸產(chǎn)生的沖擊波到達(dá)時(shí)間隨激波陣面距離的變化規(guī)律。可見(jiàn)隨著裝藥量的增加，球形TNT裝藥爆炸沖擊波到達(dá)相同位置處的時(shí)間不斷降低。不同裝藥量的沖擊波到達(dá)時(shí)間變化規(guī)律大致相同，可見(jiàn)同樣滿足相似律關(guān)系。

（5）沖擊波超壓沖量隨激波陣面距離變化的曲線

圖8 不同裝藥量球形TNT裝藥爆炸峰值動(dòng)壓隨激波陣面距離變化的曲線

圖9 不同裝藥量球形TNT裝藥爆炸沖擊波到達(dá)時(shí)間隨激波陣面距離變化的曲線

圖10為不同裝藥量球形TNT裝藥爆炸沖擊波超壓沖量隨激波陣面距離的變化曲線。從圖中可知，對(duì)于同一裝藥量，球形TNT裝藥隨著激波陣面距離的不斷增大，沖量先逐漸降低，當(dāng)?shù)竭_(dá)一定距離后，在此區(qū)域會(huì)出現(xiàn)一個(gè)平緩甚至稍微增加的階段，之后再開(kāi)始衰減，這是由于二次激波追趕上了主激波使得正壓作用沖量增加造成的；另一方面，隨著裝藥量的增加，相同位置處的沖量也逐漸增加。

（6）沖擊波的正壓作用時(shí)間隨激波陣面距離變化的曲線

圖11是不同裝藥量的球形TNT裝藥空中爆炸沖擊波的正壓作用時(shí)間隨激波陣面距離變化的曲線。從圖中可知，對(duì)于同一裝藥量，爆炸沖擊波隨激波陣面距離的變化規(guī)律與圖10類(lèi)似，正壓作用時(shí)間也存在一個(gè)轉(zhuǎn)折點(diǎn)，它先隨距離的增加而降低，但在某個(gè)位置處又開(kāi)始隨距離的增大而增大。這主要是因?yàn)榇嬖诙渭げǖ脑?，隨著沖擊波不斷向外傳播，主激波的峰值超壓不斷降低，且主激波越晚到達(dá)的點(diǎn)的正壓作用時(shí)間就會(huì)越少，但是由于二次激波的不斷追趕，當(dāng)?shù)竭_(dá)某個(gè)位置后，主激波的正壓還未衰減到負(fù)壓，而二次激波就追趕上來(lái)，從而使得正壓作用持續(xù)時(shí)間增加。對(duì)于1～50 kg的TNT裝藥，轉(zhuǎn)折點(diǎn)的位置彼此相差不大，基本在4 m左右；而對(duì)于100～1 000 kg的大裝藥量藥球，轉(zhuǎn)折點(diǎn)離爆心的距離隨裝藥量的增大而增大。另外，對(duì)于裝藥量為1～50 kg的藥球，在固定的位置點(diǎn)，轉(zhuǎn)折點(diǎn)前后均有正壓作用時(shí)間隨裝藥量的增大而增大的趨勢(shì)；但對(duì)于裝藥量為100～1 000 kg的藥球，在轉(zhuǎn)折點(diǎn)前，正壓作用時(shí)間同樣隨裝藥量的增大而增大，但是在轉(zhuǎn)折點(diǎn)后則相反，正壓作用時(shí)間隨裝藥量的增加反而降低（見(jiàn)圖11中的箭頭顯示）。

圖10 不同裝藥量球形TNT裝藥爆炸沖擊波超壓沖量隨激波陣面距離變化的曲線

三、結(jié)論

由于空氣中爆炸的近區(qū)空氣與理想氣體之間存在較大差異，因此，本文采用真實(shí)空氣狀態(tài)方程來(lái)描述爆炸近區(qū)流場(chǎng)的特性?；诒疚乃⒌睦碚?，結(jié)合人工粘性技術(shù)和激波裝配技術(shù)，對(duì)TNT裝藥真實(shí)空氣中爆炸流場(chǎng)參數(shù)進(jìn)行了計(jì)算，結(jié)果與已有文獻(xiàn)研究結(jié)果相吻合，并得到以下結(jié)論：

圖11 不同裝藥量球形TNT裝藥爆炸沖擊波正壓作用時(shí)間隨激波陣面距離變化的曲線

（1）TNT裝藥爆炸之后會(huì)產(chǎn)生一個(gè)初始超壓為60 MPa的強(qiáng)沖擊波（圖5中無(wú)量綱波陣面位置λs= 1對(duì)應(yīng)的峰值超壓），強(qiáng)度隨沖擊波的不斷傳播而逐漸衰減，并且在主激波之后存在一個(gè)不斷向外傳播的接觸間斷面。

（2）爆炸近區(qū)最典型的物理特征是存在一個(gè)爆轟產(chǎn)物區(qū)，在半徑小于3倍裝藥半徑的區(qū)域內(nèi)，壓力強(qiáng)度遠(yuǎn)遠(yuǎn)大于主激波的壓力值。

（3）當(dāng)爆轟波到達(dá)裝藥與空氣的接觸面時(shí)，在形成向外傳播的空氣沖擊波的同時(shí)，還會(huì)形成一個(gè)向內(nèi)傳播的激波，向內(nèi)傳播的激波在爆心處發(fā)生碰撞后會(huì)得到加強(qiáng)，然后繼續(xù)向外傳播，在到達(dá)接觸面后又形成一個(gè)新的向內(nèi)傳播和向外傳播的激波，從而形成定點(diǎn)超壓波形負(fù)壓區(qū)的二次、三次激波。

由此可見(jiàn)，本文提出的計(jì)算分析方法可以很好地模擬和描述TNT裝藥爆炸近區(qū)在真實(shí)空氣中的傳播規(guī)律。本文的研究成果可以對(duì)爆炸力學(xué)和防護(hù)工程研究起到重要的參考作用。

[1] 盧芳云. 一維不定常流體動(dòng)力學(xué)教程 [M]. 北京∶ 科學(xué)出版社,2006.

Lu F Y. Tutorial of one-dimensional unsteady fl ow dynamics [M].Beijing∶ China Science Publishing &Media Ltd., 2006.

[2] Taylor G. The formation of a blast wave by a very intense explosion. I. Theoretical discussion [J]. Proceeding of the Royal Society of London. 1950 (201)∶ 159–174.

[3] Von Neumann J. The point source solution [R]. Washington, D.C.∶National Defense Research Committee, 1947.

[4] Bach G G, Lee J H S. An analytical solution for blast waves [J].AIAA Journal, 1970 (8)∶ 271–275.

[5] Sedov L I. Similarity and dimensional methods in mechanics [M].New York∶ CRC Press, 1965.

[6] 湯文輝. 沖擊波物理 [M]. 北京∶ 科學(xué)出版社, 2011.

Tang W H. Shock wave physics [M]. Beijing∶ China Science Publishing &Media Ltd., 2011.

[7] Л.П.奧爾連科. 爆炸物理學(xué)(原書(shū)第三版) [M]. 孫承緯譯. 北京∶ 科學(xué)出版社, 2011.

Л.П.ОРЛЕНКО. Explosion physics [M]. Translated by Sun C W. Beijing∶ China Science Publishing &Media Ltd., 2011.

[8] 喬登江. 核爆炸物理概論 [M]. 北京∶ 原子能出版社, 1988.Qiao D J. Introduction to nuclear explosion physics [M]. Beijing∶Atomic energy Press, 1988.

[9] 周豐峻, 陳葉青, 任輝啟. 爆炸近區(qū)空氣沖擊波規(guī)則反射和非規(guī)則反射 [C].應(yīng)用力學(xué)進(jìn)展論文集. 北京∶ 中國(guó)力學(xué)學(xué)會(huì), 1999∶1–12.

Zhou F J, Chen Y Q, Ren H Q. The regular reflection and nonregular re fl ection of shockwave in the close-in fi eld of TNT charge explosion [C]. Proceedings of applied mechanics progress. Beijing∶The Chinese Society of Theoretical and Applied Mechanics, 1999∶1–12.

[10] Brode H L. Blast wave from a spherical charge [J]. Physics of Fluids, 1959, 2 (2)∶ 217–229.

[11] Baker W E. Explosives in air [M]. Austin∶ University of Texas Press, 1973.

[12] Friedman M P. A simplified analysis of spherical or cylindrical blast waves [J]. Journal of Fluid Mechanics, 1961, 11 (1)∶ 1–15.