钡和铯释放的电离层扰动效应对比

朱肖丽 胡耀垓 赵正予 张援农

(武汉大学电子信息学院,电离层实验室,武汉 430079)

碱金属或碱土金属在电离层释放后,迅速在太阳辐射作用下发生光电离,产生正离子和电子,形成人工等离子体云团.本文基于三维双成分流体模型,考虑释放区域水平风场的影响,探讨了钡和铯在电离层释放后的时空演化规律,并对钡和铯的电离层扰动效应进行了对比.模拟结果表明,不考虑中性风场时,生成的等离子体云团逐渐沿磁场被拉伸成椭球形结构,同时,膨胀的等离子体云会推开背景氧离子,在释放中心形成氧离子密度空洞,并在两侧产生两个对称的密度尖峰;水平风场的存在会使得生成的离子云逆风侧的密度梯度变陡,释放物质对背景氧离子的扰动也更大;对比钡与铯的释放结果发现,由于铯的扩散系数较小,钡云的膨胀更为迅速,Ba+云团的覆盖区域更广;而由于光电离率较大,释放相同质量下铯的离子产率更高;此外,Cs+的扫雪机效应比Ba+扫雪机更强,氧离子密度空穴和凸起处的扰动也更大.

1 引 言

20世纪60年代,人类首次发现火箭发射期间对电离层电子密度的扰动现象.此后,美欧等国开展了大量化学物质主动释放扰动电离层的试验[1−4],包括电离层电子密度耗空类物质和电子密度增强类物质,其中,电离层电子密度增强试验通常的释放物质为碱金属或碱土金属,如 Ba,Cs,Li,Na,Ca等.

通过空间物理主动释放试验、实验室研究以及数值模拟研究等对人工等离子体云团释放后的演化特性进行研究已经持续了很多年[5−8].早在1967年,Haerendel等[9]就通过简化的低密度扰动模型对人工等离子体云进行了初步定性探讨.1988 年,Schunk和Szuszczewicz[10]首先利用一维Vlasov-Poisson模型研究了高密度 Ba+,Li+和Ba+-Li+等离子体云向低密度O+背景等离子体的无碰撞扩展,以阐明早期等离子体膨胀的离子特征.随后,Mitchell等[11]利用二维静电模型对等离子体云沿磁场方向和垂直磁场的运动进行了理论和数值研究;为了全面地描述等离子体的膨胀和三维运动,更为精细的三维模型被构建,以研究背景中性风场、电磁场作用、粒子间碰撞效应和惯性作用影响下等离子体云团的膨胀特性[12−15].

化学释放人工扰动电离层的另一类物质为中性气体,从20世纪60年代开始,国外开展了大量中性物质释放形成电离层电子密度“空洞”的空间主动试验[16−18].基于试验观测结果,很多学者也从理论和数值模拟的角度对中性气体释放扰动电离层的机理及一些观测到的效应进行了理论解释[19−22].近年来,国内也开展了中性物质释放人工干扰电离层的相关研究.黄文耿和古士芬[23]、黄勇等[24]、胡耀垓等[25,26]和汪四成等[27]基于二维动力学模型,综合考虑了中性气体热扩散、等离子体双极扩散、以及离子化学反应等过程,数值模拟了多种中性气体释放后对电离层的扰动特性,并比较了不同释放条件下的电离层扰动结果;赵海生等[28]建立了更为精细的模型,考虑了热层风场、释放物初速度及运载器飞行姿态等参量对中性扩散过程的影响,并引入地磁倾角和沿场扩散项,将二维动力学模型扩展到了三维.

相比于中性物质的释放,国内对碱土金属类物质释放的研究较少.Li和Xu[29]研究了在扩散、碰撞、氧化和光电离作用下不同高度中性钡云的自扩散过程.胡耀垓等[30]在中性钡云自扩散模型的基础上,考虑钡原子的氧化和光电离损耗,探讨了释放早期中性钡云形态、亮度分布以及释放区域电子密度分布.谢良海[31]利用美国BATS-R-US代码,基于多成分磁流体模型,模拟了钡的释放效应及其与环境等离子体的相互作用过程.

开展空间化学物质主动释放试验,观测释放物质与背景电离层的相互作用过程,对于相关空间物理问题(如磁重联、等离子体不稳定性、等离子体波等)的研究,高层大气风场和电磁场的测量,实现人工电离层变态,影响短波通信及卫星通信等都具有非常重要的意义[1,4,7,9,29,32].金属钡由于质量轻、电离电位低、易气化以及便于观测等优点,成为空间物理主动试验中最常用的释放物.但是,钡在电离层的电离依赖于太阳光的作用,在黑暗条件下无法进行光电离,而铯作为碱金属族中电离势最低的元素,除了光电离外,铯在无光环境下也能通过自身热电离产生电子,打破了光照条件对释放试验的限制,早期的空间主动试验也常用铯作为释放物[33−35].相比钡而言,铯释放的模拟和试验研究都尚不充分,也鲜见报道.开展钡和铯释放效应的对比研究,比较不同人工等离子体云团的演化规律及其对背景电离层的扰动特性,对于空间主动试验中释放物的选择具有一定的指导意义.

本文基于碱金属在电离层释放的三维双成分流体模型,从粒子的连续性方程、动量方程以及电流平衡方程出发,考虑了背景电磁场力、各种粒子成分间的碰撞及背景水平风场作用,模拟了钡和铯释放后等离子体云团和背景主要粒子的时空演化规律.

2 理论模型

由于释放物质光电离产生的离子和背景电离层离子成分之间具有不同的空间分布,我们将不同离子分开处理.此外,在电离层F层,背景离子成分主要为氧离子,因此这里考虑金属离子和氧离子的双成分模型.忽略由于化学反应引起的背景氧离子的密度变化,并假设离子云的产生仅来源于光电离,释放区域等离子体满足的一般方程如下:

其中下标1,2,e和s分别表示氧离子、金属离子、电子和释放的中性物质.n和u分别表示各粒子的数密度和漂移速度.(1)式和(2)式是离子的连续性方程,σ 是离子云的光电离产生率.假设光电离产生率为离子的净生成速率,即忽略电荷交换和复合效应.在电离层F区,准中性条件是满足的(方程 (3)).离子的动量方程由 (4)式给出,其中,nj,mj,uj,ej和pj分别表示离子 j的数密度、质量、漂移速度、电荷量和压强;Dj/Dt=∂/∂t+uj·∇表示迁移导数;E和 B是电场和磁场;δ Mj/δt 代表动量变化;∇ pj=∇kTjnj为压强梯度项;J 代表总电流密度.

氧离子和金属离子的动量变化分别为:

νjα是粒子间的碰撞频率(α,s表示带电粒子和中性粒子),us是中性原子的速度.除了由于碰撞带来的动量变化,离子2还应包括源项引起的动量变化.由于地磁场的束缚作用,在垂直于磁场的方向上带电粒子之间的速度差很小,与之有关的碰撞项可以忽略[14],同时忽略电子惯性,将离子惯性项作为小的修正项,并假设地磁场恒定,求解动量方程可以得到各粒子漂移速度的表达式[14]:

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

电场由公式 E=E0-∇φ 计算,其中 E0是背景恒定电场,φ是扰动静电势,定义有效电势Φ=φ-(kTe/e)ln(ne),根据电荷守恒,可以得到关于有效电势的方程如下:

其中σ⊥为 Pedersen 电导 率,σ//为平行电导率,Q为总的电荷量[14].

释放中性云团粒子的密度和速度由以下方程给出[27]:

其中 ns是中性云的数密度;us是中性云团速度;σ表示光电离率;kc是氧化反应速率(通常在电离层、磁层高度处 σ ≫kcnO2,可不考虑氧化损耗);nO2是氧气数密度;D=kT/(mνD) 是扩散系数,νD为与扩散相关的碰撞频率;(x0,y0,z0) 为云团初始释放中心;ux,uy和uz分别是云团中心在三个方向上的速度分量;r0是云团初始半径;N0是释放粒子总数;ε 表示形状因子.扩散系数D可由以下公式计算[36]:

其中,nO为氧原子数密度,为约化质量,ms,mO,rs和rO分别为金属原子和氧原子的质量和原子半径.在电离层F区,中性云会在背景粒子的碰撞作用影响下做减速运动,其速度随时间的变化可近似表示为

这里 u0为云团初始释放速度,un为中性风场速度,νs为云团减速运动的阻尼系数,可根据弹性碰撞理论求得[37]:

考虑背景风场作用,对于有初速的释放情形,叠加的移动量如下:

对上述(1)—(17)式进行数值求解,即可得到任意时刻释放区域粒子的空间分布.

3 数值模拟结果及讨论

3.1 仿真算法设计

利用有限差分方法求解连续性方程(1)和(2),方程(12)是椭圆方程,利用交替方向隐式算法求解.时间步长取 ∇ t=0.01 s,空间步长∇x=∇y=∇z=1 km,基于MATLAB软件实现了人工等离子体云演化的模拟程序,数值仿真算法流程如图1所示.

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

背景大气密度、电离层粒子数密度、温度、磁场强度等初始条件可由MSIS-E-90大气模型、电离层IRI-2016模型和磁场模型DGRF/IGRF 1945-2020获得.所有边界条件均采用等值外推形式.主要仿真参数如表1所列.

表1 主要仿真参数表Table 1.The main simulation parameters.

3.2 数值模拟结果及讨论

为探讨金属原子在电离层释放后产生的扰动特性和释放区域主要粒子的时空演化规律,分别模拟了10 kg钡和铯在300 km高度处的释放,释放云团的初始特征半径均取 2 km,计算域为 X=Y=Z=[–25 25]km.所有的剖面图都是经过释放中心的平面,图中的等值图代表粒子数密度分布,带电粒子的速度用矢量场表示(图中黑色箭头).

3.2.1 无中性风场时

图2和图3描述了 300 km高度释放 10 kg Ba 原子,释放后 5,30和200 s钡离子及背景氧离子的数密度分布.不考虑背景风速时,钡离子云的早期密度分布为球对称的,由于磁场的存在,钡离子云在垂直磁场方向上的膨胀受到束缚(图2),而在沿磁场方向,由于钡离子云在密度梯度作用下的运动不受限制,离子云团逐渐沿着磁场方向被拉伸,逐渐变成椭球状结构 (图3).由于碰撞作用,钡离子和氧离子的动量相互耦合,钡离子平行于磁场的动能传递给氧离子,促使氧离子沿着磁场向两侧运动,不断将中心处的氧离子往外传送,形成氧离子密度“空穴”;O+到达两侧后,又受到背景热压梯度的作用,阻止了其向外的膨胀,最终造成氧离子在两侧的堆积,产生两个密度尖峰,即所谓的“扫雪机效应”[14],氧离子空穴处的形态分布与钡离子云团形态是一致的.5 s时,空穴区的氧离子密度比背景降低了27.6%,密度凸起处比背景高14.3%,钡离子的数密度峰值达到 1.332×107cm–3.

图2 无中性风场时,300 km 高度释放 10 kg 钡后钡离子和氧离子的离子数密度分布 (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 无中性风场时,300 km 高度释放 10 kg 钡后钡离子和氧离子的粒子数密度分布 (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 考虑中性风场的情况

图4给出了x方向存在1 km/s的中性风时,300 km高度释放10 kg钡的模拟结果(x-z平面),释放中心位于 (–15,0,0) km 处,由于带电粒子的密度分布在空间上是柱对称的,这里不再给出x-y平面的图像.与无风场的结果相比,中性风场的存在使得释放的中性原子有了一个水平方向的运动,生成的等离子体云团不再是对称的椭球状结构,而是一个逆风侧的密度梯度“较陡”,顺风侧的密度梯度较为“平缓”的不对称结构.比较图3和图4可以发现,中性风场存在时,释放相同质量钡原子在早期对背景氧离子的密度扰动更大,由于中性云的移动,扰动区域也有所增加.5 s 时,空穴区的氧离子密度比背景降低30.7%,密度凸起处比背景高17.7%,由于分布区域变广,钡离子的数密度峰值略微降低,为 1.223×107cm–3.

图4 存在 x 方向大小为 1 km/s的中性风时,300 km 高度释放 10 kg 钡后钡离子和氧离子的粒子数密度分布 (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中右侧的绿色球表示中性钡原子,蓝色部分代表钡离子,虽然钡离子的运动在垂直磁场方向上受到磁场的束缚,但中性钡原子不受磁场力的约束,可以自由地穿过磁场线,因此钡离子云团与中性云团会慢慢分离,如图5 所示,30 s时离子云和中性云的分离已经较为明显.中性钡原子在背景中性风的作用下向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.

的离子结构(蓝色部分),较早产生的钡离子已经沿磁场方向被拉伸(蓝色球左侧沿z方向拉伸的区域).

3.2.3 铯释放结果

10 kg Cs原子释放后 5,30和200 s铯离子及背景氧离子密度分布的演化过程如图6所示,除释放物质外,其他释放参数与 3.2.2节中一致.虽然铯由于电离势较低容易发生热电离,但在300 km高度及有日照的条件下,其自身的热电离与光电离相比是微不足道的,因此这里没有对铯的热电离过程进行探讨.定性地说,Cs+和 Ba+的膨胀特性以及对背景O+的扰动效应是类似的,虽然铯比钡原子的质量小,但铯原子的极化率较钡原子要大,因而铯的扩散系数较小,钡云的膨胀更为迅速,钡离子云团的覆盖区域更广;同时,由于铯的光电离率较大,释放相同质量下铯的离子产率更高,5 s时铯离子云的峰值数密度达到 2.248×107cm–3,是相同条件下产生的钡离子云密度的近两倍;此外,Cs+-O+的碰撞频率较 Ba+-O+更大,Cs+的扫雪机效应比Ba+扫雪机更强,氧离子密度空穴和凸起更大,30 s时钡和铯释放产生的背景氧离子的最大扰动分别为74.2%和75.1%.

图7给出了背景中性风场存在下10 kg铯和钡在300 km高度释放后,生成的等离子体云团密度的最大值以及背景氧离子密度凸起最高点的值随时间的变化情况对比.释放初期,伴随着持续的光电离,钡离子和铯离子数密度持续增加,之后由于扩散作用逐渐下降,铯离子的数密度最大值比钡离子大,但随时间的变化趋势是一致的;随着等离子体云团密度的增加,增强的密度梯度使得越来越多的氧离子被推开,两侧密度凸起处的值不断增大,随后,密度梯度开始下降,氧离子扰动也逐渐回复.

图6 存在 x 方向大小为 1 km/s 的中性风时,300 km 高度释放 10 kg 铯的粒子数密度分布 (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 结 论

采用三维双成分的流体模型研究了电离层人工等离子体云的演化,分别给出了钡在有无中性风场时的模拟结果,并对存在背景风场时钡和铯释放的电离层扰动效应进行了对比.基于已有的模拟结果,得到以下结论:1)钡和铯在电离层释放后迅速向四周膨胀,伴随着太阳辐射作用下的光电离,生成的等离子体云在磁场的束缚下逐渐由内向外呈椭球形结构,同时,由于离子间的动量耦合,膨胀的等离子体云会推开背景氧离子,在释放中心形成氧离子密度空洞,并在两侧产生两个氧离子密度尖峰,背景O+空洞的形状类似于生成的人工等离子体云团的形状;2)在不存在中性风的情况下,等离子体云团以平行磁场的运动为主,考虑背景水平风场时,云和背景受扰区域除了沿磁场方向的运动外,还会有沿着风速方向的移动,逆风侧的等离子体云团密度梯度变陡,而顺风侧密度梯度较为平缓,虽然离子云在横向上的运动会受到磁场的约束,但中性云团可以自由地穿过磁场线,因此离子云与中性云会慢慢分离;3)对比钡和铯的模拟结果发现,定性地说,Cs+和 Ba+的膨胀特性以及对背景O+的扰动效应是类似的,由于铯的扩散系数较小,钡云的膨胀更为迅速,钡离子云团的覆盖区域更广,同时,由于铯的光电离率较大,释放相同质量下铯的离子产率更高,此外,Cs+-O+的碰撞频率较Ba+-O+更大,Cs+的扫雪机效应比Ba+扫雪机更强,造成的氧离子密度扰动也更大.

本文仅利用简化的模型进行了研究,在下一步的工作中,可以综合考虑地磁倾角、磁偏角、背景电离层密度梯度、温度梯度等因素对模型加以完善;此外,对不同释放条件以及不同地球物理条件包括不同太阳和地磁活动、季节以及地域等条件下人工等离子体云团的形成和演化也可以进行更加深入、系统的研究.