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

趙宏濤，栗建橋，馬天寶

(北京理工大學(xué) 爆炸科學(xué)與技術(shù)國家重點(diǎn)實(shí)驗室， 北京 100081)

1 引言

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

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

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

2 理論模型

2.1 爆炸的電離模型

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

(1)

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

可以得出電導(dǎo)率計算公式為:

(2)

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

2.2 MHD模型

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

(3)

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

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

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

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

(4)

因為在前面的計算中已經(jīng)得到了磁感應(yīng)強(qiáng)度，僅考慮無耗無源情況，故需要考慮的Maxwell方程變?yōu)椋?/p>

(5)

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

3 數(shù)值模擬

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

(6)

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

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

(7)

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

(8)

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

(9)

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

對比前后2個時間步內(nèi)得到的磁感應(yīng)強(qiáng)度，可得到磁感應(yīng)強(qiáng)度的變化率，在網(wǎng)格上進(jìn)行積分為：

(10)

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

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

p=(γ-1)ρe

(11)

炸藥和空氣的參數(shù)如表1所示，其中ρ0為初始密度，e0為初始比內(nèi)能，γ為比熱容比，T0為初始溫度，μe為磁導(dǎo)率。

表1 數(shù)值模擬材料參數(shù)

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

表2 炸藥質(zhì)量及尺寸

圖1 數(shù)值模擬計算模型示意圖

4 數(shù)值模擬結(jié)果分析

4.1 電磁場變化情況

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

圖2 電場強(qiáng)度理論、數(shù)值模擬和實(shí)驗結(jié)果曲線

圖3給出了800 g B炸藥爆轟完全后30 μs內(nèi)x方向磁感應(yīng)強(qiáng)度Bx的演化過程。

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

圖3 不同時刻x方向磁感應(yīng)強(qiáng)度Bx數(shù)值模擬云圖

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

圖4 50 μs數(shù)值模擬云圖

圖5給出了圖1中關(guān)鍵點(diǎn)A-E上的磁感應(yīng)強(qiáng)度Bx隨時間和距離變化的情況。關(guān)鍵點(diǎn)A為爆心處，關(guān)鍵點(diǎn)E離爆心最遠(yuǎn)。

圖5 各關(guān)鍵點(diǎn)處磁場擾動曲線

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

圖6 關(guān)鍵點(diǎn)B處的電磁場變化曲線

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

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

4.2 炸藥當(dāng)量對爆炸產(chǎn)生的電磁場的影響

不同當(dāng)量炸藥爆炸產(chǎn)生磁場的數(shù)值模擬結(jié)果如圖7所示。不同當(dāng)量炸藥爆炸產(chǎn)生電場的數(shù)值模擬結(jié)果如圖8所示。

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

圖7(b)是1.5 m處磁感應(yīng)強(qiáng)度峰值隨炸藥當(dāng)量之間的關(guān)系。取擬合關(guān)系式為：

Bx=a+b×wm

(12)

從擬合結(jié)果來看，當(dāng)量指數(shù)m為0.633±0.067，其相關(guān)性為0.932。在此擬合基礎(chǔ)上，可認(rèn)為磁場擾動峰值與炸藥當(dāng)量的2/3次方成正比。

圖7 不同當(dāng)量炸藥爆炸產(chǎn)生磁場的數(shù)值模擬曲線

圖8 不同當(dāng)量炸藥爆炸產(chǎn)生電場的數(shù)值模擬曲線

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

Ez=a+b×wn

(13)

從擬合結(jié)果來看，當(dāng)量指數(shù)n為0.341±0.031，其相關(guān)性為0.918。在此擬合基礎(chǔ)上，可認(rèn)為磁場擾動峰值與炸藥當(dāng)量的1/3次方成正比，這也和文獻(xiàn)[9]中的實(shí)驗結(jié)果相一致。

5 結(jié)論

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

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

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