付梅艷, 張茂鈺

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

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

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

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

1 理論模型

1.1 麥克斯韋方程組

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

(1)

(2)

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

1.2 康普頓電流

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

1.2.1 康普頓電子的產生

高空核爆炸時，產生大量的γ光子。在分析γ光子在大氣中輻射傳輸并產生康普頓電子的過程時，采取2個假設近似：1)每一個γ光子發(fā)生康普頓散射只散射出一個康普頓電子。不再考慮散射γ光子。2)γ光子與周圍大氣發(fā)生康普頓散射作用的概率是γ光子在當地的自由程的倒數，即γ光子在一個自由程內均勻產生出康普頓電子。

爆點處，γ光子的產生率為[5]

(3)

(4)

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

(5)

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

(6)

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

1.2.2 康普頓電子的運動

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

(7)

(8)

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

1.2.3 康普頓電子的消失

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

(9)

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

1.3 電導率

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

(10)

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

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

(11)

2 計算模型

2.1 宏粒子模型

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

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

2.2 電磁場方程的計算

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

(12)

(13)

(14)

(15)

(16)

(17)

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

左邊界條件為

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

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

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

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

3 算例

3.1 計算條件

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

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

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

3.2 計算結果

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

圖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峰值電場強度小于非自洽處理方法的計算結果，但兩種方法計算的第2、第3峰值電場強度結果相當； 2)用自洽處理方法計算出來的第1峰值電場與第2峰值電場之間的時間間隔小于用非自洽處理方法計算出來的量，會使第1波峰電場與第2峰值電場之間的時間間隔縮短。但從第2波峰往后的波峰和波峰之間的時間間隔，兩種方法的計算結果幾乎相同。

4 討論

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

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

[1]王建國, 牛勝利, 張殿輝, 等. 高空核爆炸效應參數手冊[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]// 核爆電磁脈沖技術交流會文集, 北京: 國防科工委情報資料研究所, 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]. 現代應用物理, 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]. 現代應用物理, 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.)