早期高空核電磁脈沖的自洽粒子模擬

付梅艷, 張茂鈺

(1. 北京大學(xué) 數(shù)學(xué)科學(xué)院, 北京 100871; 2. 西北核技術(shù)研究所， 西安 710024)

早期高空核電磁脈沖(early-time high altitude electromagnetic pulse, E1-HEMP)是一種強度大、頻率高、頻譜寬的電磁脈沖，危害性很大[1]。20世紀(jì)60年代，Karzas和Latter討論了其產(chǎn)生機理[2-3]，認為在其產(chǎn)生過程中有一小部分核爆炸能量，經(jīng)過幾個中間步驟轉(zhuǎn)換為射頻電磁能。轉(zhuǎn)換的第1步是核爆炸釋放出γ射線；第2步是γ射線與大氣或其他物質(zhì)相互作用。事實上，高空核爆時，釋放大量的瞬發(fā)γ光子(簡稱γ光子)。γ光子在向外輻射的過程中與周圍稀薄的大氣相互作用散射出康普頓電子。散射出的康普頓電子以接近光速的速度徑向向外運動，形成康普頓電流。康普頓電子在向外運動的過程中，受到地磁場的作用，運動軌跡發(fā)生偏轉(zhuǎn)。因此，所形成的康普頓電流除了徑向分量外，還有明顯的橫向分量，從而激勵出強電磁脈沖。同時，康普頓電子在大氣中遷移時，會不斷與空氣發(fā)生碰撞損失能量，使空氣發(fā)生電離產(chǎn)生次級電子-離子對。康普頓電子受到電磁脈沖場的洛倫茲力的作用，形成新的康普頓電流；次級電子-離子對受到電磁脈沖場中電場力的作用發(fā)生漂移，形成漂移電流，它將抑制場強的增長，與電磁場方程中的電導(dǎo)率相對應(yīng)。這些電流源作為源項，激發(fā)新的電磁脈沖場?？梢姡珽1-HEMP的產(chǎn)生是電磁脈沖場與康普頓電子相互作用的自洽過程。

對E1-HEMP產(chǎn)生過程的數(shù)值模擬，過去常采用非自洽的處理方法[4-5]，即僅考慮地磁場對康普頓電子的作用，不考慮E1-HEMP對康普頓電子的影響。在一定理論分析和假設(shè)近似的基礎(chǔ)上，用這種處理方法可以較為簡單和容易地得到康普頓電流源和電導(dǎo)率，進而求解麥克斯韋方程獲得電磁脈沖場。Karzas和Latter認為E1-HEMP對康普頓電子運動的作用發(fā)生在電流脈沖后期，只會影響E1-HEMP的時域波形，而不會顯著改變發(fā)生在早期的、備受關(guān)注的峰值電場強度[3]。美國空軍武器實驗室在20世紀(jì)70年代曾編制過代碼CHEMP[6]，用來檢驗E1-HEMP對康普頓電子的影響，但是沒有明確的結(jié)論。王建國等通過比較帶電粒子受到電場力和地磁場力大小，說明了E1-HEMP對康普頓電子影響的重要性[1]。程引會等利用等效電導(dǎo)率的概念分析了自洽效應(yīng)，結(jié)果顯示考慮自洽效應(yīng)可使電場強度比不考慮自洽效應(yīng)得到的電場強度減小約10%[7]。

本文針對一維近似模型，采用全電磁粒子模擬方法[8-9]，數(shù)值模擬了E1-HEMP產(chǎn)生的自洽過程，介紹了理論及計算模型，并給出了算例。

1 理論模型

1.1 麥克斯韋方程組

國際單位制下Maxwell方程組的2個旋度方程為

(1)

(2)

其中，E，B，J分別為電場強度、磁感應(yīng)強度和電流密度，都是位置和時間的函數(shù)；ε0和μ0分別為真空介電常數(shù)和真空磁導(dǎo)率；σ為電導(dǎo)率。

1.2 康普頓電流

康普頓電流是由康普頓電子的運動形成的。求解康普頓電流的時空分布，需要知道康普頓電子數(shù)密度的時空分布和電子的速度。

1.2.1 康普頓電子的產(chǎn)生

高空核爆炸時，產(chǎn)生大量的γ光子。在分析γ光子在大氣中輻射傳輸并產(chǎn)生康普頓電子的過程時，采取2個假設(shè)近似：1)每一個γ光子發(fā)生康普頓散射只散射出一個康普頓電子。不再考慮散射γ光子。2)γ光子與周圍大氣發(fā)生康普頓散射作用的概率是γ光子在當(dāng)?shù)氐淖杂沙痰牡箶?shù)，即γ光子在一個自由程內(nèi)均勻產(chǎn)生出康普頓電子。

爆點處，γ光子的產(chǎn)生率為[5]

(3)

(4)

根據(jù)射線傳輸理論，在爆點產(chǎn)生的γ 光子將以光速c在大氣中傳輸，且在傳輸過程中γ 光子強度有擴散和吸收衰減[10]。從t時刻開始，在dt時間內(nèi)，到達距離爆心r處的γ 光子數(shù)dNγ(r,t)為

(5)

γ光子在r處產(chǎn)生的康普頓電子數(shù)dNpri(r,t)為

(6)

其中，λ(r)為γ光子在r處的自由程。

1.2.2 康普頓電子的運動

在電磁場和地磁場中，康普頓電子受牛頓-洛倫茲力的作用，速度和位置隨時間發(fā)生變化。當(dāng)康普頓電子的速度與光速相差不大時，其運動方程為相對論牛頓-洛倫茲力方程：

(7)

(8)

其中，m0是康普頓電子的靜止質(zhì)量；v是康普頓電子的速度；γ是相對論因子；B0是地磁場的磁感應(yīng)強度。

1.2.3 康普頓電子的消失

康普頓電子在大氣中遷移與周圍空氣分子、原子不斷發(fā)生碰撞，損失能量。類似γ光子的吸收衰減，康普頓電子數(shù)按指數(shù)規(guī)律衰減：

(9)

其中，L1是康普頓電子開始衰減的位置，即產(chǎn)生的位置；N1是在L1處的康普頓電子數(shù)；L2是考察點位置；N2是在L2處的康普頓電子數(shù)；Re(r)是康普頓電子在r處的最大射程。

1.3 電導(dǎo)率

康普頓電子與大氣分子碰撞時，損失自身能量的同時使空氣分子發(fā)生電離，產(chǎn)生大量的次級電子和正離子?？灯疹D電子平均損失能量33 eV，可產(chǎn)生一對次級電子和離子。次級電子和離子的貢獻是使空氣的電導(dǎo)率σ大大增加。σ可以通過求解空氣化學(xué)方程得到[4]，也可通過分析等離子體參數(shù)得到[11]，計算公式為

(10)

其中，nsec(t)表示二次電子的數(shù)密度；eq表示二次電子電荷的絕對值；msec表示二次電子的靜止質(zhì)量；υ為二次電子的碰撞頻率。

類似于γ光子產(chǎn)生康普頓電子的過程，假設(shè)二次電子由康普頓電子在一個射程內(nèi)均勻產(chǎn)生，可得到二次電子。t時刻，在r處康普頓電子產(chǎn)生二次電子的速率為

(11)

2 計算模型

2.1 宏粒子模型

宏粒子是全電磁粒子模擬方法中的基本概念之一，是指由相空間中一些帶電粒子形成的、具有一定形狀和空間分布的“超級粒子”。宏粒子中所有電子的速度相同。假設(shè)宏粒子包含n個帶電粒子，則宏粒子的電荷為ne，質(zhì)量為nm，其中，e為電子的電荷，m為電子的質(zhì)量。這里，帶電粒子是指康普頓電子。

按照一定的標(biāo)準(zhǔn)，把每一時刻、每個空間位置處新產(chǎn)生的康普頓電子設(shè)置為一定數(shù)目的宏粒子數(shù)。例如，將100個康普頓電子視為一個宏粒子。對計算區(qū)域內(nèi)所有宏粒子的位置、速度及包含的康普頓電子數(shù)進行疊加，即可得到它們產(chǎn)生的初級電流密度。同時，設(shè)置一定的閾值，當(dāng)宏粒子所包含的康普頓電子數(shù)小于該閾值時，則認為宏粒子已經(jīng)衰減殆盡。

2.2 電磁場方程的計算

將式(1)、式(2)在球坐標(biāo)系下展開，略去θ、φ的偏微商項，三維問題可以簡化為一維問題[5]，得到電磁場方程為

(12)

(13)

(14)

(15)

(16)

(17)

采用時域有限差分方法求解上述方程，取Mur一階近似邊界條件，右邊界條件為

左邊界條件為

2.3 相對論牛頓-洛倫茲力方程的求解

采用Boris提出的“半加速-旋轉(zhuǎn)-半加速”方法中的修正公式，求解相對論牛頓-洛倫茲力方程，可得到宏粒子的速度和位置[12]。

2.4 電流密度和電磁場量的分配

為了避免“自力”的產(chǎn)生，電流的分配方法和電磁場量的分配方法一致。本文采用Eastwood提出的方法[13]，即根據(jù)宏粒子運動的起點和終點，計算電流流量，如果宏粒子的運動超過一個網(wǎng)格，則把整個運動分解成一系列網(wǎng)格內(nèi)的運動，然后，把每一網(wǎng)格內(nèi)的電流密度分配到網(wǎng)格線上。對所有宏粒子都以同樣的方式進行處理。

3 算例

3.1 計算條件

選取爆高為100 km，當(dāng)量1 000 kt。假設(shè)1 kt當(dāng)量產(chǎn)生的能量為2.61×1025MeV，爆炸當(dāng)量的0.3%以能量為1.5 MeV的瞬發(fā)γ光子釋放。瞬發(fā)γ源的時間變化函數(shù)取雙指數(shù)函數(shù)，相應(yīng)參數(shù)分別為4.0×107s-1和4.76×108s-1?？灯疹D電子的初速度為v=(vr,vθ,vφ)=(βc, 0, 0), 其中，常數(shù)β=0.914，c為光速。地磁場的簡化近似及坐標(biāo)系的選取方法與文獻[14]一致。

一維計算區(qū)域的方向位于距離爆心投影點正南約170 km處的點與爆心之間的連線方向；計算范圍為rbegin≤r≤rend，其中rbegin為計算區(qū)域內(nèi)半徑，rend為計算區(qū)域外半徑。這里取rbegin為69.282 km，rend為92.376 km，對應(yīng)的高度范圍hbegin為40 km，hend為20 km?？臻g網(wǎng)格步長dr取0.5 m，時間步長取1.0×10-10s，滿足Courant-Friedrich-Levy條件的要求。給出了三個觀測點P1,P2,P3上的計算結(jié)果，它們距離計算區(qū)域內(nèi)半徑分別為50, 500, 1 000 m。計算區(qū)域及觀測點在整個空間的位置和方向示意圖，如圖1所示，圖中粗線部分即為計算區(qū)域。

圖1 計算區(qū)域示意圖Fig.1 Diagram of computational region

3.2 計算結(jié)果

圖2—圖10為P1,P2,P3處電場分量的時域波形圖，并與非自洽處理方法的計算結(jié)果進行了比對。

圖2 觀測點P1的Er時域波形Fig.2 Waveform of Er at test point P1

圖3 觀測點P2的Er時域波形Fig.3 Waveform of Er at test point P2

圖4 觀測點P3的Er時域波形Fig.4 Waveform of Er at test point P3

圖5 觀測點P1的Eθ時域波形Fig.5 Waveform of Eθ at test point P1

圖6 觀測點P2的Eθ時域波形Fig.6 Waveform of Eθ at test point P2

圖7 觀測點P3的Eθ時域波形Fig.7 Waveform of Eθ at test point P3

圖8 觀測點P1上的Eφ時域波形圖Fig.8 Waveform of Eφ at test point P1

圖9 觀測點P2的Eφ時域波形Fig.9 Waveform of Eφ at test point P2

圖10 觀測點P3的Eφ時域波形Fig.10 Waveform of Eφ at test point P3

從圖中可以看出，用自洽處理方法計算出來的電場時域波形和用非自洽處理方法計算出來的電場時域波形整體形狀幾乎一致，到達第1峰值電場強度的時間幾乎相同，差別在于：1)用自洽處理方法計算出來的第1峰值電場強度小于非自洽處理方法的計算結(jié)果，但兩種方法計算的第2、第3峰值電場強度結(jié)果相當(dāng)； 2)用自洽處理方法計算出來的第1峰值電場與第2峰值電場之間的時間間隔小于用非自洽處理方法計算出來的量，會使第1波峰電場與第2峰值電場之間的時間間隔縮短。但從第2波峰往后的波峰和波峰之間的時間間隔，兩種方法的計算結(jié)果幾乎相同。

4 討論

本文采用全電磁粒子模擬方法和一維近似模型，數(shù)值模擬了電磁脈沖場與康普頓電子相互作用產(chǎn)生E1-HEMP的過程。從源區(qū)場的計算結(jié)果可知，考慮自洽效應(yīng)計算出來的第1峰值電場強度會略微減小，到達第2、第3等峰值電場強度的時間變快了，波形形狀發(fā)生了改變，這和Karzas等對自洽作用的估計[3]基本一致。

關(guān)注E1-HEMP，很多時候是關(guān)心它在地面上的時空分布，而要想得到地面上的輻射場，就需要知道到達源區(qū)下邊界處的輻射場[11]。E1-HEMP的源區(qū)高度范圍達20～30 km，而全電磁粒子模擬方法中計算電磁場的FDTD方法對網(wǎng)格步長、時間步長有嚴格的限制，再加上對宏粒子的處理非常費時，所以在有限時間內(nèi)，只能得到少數(shù)幾個觀察點上的電場時域波形，無法得到整個源區(qū)內(nèi)輻射場的時空分布，進而無法得到源區(qū)外輻射場的時域變化。下一步，擬在推遲時間變換下，進一步研究E1-HEMP產(chǎn)生的自洽過程。

[1]王建國, 牛勝利, 張殿輝, 等. 高空核爆炸效應(yīng)參數(shù)手冊[M]. 北京: 原子能出版社, 2010: 145-148. (WANG Jian-guo, NIU Sheng-li, ZHANG Dian-hui, et al. Parameter Handbook of High Altitude Nuclear Detonation Effects [M]. Beijing: Atomic Energy Press, 2010: 145-148.)

[2]KARZAS W J, LATTER R. Electromagnetic radiation from a nuclear explosion in space [J]. Phys Rev, 1962, 126(6): 1 919-1 962.

[3]KARZAS W J, LATTER R. Detection of the electromagnetic radiation from nuclear explosions in space [J]. Phys Rev, 1963, 137(5B): 1 369-1 378.

[4]LONGMIRE C L. On the electromagnetic pulse produced by nuclear explosions [J]. IEEE Trans Antennas Propag, 1978, 26(1): 3-13.

[5]陳雨生. 核爆炸電磁脈沖研究[C]// 核爆電磁脈沖技術(shù)交流會文集, 北京: 國防科工委情報資料研究所, 1980: 1-35. (CHEN Yu-sheng. Research on NEMP technical communion conference corpus on NEMP[C]// Beijing: National Deference Scientific and Engineer Committee, 1980: 1-35.)

[6]WITTER L A, CANAVAN G H, BRAV J E. CHEMP: A code for calculation of high-altitude EMP[R]. AD-783239, Air force Weapons Lab., 1974.

[7]程引會, 李進璽, 馬良, 等. 高空電磁脈沖環(huán)境計算中的自洽方法[J]. 計算物理, 2017, 34(4): 403-408. (CHENG Yin-hui, LI Jin-xi, MA Liang, et al. Self-consistent method in calculation of high-altitude electromagnetic pulse environment [J]. Chinese Journal of Computational Physics, 2017, 34(4): 403-408.)

[8]DAWSON J M. Particle simulation of plasmas [J]. Review of Modern Physics, 1983, 55(2): 403-447.

[9] BIRDSAL C K, LANGDON A B. Plasma Physics via Computer Simulation [M]. New York: Adam Hilger, 1991: 55-79.

[10] BYRD R L. Atmospheric transport of neutrons and gamma rays from a high-altitude nuclear detonation[R]. LA-12962-MS. Los Alamos National Laboratory, 1995.

[11] MENG C. Numerical simulation of HEMP environment [J]. IEEE Trans Electromagn Compat, 2013, 55(3): 440-455.

[12] 王建國. 真空電子器件的粒子模擬方法[J]. 現(xiàn)代應(yīng)用物理, 2013, 4(3): 251-262. (WANG Jian-guo. Particle simulation method of vacuum electronic devices[J]. Modern Applied Physics, 2013, 4(3): 251-262.)

[13] EASTWOOD J W. The virtual particle electromagnetic particle-mesh method[J]. Computer Physics Communications, 1991, 64(2): 252-266.

[14] 付梅艷, 張茂鈺. 用高頻近似模型和入射-出射波模型對早期高空核電磁脈沖場的計算比對[J]. 現(xiàn)代應(yīng)用物理, 2015, 6(2): 107-112. (FU Mei-yan, ZHANG Mao-yu. Comparison of two simplified model of early-time high altitude electromagnetic pulse field equation [J]. Modern Applied Physics, 2015, 6(2): 107-112.)