炸藥爆炸過程中電磁輻射數值模擬研究

趙宏濤，栗建橋，馬天寶

(北京理工大學 爆炸科學與技術國家重點實驗室， 北京 100081)

1 引言

炸藥爆炸伴隨著電磁輻射現象的產生，這種現象不僅會對爆炸場的信號測量產生干擾，嚴重時還會對測量系統和設備產生破壞。戴晴等[1]用150 g鋁鎂工質進行實驗，其爆炸產生的電磁輻射使5 m外的計算機死機。常規炸藥爆炸產生電磁信號這一現象最早發現于1954年，由Kolsky發表于Nature[2]，之后不斷有學者進行相關的實驗研究。1958年，Cook[3]針對常規炸藥爆炸產生電磁波這一現象，提出了一種解釋機理，他認為是導電的爆轟產物在地球背景電磁場中產生了極化，運動至接地后放電產生電磁波。Boronin等[4-5]進行了針對性實驗，來驗證這個解釋。在之后仍然有大量的研究，主要集中在爆炸產生的電磁信號的測量和豐富實驗成果上。2001年Fine[6]提出了更為直觀的理論解釋，他認為爆炸產生的等離子體在激波前沿沖擊波厚度內被加速引起了電磁輻射，給出了詳細的計算公式，可以用來估算爆炸產生電磁場強度的量級，并給出了大量的實驗數據。2014年Kuhl[7]提出了一種新的理論模型，Kuhl等認為等離子體在背景磁場中的集群運動會產生電流，而這個電流會導致周圍空間電磁場的變化，他們根據這一理論模型建立了一維爆炸問題模型，模擬了球形爆轟產生的電磁輻射現象。

國內也在這方面進行了一系列研究，但相比于國外的研究，國內的研究起步較晚，研究也主要集中在實驗方面。1997年，陳生玉等[8]測量了帶殼裝藥爆炸引起的電磁輻射，用量綱分析法對實驗測得的數據進行分析，得出了電磁脈沖的幅值隨裝藥量增大而增大等結論，分析了電場強度和各個爆炸參數之間的關系。2011年，曹景陽等[9]對聚能炸藥索爆炸引起的電磁干擾進行了測量，分析實驗結果得到爆炸索爆炸后的數十毫秒內都會有明顯的電磁脈沖出現，其頻率主要集中在兆赫茲量級，在單個頻點上的電場峰值可達數個V/m。王長利等[10]對B炸藥和梯黑鋁炸藥爆炸產生的電磁效應進行測量，結果表明炸藥爆炸產生的電磁輻射信號強度與炸藥當量的1/3次方成正比，隨距離增加而衰減。幾種典型炸藥的爆炸電磁輻射信號主要頻率在100 MHz以內，同一種炸藥信號頻譜特征與當量、爆心距無關，不同種類炸藥產生信號頻譜不同。陳鴻等[11]對RDX基含鋁炸藥爆炸電磁輻射信號特性進行了實驗研究。結果表明，RDX以及RDX基含鋁炸藥爆炸電磁輻射信號發生時刻與起爆時刻相比有明顯延遲。距爆心2 m處，爆炸電磁輻射信號強度在1.87～15.20 V/m范圍內，隨距離的增加而衰減。目前國內外對炸藥爆炸產生電磁輻射的認識仍然不夠清晰，對炸藥爆炸過程中電磁輻射的產生機理和影響因素缺乏全面的了解。因此本研究能加深對炸藥爆炸過程中產生電磁輻射問題的認識，揭示該過程中電磁輻射產生機理與時空演化規律，并在爆炸測試電磁防護方面有一定的工程應用價值。

由于國內外在炸藥爆炸產生電磁輻射領域內的數值模擬研究很少，因此對其進行數值模擬研究是十分必要的。本文將爆炸流體計算方法、熱平衡電離方程、非理想磁流體動力學模型和麥克斯韋方程相結合，通過基于AUSM+-up通量分裂方法的有限體積法進行炸藥爆炸產生電磁輻射數值模擬，并針對不同當量的裝藥進行了數值模擬研究，給出了當量對炸藥爆炸產生電磁輻射的影響規律。

2 理論模型

2.1 爆炸的電離模型

從Lee[12]和栗建橋[13]等的研究可知，炸藥爆炸過程中產生的等離子體處于熱平衡狀態，因此選擇熱平衡電離模型作為炸藥爆炸產生等離子體的模型。在考慮了體系勢能、溫度和電離程度后，文獻[12]給出了炸藥爆炸產生等離子體的電離模型為：

(1)

式中：ne為電子的數密度；me為電子質量；k為玻爾茲曼常數；T為氣體溫度；h為普朗克常數；p為各個組分；Ip為組分p的一階電離能。

可以得出電導率計算公式為:

(2)

式中：q為原電荷數；d為中性粒子的有效直徑；n為電離氣體所有粒子的數密度。

2.2 MHD模型

相關文獻[7]分析認為，炸藥爆炸產生的電磁輻射是由高溫高壓的導電爆轟產物在自然磁場中運動而產生的，因此采用磁流體動力學(MHD)模型來描述等離子體的運動和電磁場的演化，以實現對爆炸過程中電磁輻射的數值模擬。結合磁擴散后的MHD控制方程如式(3)所示。

(3)

式中：ρ為流體密度；v為速度；B為磁感應強度；e為比內能；σ為電導率；ε為介電常數；μ為磁導率；p為壓強。從式(3)可以看出，電導率σ對最后一個方程有很大的影響。當電導率很高時，磁感應強度的時間和空間二階導數都很小，磁場完全凍結在等離子體中，并隨著等離子體的運動而運動，這對應著理想情況下的磁流體動力學模型。當等離子體的電導率和介電常數相當時，不僅要考慮到磁場隨等離子體運動產生的變化，還需要考慮到磁擴散效應，電磁波效應會比較顯著。

從上述分析可以看出，炸藥爆炸過程中，爆轟產物電離產生的等離子體在地磁場中運動是爆炸過程中電磁輻射的根源。而等離子體的電導率是影響爆炸過程中電磁輻射特性的唯一因素，因此將2.1節得到的等離子體電導率代入全MHD模型，求解即可得到爆炸過程中的電磁輻射特性。由于爆炸過程中電導率空間分布不均勻，需要將其分為理想磁流體和非理想磁流體2種情況進行分析。

在理想磁流體情況下，忽略磁擴散和電磁波效應，由于磁場完全被凍結在等離子體中，因此磁場演化的時間尺度和流場運動的時間尺度完全一致。只需要對控制方程(3)中的前6個公式求解，即可得到理想磁流體情況下的各個物理量的變化。

在非理想磁流體情況下，由于電導率相對較低時的磁擴散和電磁波效應較為顯著，需要同時考慮磁擴散、電磁波效應和求解穩定性條件，因此求解過程中的時間步長非常短，與流場求解時的時間步長存在量級差距，為了解決這個問題，求解過程被拆成2步。首先計算理想磁流體控制方程，求得并更新所有求解變量，再通過隱式時間積分法，在這個時間步內應用隱式算法求解方程(4)，之后再單獨更新磁感應強度。

(4)

因為在前面的計算中已經得到了磁感應強度，僅考慮無耗無源情況，故需要考慮的Maxwell方程變為：

(5)

式中：E為電場強度。只需要求解方程(5)即可得電場強度。

3 數值模擬

隨著仿真技術的發展，爆炸與高速沖擊等強間斷非線性問題出現了很多解決方法[14-16]。本文采用基于AUSM+-up通量分裂方法[18]的有限體積法對全MHD控制方程進行計算，給出通量描述如下：

(6)

式中：Fc為守恒通量；Fp為壓力通量；FM為MHD通量；ρ為流體密度；U為速度；μe為磁導率，其余各量在式(3)中已給出計算公式。

AUSM+-up通量分裂只涉及理想磁流體方程的求解，對于電磁波和磁擴散的求解在后面用隱式方法給出。對于理想磁流體，將通量分解為守恒部分Fc、壓力部分Fp和MHD部分FM，對于前兩部分，采用標準的AUSM+-up通量分裂方法即可，對于FM，可按照粘性通量格式進行處理，即：

(7)

式中：下標g表示單元界面通量；i表示界面左單元；j表示界面右單元，最終界面通量為：

(8)

對于磁擴散技術過高的情況，需要在一個顯式時間積分步內，采用隱式方法求解?？紤]時間離散后，給出完備的離散格式為：

(9)

歷遍所有單元界面，組成線性方程組，通過MKL中的PardisoLU分解稀疏矩陣，求解器直接求得最終的磁感應強度。

對比前后2個時間步內得到的磁感應強度，可得到磁感應強度的變化率，在網格上進行積分為：

(10)

為了盡量減少影響因素，開展了單純的空中爆炸數值模擬，忽略地面反射，炸藥為正方形。在空中爆炸的情況下，可以定量分析產生的電磁場參數和炸藥當量之間的關系。以自然磁場作為背景磁場，進行爆炸產生的電磁輻射數值模擬研究。取自然磁場磁感應強度為0.5 Gs，其余變量取國際單位制，計算時間取0.2 ms。

炸藥采用瞬時起爆模型，狀態方程采用理想氣體狀態方程，有：

p=(γ-1)ρe

(11)

炸藥和空氣的參數如表1所示，其中ρ0為初始密度，e0為初始比內能，γ為比熱容比，T0為初始溫度，μe為磁導率。

表1 數值模擬材料參數

二維的計算網格模型如圖1所示，計算域為邊長10 m的正方形，中心位置放置炸藥。不同當量(50 g、200 g、800 g和3.2 kg)的B炸藥尺寸如表2所示，網格劃分采用中心密集，邊界稀疏，保證計算能準確描述爆炸初始階段各個物理量變化率較大的情況。整體的網格采用非結構網格，在邊界處約束磁場強度為初始值，設置初始磁場為x方向，大小為0.5 Gs，對流動過程的邊界設置為流出邊界。計算過程可以獲得模型關鍵點A-F的時程曲線如圖1所示。

表2 炸藥質量及尺寸

圖1 數值模擬計算模型示意圖

4 數值模擬結果分析

4.1 電磁場變化情況

圖2(b)給出了3.2 kg B炸藥爆炸時產生在爆心處電場擾動的數值模擬結果。由于在計算中采用了瞬時爆轟模型，忽略了起爆到炸藥完全爆轟的時間，因此在模擬計算一開始便產生了電場擾動信號，考慮到炸藥的邊長為13.57 cm，其對角線的一半長度為9.6 cm，從起爆到完全爆轟約為10 μs，因此在數值模擬結果中將-10 μs作為起爆時間，信號的出現時間與文獻[7]符合較好。從整體的波形來看，數值模擬結果的波形和Kual[7]給出的理論計算結果相似，同時和任會蘭[17]給出的實驗結果在電場脈沖持續時間和波形上相似，由此可以證明本文中的數值模擬結果的可靠性。

圖2 電場強度理論、數值模擬和實驗結果曲線

圖3給出了800 g B炸藥爆轟完全后30 μs內x方向磁感應強度Bx的演化過程。

從圖3數值模擬的結果來看，在炸藥爆轟完全瞬間，由于等離子體的運動造成了磁場的劇烈變化，爆炸中心處的磁場出現了明顯的下降，隨著爆轟產物的運動，在y方向上形成了2個明顯的磁場峰值區域，并隨著時間的推移，這2個磁場峰值區域之間的距離逐漸增大。這主要是因為爆轟產物運動初期，由于高溫高壓爆轟產物能量密度高，形成的等離子體電導率高，處于理想磁流體情況下，磁場被凍結在等離子體內，與等離子體即爆轟產物一起運動，因此在y方向上產生了2個明顯的峰值區域，并隨著爆轟產物一起運動。

圖3 不同時刻x方向磁感應強度Bx數值模擬云圖

圖4則給出了800 g B炸藥爆轟完全后50 μs時x方向的流動速度、電導率、比內能及x方向磁感應強度Bx的空間分布。從上一節的理論可知，電導率的分布主要取決于流場中的能量分布，從圖4(b)、圖4(c)中也可以看出，電導率的峰值主要出現在流場速度和比內能最大的地方，這個結果與理論互相印證。對比圖4(a)、圖4(b)和圖4(c)，可以發現：電導率與比內能高的地方正好是爆轟產物運動的前沿，在這些地方，由于高電導率導致磁擴散度很小，磁場被凍結在等離子體內，磁場中的擾動可以隨時間積累起來。而其他地方電導率極低，磁擴散度很大，磁場很難形成明顯的積累，產生的磁場擾動以很快的速度向周圍擴散，從而難以形成明顯的磁場擾動積累。

圖4 50 μs數值模擬云圖

圖5給出了圖1中關鍵點A-E上的磁感應強度Bx隨時間和距離變化的情況。關鍵點A為爆心處，關鍵點E離爆心最遠。

圖5 各關鍵點處磁場擾動曲線

從圖5(a)可以看出，炸藥爆炸后磁感應強度先增長，在到達峰值時刻后，磁場擾動逐漸降低，到后來趨于平衡。這是因為在爆轟產物運動初期，由于其能量密度高、電導率高、磁擴散度低，磁場被凍結在爆轟產物內，產生的磁場擾動可隨時間進行積累。而在爆轟產物膨脹到一定程度后，由于能量密度的下降導致電導率顯著降低，從而引起磁擴散度升高，磁場積累的擾動逐漸向外擴散出去。從圖5(b)可以看出，關鍵點離爆心越遠，磁場擾動峰值越小，這符合電磁場擴散中的能量守恒原理。圖6為800 g B 炸藥爆炸后關鍵點B處的電場和磁場變化情況。

圖6 關鍵點B處的電磁場變化曲線

從圖6可以看出，磁感應強度先上升后下降并逐漸趨于平衡。電場強度先上升后下降，下降至零后有向負向的增長，增長至第2個峰值后逐漸減少并趨于平衡。

電場峰值出現在磁場上升沿斜率最大處，即磁場變化最劇烈的時刻，而當磁場峰值出現后，電場正好有正轉負。這符合感生電場與磁場之間的規律。由于磁場的變化為上升沿陡，峰值之后的下降速度比上升慢，因此導致電場在負向的峰值小于正向。隨著時間的變化，電場和磁場都最終趨于平衡。

4.2 炸藥當量對爆炸產生的電磁場的影響

不同當量炸藥爆炸產生磁場的數值模擬結果如圖7所示。不同當量炸藥爆炸產生電場的數值模擬結果如圖8所示。

圖7(a)為對不同當量(50 g、200 g、800 g、3.2 kg)炸藥爆炸在1.5 m處產生磁場的數值模擬結果。從結果來看，炸藥當量越大，產生的磁場擾動峰值越大，脈沖持續時間越長，磁場擾動的峰值時間也有所推遲。這主要是由于當藥量增大時，爆炸氣體產物維持高能量狀態的時間也相應增大，電導率維持在較高水平的時間加長，使產生的磁場擾動積累時間增加，相應的磁場峰值也會增加，磁場向外擴散的時間增大。

圖7(b)是1.5 m處磁感應強度峰值隨炸藥當量之間的關系。取擬合關系式為：

Bx=a+b×wm

(12)

從擬合結果來看，當量指數m為0.633±0.067，其相關性為0.932。在此擬合基礎上，可認為磁場擾動峰值與炸藥當量的2/3次方成正比。

圖7 不同當量炸藥爆炸產生磁場的數值模擬曲線

圖8 不同當量炸藥爆炸產生電場的數值模擬曲線

圖8(a)為對不同當量(50 g、200 g、800 g、3.2 kg)B炸藥爆炸在1.5 m處產生的電場擾動進行數值模擬的結果。從數值模擬的結果來看，電場與磁場擾動隨炸藥當量的變化趨勢相同，即當量越大，電場擾動峰值越大，峰值到達時間略有推遲，具體原因分析與磁場相似，在此不再贅述。取擬合公式為：

Ez=a+b×wn

(13)

從擬合結果來看，當量指數n為0.341±0.031，其相關性為0.918。在此擬合基礎上，可認為磁場擾動峰值與炸藥當量的1/3次方成正比，這也和文獻[9]中的實驗結果相一致。

5 結論

1) 基于本文提出的磁流體動力學與麥克斯韋方程結合的方法，通過數值求解可以獲得炸藥爆炸過程中電磁輻射的結果，具有參考價值。

2) 常規炸藥爆炸會產生電磁輻射，其中磁場信號只有一個峰值，上升沿比下降沿陡。電場信號有2個峰值，正向的峰值大于負向的峰值。電場信號和磁場信號持續時間相當。電場信號峰值出現時間為磁場信號上升沿中部，為整個爆轟過程中磁場變化最劇烈的時刻。

3) 對于同一種炸藥而言，當炸藥當量增大時，其爆炸產生磁信號峰值、電信號峰值、信號持續時間和信號峰值到達時間都隨之增大，磁信號峰值與當量的2/3次方成正比，電信號峰值與當量1/3次方成正比。