朱肖麗 胡耀垓 趙正予 張援農(nóng)

(武漢大學(xué)電子信息學(xué)院,電離層實驗室,武漢 430079)

堿金屬或堿土金屬在電離層釋放后,迅速在太陽輻射作用下發(fā)生光電離,產(chǎn)生正離子和電子,形成人工等離子體云團.本文基于三維雙成分流體模型,考慮釋放區(qū)域水平風(fēng)場的影響,探討了鋇和銫在電離層釋放后的時空演化規(guī)律,并對鋇和銫的電離層擾動效應(yīng)進行了對比.模擬結(jié)果表明,不考慮中性風(fēng)場時,生成的等離子體云團逐漸沿磁場被拉伸成橢球形結(jié)構(gòu),同時,膨脹的等離子體云會推開背景氧離子,在釋放中心形成氧離子密度空洞,并在兩側(cè)產(chǎn)生兩個對稱的密度尖峰;水平風(fēng)場的存在會使得生成的離子云逆風(fēng)側(cè)的密度梯度變陡,釋放物質(zhì)對背景氧離子的擾動也更大;對比鋇與銫的釋放結(jié)果發(fā)現(xiàn),由于銫的擴散系數(shù)較小,鋇云的膨脹更為迅速,Ba+云團的覆蓋區(qū)域更廣;而由于光電離率較大,釋放相同質(zhì)量下銫的離子產(chǎn)率更高;此外,Cs+的掃雪機效應(yīng)比Ba+掃雪機更強,氧離子密度空穴和凸起處的擾動也更大.

1 引 言

20世紀(jì)60年代,人類首次發(fā)現(xiàn)火箭發(fā)射期間對電離層電子密度的擾動現(xiàn)象.此后,美歐等國開展了大量化學(xué)物質(zhì)主動釋放擾動電離層的試驗[1?4],包括電離層電子密度耗空類物質(zhì)和電子密度增強類物質(zhì),其中,電離層電子密度增強試驗通常的釋放物質(zhì)為堿金屬或堿土金屬,如 Ba,Cs,Li,Na,Ca等.

通過空間物理主動釋放試驗、實驗室研究以及數(shù)值模擬研究等對人工等離子體云團釋放后的演化特性進行研究已經(jīng)持續(xù)了很多年[5?8].早在1967年,Haerendel等[9]就通過簡化的低密度擾動模型對人工等離子體云進行了初步定性探討.1988 年,Schunk和Szuszczewicz[10]首先利用一維Vlasov-Poisson模型研究了高密度 Ba+,Li+和Ba+-Li+等離子體云向低密度O+背景等離子體的無碰撞擴展,以闡明早期等離子體膨脹的離子特征.隨后,Mitchell等[11]利用二維靜電模型對等離子體云沿磁場方向和垂直磁場的運動進行了理論和數(shù)值研究;為了全面地描述等離子體的膨脹和三維運動,更為精細的三維模型被構(gòu)建,以研究背景中性風(fēng)場、電磁場作用、粒子間碰撞效應(yīng)和慣性作用影響下等離子體云團的膨脹特性[12?15].

化學(xué)釋放人工擾動電離層的另一類物質(zhì)為中性氣體,從20世紀(jì)60年代開始,國外開展了大量中性物質(zhì)釋放形成電離層電子密度“空洞”的空間主動試驗[16?18].基于試驗觀測結(jié)果,很多學(xué)者也從理論和數(shù)值模擬的角度對中性氣體釋放擾動電離層的機理及一些觀測到的效應(yīng)進行了理論解釋[19?22].近年來,國內(nèi)也開展了中性物質(zhì)釋放人工干擾電離層的相關(guān)研究.黃文耿和古士芬[23]、黃勇等[24]、胡耀垓等[25,26]和汪四成等[27]基于二維動力學(xué)模型,綜合考慮了中性氣體熱擴散、等離子體雙極擴散、以及離子化學(xué)反應(yīng)等過程,數(shù)值模擬了多種中性氣體釋放后對電離層的擾動特性,并比較了不同釋放條件下的電離層擾動結(jié)果;趙海生等[28]建立了更為精細的模型,考慮了熱層風(fēng)場、釋放物初速度及運載器飛行姿態(tài)等參量對中性擴散過程的影響,并引入地磁傾角和沿場擴散項,將二維動力學(xué)模型擴展到了三維.

相比于中性物質(zhì)的釋放,國內(nèi)對堿土金屬類物質(zhì)釋放的研究較少.Li和Xu[29]研究了在擴散、碰撞、氧化和光電離作用下不同高度中性鋇云的自擴散過程.胡耀垓等[30]在中性鋇云自擴散模型的基礎(chǔ)上,考慮鋇原子的氧化和光電離損耗,探討了釋放早期中性鋇云形態(tài)、亮度分布以及釋放區(qū)域電子密度分布.謝良海[31]利用美國BATS-R-US代碼,基于多成分磁流體模型,模擬了鋇的釋放效應(yīng)及其與環(huán)境等離子體的相互作用過程.

開展空間化學(xué)物質(zhì)主動釋放試驗,觀測釋放物質(zhì)與背景電離層的相互作用過程,對于相關(guān)空間物理問題(如磁重聯(lián)、等離子體不穩(wěn)定性、等離子體波等)的研究,高層大氣風(fēng)場和電磁場的測量,實現(xiàn)人工電離層變態(tài),影響短波通信及衛(wèi)星通信等都具有非常重要的意義[1,4,7,9,29,32].金屬鋇由于質(zhì)量輕、電離電位低、易氣化以及便于觀測等優(yōu)點,成為空間物理主動試驗中最常用的釋放物.但是,鋇在電離層的電離依賴于太陽光的作用,在黑暗條件下無法進行光電離,而銫作為堿金屬族中電離勢最低的元素,除了光電離外,銫在無光環(huán)境下也能通過自身熱電離產(chǎn)生電子,打破了光照條件對釋放試驗的限制,早期的空間主動試驗也常用銫作為釋放物[33?35].相比鋇而言,銫釋放的模擬和試驗研究都尚不充分,也鮮見報道.開展鋇和銫釋放效應(yīng)的對比研究,比較不同人工等離子體云團的演化規(guī)律及其對背景電離層的擾動特性,對于空間主動試驗中釋放物的選擇具有一定的指導(dǎo)意義.

本文基于堿金屬在電離層釋放的三維雙成分流體模型,從粒子的連續(xù)性方程、動量方程以及電流平衡方程出發(fā),考慮了背景電磁場力、各種粒子成分間的碰撞及背景水平風(fēng)場作用,模擬了鋇和銫釋放后等離子體云團和背景主要粒子的時空演化規(guī)律.

2 理論模型

由于釋放物質(zhì)光電離產(chǎn)生的離子和背景電離層離子成分之間具有不同的空間分布,我們將不同離子分開處理.此外,在電離層F層,背景離子成分主要為氧離子,因此這里考慮金屬離子和氧離子的雙成分模型.忽略由于化學(xué)反應(yīng)引起的背景氧離子的密度變化,并假設(shè)離子云的產(chǎn)生僅來源于光電離,釋放區(qū)域等離子體滿足的一般方程如下:

其中下標(biāo)1,2,e和s分別表示氧離子、金屬離子、電子和釋放的中性物質(zhì).n和u分別表示各粒子的數(shù)密度和漂移速度.(1)式和(2)式是離子的連續(xù)性方程,σ 是離子云的光電離產(chǎn)生率.假設(shè)光電離產(chǎn)生率為離子的凈生成速率,即忽略電荷交換和復(fù)合效應(yīng).在電離層F區(qū),準(zhǔn)中性條件是滿足的(方程 (3)).離子的動量方程由 (4)式給出,其中,nj,mj,uj,ej和pj分別表示離子 j的數(shù)密度、質(zhì)量、漂移速度、電荷量和壓強;Dj/Dt=?/?t+uj·?表示遷移導(dǎo)數(shù);E和 B是電場和磁場;δ Mj/δt 代表動量變化;? pj=?kTjnj為壓強梯度項;J 代表總電流密度.

氧離子和金屬離子的動量變化分別為:

νjα是粒子間的碰撞頻率(α,s表示帶電粒子和中性粒子),us是中性原子的速度.除了由于碰撞帶來的動量變化,離子2還應(yīng)包括源項引起的動量變化.由于地磁場的束縛作用,在垂直于磁場的方向上帶電粒子之間的速度差很小,與之有關(guān)的碰撞項可以忽略[14],同時忽略電子慣性,將離子慣性項作為小的修正項,并假設(shè)地磁場恒定,求解動量方程可以得到各粒子漂移速度的表達式[14]:

其中 z=B/|B| 表示沿磁場方向的單位向量.Ωj=ejB/mj是回旋頻率.aj和sj的表達式如下:

電場由公式 E=E0-?φ 計算,其中 E0是背景恒定電場,φ是擾動靜電勢,定義有效電勢Φ=φ-(kTe/e)ln(ne),根據(jù)電荷守恒,可以得到關(guān)于有效電勢的方程如下:

其中σ⊥為 Pedersen 電導(dǎo) 率,σ//為平行電導(dǎo)率,Q為總的電荷量[14].

釋放中性云團粒子的密度和速度由以下方程給出[27]:

其中 ns是中性云的數(shù)密度;us是中性云團速度;σ表示光電離率;kc是氧化反應(yīng)速率(通常在電離層、磁層高度處 σ ?kcnO2,可不考慮氧化損耗);nO2是氧氣數(shù)密度;D=kT/(mνD) 是擴散系數(shù),νD為與擴散相關(guān)的碰撞頻率;(x0,y0,z0) 為云團初始釋放中心;ux,uy和uz分別是云團中心在三個方向上的速度分量;r0是云團初始半徑;N0是釋放粒子總數(shù);ε 表示形狀因子.擴散系數(shù)D可由以下公式計算[36]:

其中,nO為氧原子數(shù)密度,為約化質(zhì)量,ms,mO,rs和rO分別為金屬原子和氧原子的質(zhì)量和原子半徑.在電離層F區(qū),中性云會在背景粒子的碰撞作用影響下做減速運動,其速度隨時間的變化可近似表示為

這里 u0為云團初始釋放速度,un為中性風(fēng)場速度,νs為云團減速運動的阻尼系數(shù),可根據(jù)彈性碰撞理論求得[37]:

考慮背景風(fēng)場作用,對于有初速的釋放情形,疊加的移動量如下:

對上述(1)—(17)式進行數(shù)值求解,即可得到任意時刻釋放區(qū)域粒子的空間分布.

3 數(shù)值模擬結(jié)果及討論

3.1 仿真算法設(shè)計

利用有限差分方法求解連續(xù)性方程(1)和(2),方程(12)是橢圓方程,利用交替方向隱式算法求解.時間步長取 ? t=0.01 s,空間步長?x=?y=?z=1 km,基于MATLAB軟件實現(xiàn)了人工等離子體云演化的模擬程序,數(shù)值仿真算法流程如圖1所示.

圖1 仿真算法流程圖Fig.1.Flow chart of simulation algorithm.

背景大氣密度、電離層粒子數(shù)密度、溫度、磁場強度等初始條件可由MSIS-E-90大氣模型、電離層IRI-2016模型和磁場模型DGRF/IGRF 1945-2020獲得.所有邊界條件均采用等值外推形式.主要仿真參數(shù)如表1所列.

表1 主要仿真參數(shù)表Table 1.The main simulation parameters.

3.2 數(shù)值模擬結(jié)果及討論

為探討金屬原子在電離層釋放后產(chǎn)生的擾動特性和釋放區(qū)域主要粒子的時空演化規(guī)律,分別模擬了10 kg鋇和銫在300 km高度處的釋放,釋放云團的初始特征半徑均取 2 km,計算域為 X=Y=Z=[–25 25]km.所有的剖面圖都是經(jīng)過釋放中心的平面,圖中的等值圖代表粒子數(shù)密度分布,帶電粒子的速度用矢量場表示(圖中黑色箭頭).

3.2.1 無中性風(fēng)場時

圖2和圖3描述了 300 km高度釋放 10 kg Ba 原子,釋放后 5,30和200 s鋇離子及背景氧離子的數(shù)密度分布.不考慮背景風(fēng)速時,鋇離子云的早期密度分布為球?qū)ΨQ的,由于磁場的存在,鋇離子云在垂直磁場方向上的膨脹受到束縛(圖2),而在沿磁場方向,由于鋇離子云在密度梯度作用下的運動不受限制,離子云團逐漸沿著磁場方向被拉伸,逐漸變成橢球狀結(jié)構(gòu) (圖3).由于碰撞作用,鋇離子和氧離子的動量相互耦合,鋇離子平行于磁場的動能傳遞給氧離子,促使氧離子沿著磁場向兩側(cè)運動,不斷將中心處的氧離子往外傳送,形成氧離子密度“空穴”;O+到達兩側(cè)后,又受到背景熱壓梯度的作用,阻止了其向外的膨脹,最終造成氧離子在兩側(cè)的堆積,產(chǎn)生兩個密度尖峰,即所謂的“掃雪機效應(yīng)”[14],氧離子空穴處的形態(tài)分布與鋇離子云團形態(tài)是一致的.5 s時,空穴區(qū)的氧離子密度比背景降低了27.6%,密度凸起處比背景高14.3%,鋇離子的數(shù)密度峰值達到 1.332×107cm–3.

圖2 無中性風(fēng)場時,300 km 高度釋放 10 kg 鋇后鋇離子和氧離子的離子數(shù)密度分布 (x-y 平面) (a) Ba+,t=5 s;(b) Ba+,t=30 s;(c) Ba+,t=200 s;(d) O+,t=5 s;(e) O+,t=30 s;(f) O+,t=200 sFig.2.Density distribution of Ba+ and O+ (in x-y plane) after 10 kg barium released at 300 km while no neutral wind is considered:(a) Ba+,t=5 s;(b) Ba+,t=30 s;(c) Ba+,t=200 s;(d) O+,t=5 s;(e) O+,t=30 s;(f) O+,t=200 s.

圖3 無中性風(fēng)場時,300 km 高度釋放 10 kg 鋇后鋇離子和氧離子的粒子數(shù)密度分布 (x-z 平面) (a) O+,t=5 s;(b) O+,t=30 s;(c) O+,t=200 s;(d) Ba+,t=5 s;(e) Ba+,t=30 s;(f) Ba+,t=200 sFig.3.Density distribution of Ba+ and O+ (in x-z plane) after 10 kg barium released at 300 km while no neutral wind is considered:(a) O+,t=5 s;(b) O+,t=30 s;(c) O+,t=200 s;(d) Ba+,t=5 s;(e) Ba+,t=30 s;(f) Ba+,t=200 s.

3.2.2 考慮中性風(fēng)場的情況

圖4給出了x方向存在1 km/s的中性風(fēng)時,300 km高度釋放10 kg鋇的模擬結(jié)果(x-z平面),釋放中心位于 (–15,0,0) km 處,由于帶電粒子的密度分布在空間上是柱對稱的,這里不再給出x-y平面的圖像.與無風(fēng)場的結(jié)果相比,中性風(fēng)場的存在使得釋放的中性原子有了一個水平方向的運動,生成的等離子體云團不再是對稱的橢球狀結(jié)構(gòu),而是一個逆風(fēng)側(cè)的密度梯度“較陡”,順風(fēng)側(cè)的密度梯度較為“平緩”的不對稱結(jié)構(gòu).比較圖3和圖4可以發(fā)現(xiàn),中性風(fēng)場存在時,釋放相同質(zhì)量鋇原子在早期對背景氧離子的密度擾動更大,由于中性云的移動,擾動區(qū)域也有所增加.5 s 時,空穴區(qū)的氧離子密度比背景降低30.7%,密度凸起處比背景高17.7%,由于分布區(qū)域變廣,鋇離子的數(shù)密度峰值略微降低,為 1.223×107cm–3.

圖4 存在 x 方向大小為 1 km/s的中性風(fēng)時,300 km 高度釋放 10 kg 鋇后鋇離子和氧離子的粒子數(shù)密度分布 (x-z 平面) (a) O+,t=5 s;(b) O+,t=30 s;(c) O+,t=200 s;(d) Ba+,t=5 s;(e) Ba+,t=30 s;(f) Ba+,t=200 sFig.4.Density distribution of Ba+ and O+ (in x-z plane) after 10 kg barium released at 300 km with a neutral wind of 1 km/s in the x direction:(a) O+,t=5 s;(b) O+,t=30 s;(c) O+,t=200 s;(d) Ba+,t=5 s;(e) Ba+,t=30 s;(f) Ba+,t=200 s.

圖5中右側(cè)的綠色球表示中性鋇原子,藍色部分代表鋇離子,雖然鋇離子的運動在垂直磁場方向上受到磁場的束縛,但中性鋇原子不受磁場力的約束,可以自由地穿過磁場線,因此鋇離子云團與中性云團會慢慢分離,如圖5 所示,30 s時離子云和中性云的分離已經(jīng)較為明顯.中性鋇原子在背景中性風(fēng)的作用下向x軸正方向移動,由于移動過程中持續(xù)的光電離作用,鋇中性云運動后方形成了拉長

圖5 鋇中性云團 (綠色)和離子云團 (藍色)在釋放后30 s時的三維分布示意圖Fig.5.Three-dimensional density distribution of barium neutral cloud (green sphere) and ion cloud (blue sphere) at 30 s after release.

的離子結(jié)構(gòu)(藍色部分),較早產(chǎn)生的鋇離子已經(jīng)沿磁場方向被拉伸(藍色球左側(cè)沿z方向拉伸的區(qū)域).

3.2.3 銫釋放結(jié)果

10 kg Cs原子釋放后 5,30和200 s銫離子及背景氧離子密度分布的演化過程如圖6所示,除釋放物質(zhì)外,其他釋放參數(shù)與 3.2.2節(jié)中一致.雖然銫由于電離勢較低容易發(fā)生熱電離,但在300 km高度及有日照的條件下,其自身的熱電離與光電離相比是微不足道的,因此這里沒有對銫的熱電離過程進行探討.定性地說,Cs+和 Ba+的膨脹特性以及對背景O+的擾動效應(yīng)是類似的,雖然銫比鋇原子的質(zhì)量小,但銫原子的極化率較鋇原子要大,因而銫的擴散系數(shù)較小,鋇云的膨脹更為迅速,鋇離子云團的覆蓋區(qū)域更廣;同時,由于銫的光電離率較大,釋放相同質(zhì)量下銫的離子產(chǎn)率更高,5 s時銫離子云的峰值數(shù)密度達到 2.248×107cm–3,是相同條件下產(chǎn)生的鋇離子云密度的近兩倍;此外,Cs+-O+的碰撞頻率較 Ba+-O+更大,Cs+的掃雪機效應(yīng)比Ba+掃雪機更強,氧離子密度空穴和凸起更大,30 s時鋇和銫釋放產(chǎn)生的背景氧離子的最大擾動分別為74.2%和75.1%.

圖7給出了背景中性風(fēng)場存在下10 kg銫和鋇在300 km高度釋放后,生成的等離子體云團密度的最大值以及背景氧離子密度凸起最高點的值隨時間的變化情況對比.釋放初期,伴隨著持續(xù)的光電離,鋇離子和銫離子數(shù)密度持續(xù)增加,之后由于擴散作用逐漸下降,銫離子的數(shù)密度最大值比鋇離子大,但隨時間的變化趨勢是一致的;隨著等離子體云團密度的增加,增強的密度梯度使得越來越多的氧離子被推開,兩側(cè)密度凸起處的值不斷增大,隨后,密度梯度開始下降,氧離子擾動也逐漸回復(fù).

圖6 存在 x 方向大小為 1 km/s 的中性風(fēng)時,300 km 高度釋放 10 kg 銫的粒子數(shù)密度分布 (x-z 平面) (a) O+,t=5 s;(b) O+,t=30 s;(c) O+,t=200 s;(d) Cs+,t=5 s;(e) Cs+,t=30 s;(f) Cs+,t=200 sFig.6.Density distribution of Cs+ and O+ (in x-z plane) after 10 kg cesium released at 300 km with a neutral wind of 1 km/s in the x direction:(a) O+,t=5 s;(b) O+,t=30 s;(c) O+,t=200 s;(d) Cs+,t=5 s;(e) Cs+,t=30 s;(f) Cs+,t=200 s.

圖7 生成的等離子體云團的密度最大值(a)和背景氧離子的最大擾動值(b)隨時間的變化Fig.7.The maximum density of artificial plasma cloud (a) and the maximum disturbance of background oxygen ion (a) versus time.

4 結(jié) 論

采用三維雙成分的流體模型研究了電離層人工等離子體云的演化,分別給出了鋇在有無中性風(fēng)場時的模擬結(jié)果,并對存在背景風(fēng)場時鋇和銫釋放的電離層擾動效應(yīng)進行了對比.基于已有的模擬結(jié)果,得到以下結(jié)論:1)鋇和銫在電離層釋放后迅速向四周膨脹,伴隨著太陽輻射作用下的光電離,生成的等離子體云在磁場的束縛下逐漸由內(nèi)向外呈橢球形結(jié)構(gòu),同時,由于離子間的動量耦合,膨脹的等離子體云會推開背景氧離子,在釋放中心形成氧離子密度空洞,并在兩側(cè)產(chǎn)生兩個氧離子密度尖峰,背景O+空洞的形狀類似于生成的人工等離子體云團的形狀;2)在不存在中性風(fēng)的情況下,等離子體云團以平行磁場的運動為主,考慮背景水平風(fēng)場時,云和背景受擾區(qū)域除了沿磁場方向的運動外,還會有沿著風(fēng)速方向的移動,逆風(fēng)側(cè)的等離子體云團密度梯度變陡,而順風(fēng)側(cè)密度梯度較為平緩,雖然離子云在橫向上的運動會受到磁場的約束,但中性云團可以自由地穿過磁場線,因此離子云與中性云會慢慢分離;3)對比鋇和銫的模擬結(jié)果發(fā)現(xiàn),定性地說,Cs+和 Ba+的膨脹特性以及對背景O+的擾動效應(yīng)是類似的,由于銫的擴散系數(shù)較小,鋇云的膨脹更為迅速,鋇離子云團的覆蓋區(qū)域更廣,同時,由于銫的光電離率較大,釋放相同質(zhì)量下銫的離子產(chǎn)率更高,此外,Cs+-O+的碰撞頻率較Ba+-O+更大,Cs+的掃雪機效應(yīng)比Ba+掃雪機更強,造成的氧離子密度擾動也更大.

本文僅利用簡化的模型進行了研究,在下一步的工作中,可以綜合考慮地磁傾角、磁偏角、背景電離層密度梯度、溫度梯度等因素對模型加以完善;此外,對不同釋放條件以及不同地球物理條件包括不同太陽和地磁活動、季節(jié)以及地域等條件下人工等離子體云團的形成和演化也可以進行更加深入、系統(tǒng)的研究.